Thursday, December 19, 2013

GRASS GIS Web Map Service with Pyramid

Lately I explored the capabilities of the Python bindings in GRASS GIS. Since I'm quite familiar with Python web frameworks it was an easy task to wrap the d.mon module and create a Web Map Service based on a GRASS location.

Here I'll show how I implemented the WMS using the Pyramid framework in a single file.

Below is the core function of the WMS, it is basically a wrapper for the d.mon module.
First of all the GRASS Python modules are imported and a new GRASS environment is started. After setting the current region from the bounding box parameters in the WMS request, a new Cairo monitor is started with d.mon. The output is written to a temporary file which will be returned.
def _grass_wms(layers=[], bbox=[], width=256, height=256):
    # Path to GRASS installation
    gisbase = "/usr/local/grass-7.0.svn"
    # Set the environment variable
    os.environ["GISBASE"] = gisbase

    # Append the GRASS Python directory to the path
    sys.path.append("%s/etc/python" % gisbase)

    # Import the GRASS scripts
    from grass.script import core as grass
    from grass.script import setup as gsetup

    # Open the mapset
    gsetup.init(gisbase,
                "/home/adrian/Data/GISBASE7",
                "Landsat_128047",
                "adrian")

    vector_layers = grass.list_strings("vect")
    raster_layers = grass.list_strings("rast")

    grass.run_command("g.region", w=bbox[0], s=bbox[1], e=bbox[2], n=bbox[3])

    # Create a temp file
    tempfile = grass.tempfile()
    filename = "%s.png" % tempfile

    grass.run_command("d.mon",
                      start="cairo",
                      width=width,
                      height=height,
                      output=filename)
    for layer in layers:
        if layer in raster_layers:
            grass.run_command("d.rast", map=layer, quiet=1)
        elif layer in vector_layers:
            grass.run_command("d.vect", map=layer, quiet=1, fcolor="0:0:255", color=None)

    grass.run_command("d.mon", stop="cairo")

    return filename

Above method is called by a Pyramid view that gets the layers, bounding box and image width and height parameters from the request and returns the image:
@view_config(route_name="wms")
def wms_view(request):

    layers = request.params.get("LAYERS", "").split(",")

    bbox = request.params.get("BBOX", "").split(",")

    width = request.params.get("WIDTH")
    height = request.params.get("HEIGHT")

    try:
        filename = _grass_wms(layers=layers,
                              bbox=bbox,
                              width=width,
                              height=height)
        f = open(filename, "r+")
        return Response(body=f.read(), content_type="image/png")
    except:
        pass

Furthermore I added a super simple index page which uses OpenLayers to show the WMS layer:
@view_config(route_name="index")
def index_view(request):
    body = """
<html xml:lang="en"
xmlns:tal="http://xml.zope.org/namespaces/tal"
xmlns="http://www.w3.org/1999/xhtml">
<head>
  <title>GRASS GIS Web Map Service</title>
  <script src="http://openlayers.org/api/OpenLayers.js" type="text/javascript"></script>
  <script type="text/javascript">
    var map, wms_layer;
    function init(){
        var map = new OpenLayers.Map("map-div",{
            projection: "EPSG:32648",
            maxExtent: new OpenLayers.Bounds(235080, 2132280, 257100, 2147400),
            numZoomLevels: 6,
        });
        var wms_layer = new OpenLayers.Layer.WMS("GRASS WMS", '/wms', {
            layers: "namngum5_patch,namngum5_bw,namngum5_smooth"
        },{
            singleTile: true,
        });
        map.addLayer(wms_layer);
        map.zoomToMaxExtent();
    }
  </script>
</head>
<body onload="init()">
    <div id="map-div" style="height: 600px; width: 800px;">
    </div>
</body>
</html>"""
    return Response(body=body, content_type="text/html", status=200)

Last but not least it needs some Pyramid / WSGI code to get the server running:
def main(global_config, ** settings):
    """ This function returns a Pyramid WSGI application.
    """
    config = Configurator(settings=settings)
    config.add_static_view('static', 'static', cache_max_age=3600)
    config.add_route("index", '/')
    config.add_route("wms", '/wms')
    config.scan()
    return config.make_wsgi_app()

Since this is a very quick implementation, there are some unsolved issues:
  • Reprojection on the fly is not supported, i.e. a WMS request must be in the same projection as the GRASS location
  • Pyramid requests are not thread-safe i.e. if there are several requests setting a different region in GRASS the monitor output gets mixed up
  • probably others as well
Last but not least two screenshots to illustrate the Web Map Service:
The index page shows the WMS layer
The same location and layers in GRASS GIS

No comments:

Post a Comment