Web Map Service with PyCairo and Pyramid

From time to time I have the chance to work with PyCairo, the Python bindings to the cairo graphic library. I appreciate this library because it allows to draw graphics primitives so super fast and thus it is very suitable to use it in web projects.
After my trials with GRASS GIS as WMS backend I thought I could create a raster web map service with improved performance by using PyCairo instead of GRASS GIS as renderer.

Zoom in - it is a super fast WMS!

Since I wanted to limit the dependencies as much as possible I had the following constraints:

• Input raster dataset must be a PNG file, since cairo natively does not read nor write another raster file format
• The input raster datasets needs to be georeferenced with an accompanying world file to avoid an additional GIS library like GDAL. A world file contains exactly the same transformation parameters that are needed to place the raster image to the cairo context thus it can directly applied as transformation matrix.

As with my GRASS GIS trials I wrapped everything in a simple Pyramid application.

To avoid to read the input raster on every request it is read at server start up together with the world file and stored in a global variable as cairo.SurfacePattern.

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()
  
    # Check pycairo capabilities
    if not (cairo.HAS_IMAGE_SURFACE and cairo.HAS_PNG_FUNCTIONS):
        raise SystemExit('cairo was not compiled with ImageSurface and PNG support')
    # Read the input files
    raster_data = cairo.ImageSurface.create_from_png(os.path.join(os.path.dirname(__file__), "data", "LC8128047_w_MASK.png"))
    world_file = open(os.path.join(os.path.dirname(__file__), "data", "LC8128047_w_MASK.wld"), 'r')
    world_file_content = world_file.read()
    
    # Create a cairo pattern from the input raster layer
    # and declare the pattern global
    global raster_pattern
    raster_pattern = cairo.SurfacePattern(raster_data)

    # Parse the world file
    lines = world_file_content.split("\n")
    georef_scale_x = float(lines[0])
    georef_scale_y = float(lines[3])
    georef_trans_x = float(lines[4])
    georef_trans_y = float(lines[5])
    # Setup a new transformation matrix for the georeferenced raster layer
    scaler = cairo.Matrix()
    scaler.scale(georef_scale_x, georef_scale_y)
    scaler.translate(georef_trans_x / georef_scale_x, georef_trans_y / georef_scale_y)
    scaler.invert()
    raster_pattern.set_matrix(scaler)
        
    # Set resampling filter
    raster_pattern.set_filter(cairo.FILTER_FAST)
    
    # Start the WSGI application
    return config.make_wsgi_app()

Having this SurfacePattern a web map service view is very easy. The request parameters need to be read and a new cairo context created and transformed. Then I just paint the SurfacePattern to the context and return the newly created in-memory image. It is not even necessary to write it to the hard disk which saves time as well.

@view_config(route_name="wms")
def wms_view(request):
    
    start_time = time()

    # Read the WMS parameters
    layers = request.params.get("LAYERS", "").split(",")
    bbox = [float(b) for b in request.params.get("BBOX", "").split(",")]
    width = int(request.params.get("WIDTH"))
    height = int(request.params.get("HEIGHT"))
    
    # Create a new cairo surface with requested size
    canvas = cairo.ImageSurface(cairo.FORMAT_ARGB32, width, height)
    ctx = cairo.Context(canvas)
    # Calculate the scale in x and y direction
    scale_x = width / (bbox[2] - bbox[0])
    scale_y = height / (bbox[1] - bbox[3])
    # Create a forward transformation matrix
    ctx_matrix = cairo.Matrix(y0=(-1) * bbox[3] * scale_y,
                          x0=(-1) * bbox[0] * scale_x,
                          yy=scale_y,
                          xx=scale_x)
    ctx.set_matrix(ctx_matrix)

    # Set the raster pattern as source and paint the canvas
    ctx.set_source(raster_pattern)
    ctx.paint()
    
    # Create a new response object
    response = Response(content_type="image/png", status="200 OK")
    try:
        # Write the canvas as png to the response body
        canvas.write_to_png(response.body_file)
        return response
    except:
        pass
    finally:
        log.debug("'cairo' took: %ss" % (time() - start_time))

Of course it is finally necessary to wrap up everything in a simple OpenLayers web map:

@view_config(route_name='index', renderer='templates/mytemplate.pt')
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>PyCairo 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(208800.0, 1967189.0, 430230.0, 2188289.0),
                numZoomLevels: 8,
            });
            var wms_layer = new OpenLayers.Layer.WMS("GRASS WMS", '/wms', {
                layers: "LC8128047_w_MASK"
            });
            map.addLayer(wms_layer);
            map.zoomToMaxExtent();
        }
    </script>
</head>
<body onload="init()">
    <div id="map-div" style="bottom: 0px; left: 0px; position: absolute; right: 0px; top: 0px;">
</div>
</body>
</html>"""

    return Response(body=body, content_type="text/html", status=200)

Conclusions: Once more PyCairo didn't disappoint! To return a single raster tile it only takes around 0.01 to 0.04 seconds which is really fast.
As with the previous GRASS GIS service reprojection on the fly is not supported.

Finally another screenshot of the resulting map with a cloudless Landsat 8 scene:

The cloudless Landsat 8 scene