Stay Connected with the Boundless Blog

Hacking with the GeoScript Library

The 3.0 release of the OpenGeo Suite added a variety of software features and improvements, including spatial processing capabilities. With scripting support, languages like JavaScript and Python could be leveraged to fit your needs. These scripting languages expose the GeoScript library, allowing users to leverage all GeoScript functions in processing scripts. While OpenGeo Suite 3.0 ships with many useful processes, the additional server-side scripting allows you to deploy many more to fit your enterprise’s specific wants and needs.

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.

Data Conversion

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='', 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]


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 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:

writeSLD(style, 'road_casing.sld')

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

# 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,%fn' % (t.west, t.south, t.east, t.north)

What’s Next?

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:


# ([-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())

Build a color map and render the dem:

from 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.