Developing a Kernel Density Service with GeoTrellis

Developing a Kernel Density Service with GeoTrellis

To go along with the recent Azavea Atlas post, we’re going to build a WMS service for the Western Jackalope using GeoTrellis.

GeoTrellis is a high performance geoprocessing engine and programming toolkit developed here at Azavea.

There’s an example showing the Western Jackalope sightings running at http://207.245.89.238/labs1-demo/index.html.

[iframe src=http://207.245.89.238/labs1-demo/index.html]

Scala Environment

We’re using the same made-up dataset that John used in the Atlas Article. The first step is getting our environment setup. If performance is key, we’ve found the Oracle JVM to be a bit faster than OpenJDK. For simplicity, we’re going to use the OpenJDK found in the standard apt repo:

[github file = “/ahinz/labs1-demo/blob/blog/INSTALL.md” start_line = “1” end_line = “2”]

To build scala projects we use “sbt”. If you don’t have it, you can install it:

[github file = “/ahinz/labs1-demo/blob/blog/INSTALL.md” start_line = “2” end_line = “3”]

To get moving we need two boilerplate files. The first is our project definition (which lives in the root in build.sbt)

[github file = “/ahinz/labs1-demo/blob/blog/build.sbt”]

The second file is needed to explain to GeoTrellis what classes we want exposed and on what port by updating the settings in the application.conf file

[github file = “/ahinz/labs1-demo/blob/blog/src/main/resources/application.conf”]

To make sure everything is looking good we’re going to setup a simple hello world service. We’re using Jersey/JAXRS for web services in this tutorial but you should feel free to use whatever.

Testing Service

[github file = “/ahinz/labs1-demo/blob/blog/src/main/scala/demo/azavea/Hello.scala”]

At this point we can do “./sbt run” and see our service startup by going to http://localhost:8888/hello and you should see “Hello World”. Jersey uses annotation to determine what methods to expose. We use the “@Path(…)” annotation to represent that path of a server (“/hello” in this example). Decorating our method with @GET exposes it via GET requests.

Formatting Our Data

Before we can analyze the data we need to grab some data and prep our environment. We use GeoTools to parse our shapefile in to GeoTrellis features.

[github file = “/ahinz/labs1-demo/blob/blog/src/main/scala/demo/azavea/Resource.scala” start_line = “1” end_line = “52”]

The last few lines set up our GeoTrellis server executor. Usually GeoTrellis reads a catalog file to determine where relevant raster data is on the server. Since we’re only generating rasters (not reading them) we used a blank catalog. The server is responsible for running operations.

The last few lines setup our color handling. All of the regular color ramps have their alpha values set to 100%. For our heat map we want low values to fade away. In the above case the lowest values are dropped completely (0x00000000) and the remaining values are set to 50% opacity.

Kernel Density

Here’s a kernel density service that we’ll be using:

[github file = “/ahinz/labs1-demo/blob/blog/src/main/scala/demo/azavea/Resource.scala” start_line = “53” end_line = “112”]

We start out grabbing our bounding box (xmin, ymin, xmax, ymax), the width and height, size of our kernel and the kernel’s spread (standard deviation in meters). “Context” is a global object in “demo.azavea” that we use to store some static stuff, like our feature and server.

GeoTrellis provides a few operations for converting strings to relevant objects:

More can be found in the “geotrellis.rest.op.string” package

A common kernel to use for Kernel Density operations is the Gaussian, which we create using CreateGaussianRaster. In our example we allow the standard deviation to be set by the service and set the amplitude to 100.

The main show is the KernelDensity operation. We’re using the feature set that we loaded from the shapefile. The second argument is a function to convert our feature value in a number. In our case each feature is a Point[Int] so we use the identity function.

This kernel density server represents everything in meters and kilometers so we need to scale the values appropriately before we dive in to the rendering (spreadOp and sizeOp).

The last 10 lines render the PNG and and return it. Since we’re generating tiles on the fly we’re not going to use quantile breaks because the tile edges won’t line up. Instead, we create our own break array and use RenderPng directly.

Resources

For more information checkout http://geotrellis.github.com/ or chat with us at #geotrellis on freenode.

The code can be found on github at https://github.com/geotrellis/labs1-demo