PyGrunn: Geoprocessing with Python - Greg Kowal¶

Tags: pygrunn, python

(One of the summaries of the one-day 2014 PyGrunn conference in Groningen in the Netherlands).

Grek Kowal will talk about geoprocessing with GDAL/OGR. Geoprocessing is a lot about map projections. The earth is round, but not quite. And “round” needs to be mapped on a flat square screen, so you need to to reproject. Lots of ways to do that.

What do people do with geoprocessing? Analyzing sensor data, calculating optimal traffic routes, flood risks, etc. Often you have to take multiple data sources and merge them somehow to end up with a proper map.

What’s there in python?

• Shapely is nice and pythonic, but it cannot handle many data sources.

• ArcPy is proprietary.

• GDAL/OGR is a real swiss army knife: useful like hell, but ugly.

• QGIS is a graphical swiss knife, very handy for visually toying with your data.

GDAL is for rasters (bitmap), OGR for vectors (point, line, polygon). They’re now distributed together. (Note: if you manage to install it on OSX you’re a king: it is hard…).

Why are they so interesting? OGR supports 78 formats from postgresql to the Czech cadastral exchange data format. GDAL in turn supports 133 formates from PNG and JPEG to a “Swedisch Grid RIK”. That’s what he means with swiss army knife. They can read everything.

As a demo, he showed SoilScan. There’s a company that drives around with quads with a sensor box in the back. They get hired by farmers to scan their land for fertility and soil type and pollution. This allows farmers to use their fertilizer more economically.

The result is a bunch of data points with x/y locations and values for particles like for instance potassium-40, uranium-238 and thorium-232. These values are a pretty good indicator for the amount of fertilizard you need. He did some preprocessing on that data with a small “R” script, driven from Python. The result (basically a list of (x, y, value) lists) got fed into GDAL for creating a raster.

But a raster bitmap takes up much space and memory. So a vector file would be nicer, also for the later user interface. He fed the GDAL result into OGR and simply told it to convert it to a polygon. The polygon takes up much less memory than the raster representation.

The data that came out overflowed the boundary of the farmer’s field a bit. So he read in a shapefile with the farmer’s field and clipped the polygon results to that shape.

In the end he showed a leaflet map with the results. Nice.