Hacking with GeoScript
You may be surprised to hear that the GeoScript project is nearly 3 years old. In this time it’s accumulated a good deal of functionality. I find myself using the library daily to perform a variety of tasks and I thought sharing some of the most common things I do with GeoScript might be helpful to you.
Recently I worked on a project that involved analyzing and improving GeoServer rendering performance with Oracle. The first step was to gather a baseline measurement of “acceptable performance”. In order to measure respective rendering time between the two databases with the same data our plan was to port all the data to PostGIS. As requested the client sent a dump of their Oracle database and we got to work on converting it:
from geoscript.workspace import Oracle, PostGIS # connect to oracle server ora = Oracle('oradb', host='10.0.1.3', user='jdeolive', 'passwd=...') # connect to local postgis pg = PostGIS('postgis') # grab all the layers from oracle and dump them into postgis for layer_name in ora.layers(): layer = ora[layer_name] pg.add(layer) ora.close() pg.close()
Voila! After some GeoScript magic we had a PostGIS database with all the client data.
One of the most common uses I find for GeoScript is to quickly prototype and create styles. Users of GeoServer know the pain of authoring styles directly in SLD; even with visual tools common tasks like creating “road casings” are often tedious. Here’s an example of some OSM styling that GeoScript made easy.
First I grabbed some OSM Shapefile data:
from geoscript.layer import Shapefile from geoscript.geom import Bounds from geoscript.render import draw # load the data shp = Shapefile('alberta_highway.shp') # the area we care about area = Bounds(-12842844.95, 6636065.22, -12841546.23, 6636833.98, 'epsg:900913') # draw it draw(shp, bounds=area, bgcolor='#efebe3')
Then I created a quick casing style with a label:
from geoscript.style import * # create the casing style = Stroke('black', 12, cap='round') style += Stroke('white', 11, cap='round').zindex(1) # add a label, following the lines label = Label('NAME', 'DejaVuSans').line(group=True, follow=True) style += label draw(shp, style, bounds=area, bgcolor='#efebe3')
Still a bit cluttered. Time to filter it down to just what I want to see.
# filter to roads we care about style = style.where("TYPE IN ('residential', 'tertiary')") draw(shp, style, bounds=area, bgcolor='#efebe3')
Better. Time to export the underlying SLD to dump into GeoServer:
The end result: road_casing.sld
Fun with Bounding Boxes
I’ve been load-testing GeoServer quite a bit lately and do so I have to generate test plans that execute a number of WMS requests. In this case I wanted the tests to simulate a typical scenario of concentrated access centering on major cities.
The first step was to grab a dataset with city limits:
from geoscript.layer import Shapefile # builtup area from shp = Shapefile('world_boundaries/builtup_area.shp') # city names to include cities = ['CALGARY', 'VANCOUVER', 'DENVER', 'CHICAGO'] # world bounds in web mercator box = Bounds(-20037508,-19929239,20037508,19929239, 'epsg:900913')
Next we’ll loop through all the cities I’m interested in and get the city bounds in Web Mercator.
for city_name in cities: # union all the polygons into one geometry city_area = reduce(lambda x,y: x.union(y), [f.geom for f in shp.features("nam = '%s'" % city_name)]) # reproject to web mercator city_bnds = Bounds(env=city_area.bounds(), prj='epsg:3395').reproject('epsg:900913')
Last we generate bounding box tiles that intersect the city limits for the zoom levels I want available:
# go from min to max zooms for z in range(8, 15): # get resolution of this zoom level res = 1.0/pow(2,z) # tile our main bbox at this resolution within the city bounds for t in box.tiles(r, city_bnds): print '%f,%f,%f,%f\n' % (t.west, t.south, t.east, t.north)
The next big thing on the GeoScript roadmap is raster data support. It’s something I’m particularly excited about and happy to say that some of us at OpenGeo have been busy working on. Here is a quick sample of what’s to come:
Grab some raster data:
from geoscript.layer import GeoTIFF dem = GeoTIFF('sfdem.tif')
Do some analysis on it:
dem.bands() # [GRAY_INDEX] dem.extrema() # ([-9.999999933815811e+36], [1840.0]) histo = dem.histogram(low=0, high=1840)
Render the histogram:
from geoscript.plot import bar data = zip(map(lambda (x,y): y, h.bins()), h.counts()) bar.xy(data).show()
Build a color map and render the dem:
from geoscript.style import * # create the intervals intervals = range(0,2000,10) # interpolate some colors colors = Color('white').interpolate(Color('black'), n=len(intervals)) # build the color map and render draw(dem, ColorMap(zip(intervals, colors)))
And there you have it, raster analysis made easy with GeoScript. Remember that these are just a few examples of what’s possible with server-side scripting. Download the OpenGeo Suite and try it out for yourselves. We love feedback so if you have any questions or comments about GeoScript, or wrote something that you want to share please let us know. Leave a comment on this post or send us an email. Thanks for reading and happy scripting.