Web Maps 2: Python Services
Thus far we have been using OpenLayers JavaScript to add mapping functionality to the page, and in separate exercises you learned how to use Geo-Webservices (using MapServer) and consume these in a OpenLayers map client. You also learned how to create Python services, and used one to include a dynamic version of a bar chart. Now we will continue to also add a Python service to create a Thematic Map of the same data of country the provinces we used in the bar chart. We could go search for an OWS service that produces what we need, but instead we are going to create our own Python service. The mechanism necessary for this task is a so-called 'server-side script', written in CGI Python. If you have never used Python to create a CGI script before, make sure to first do the exercise we have created to introduce you to that...!
You might have also done the exercise on how to use scripts with a so-called REST interface. Using REST services falls outside the scope of this exercise, but you are of course welcome to add a REST interface to your services here also...!
Vector Layers using Python services
The data we need for the map service is the same we used for the bar chart, the world.prov_pop table of the exercises database. We need to generate data that can be used as input for the OpenLayers javascript library to map this on top of the OpenStreetMap layer. The most obvious data format to use for this is GeoJSON.
GeoJSON is a format in which simple geometry is expressed in JSON syntax. GeoJSON supports the geometry types
Point, LineString, Polygon, MultiPoint, MultiLineString, and MultiPolygon. Geometric objects with additional
attributes/properties are Feature objects. Sets of features are contained by FeatureCollection objects. It is
an IEEE standard: http://geojson.org/. In Python any
JSON, so also GeoJSON, is essentially a Dictionary
. GeoJSON is supported by almost all GIS software, including OpenLayers.
The listing below shows an SQL query that can be used to generate the data of the provinces for the service as a GeoJSON output.
1 2 3 4 5 6 7 8 9 10 11 12 |
SELECT json_build_object( 'type', 'FeatureCollection', 'features', json_agg( json_build_object( 'type', 'Feature', 'properties', jsonb_set(row_to_json(prov)::jsonb, '{geom}', '0', false), 'geometry', ST_AsGeoJSON(geom)::json ) ) ) FROM (SELECT gid, admin_name, pop_admin, geom, cntry_name FROM world.prov_pop) as prov WHERE cntry_name = 'Netherlands'; |
We do not explain all the parts of the SQL here, but note that we use the PostGres built-in functions
json_build_object()
and json_agg()
to aggregate the rows selected into a valid
GeoSJON Feature Collection.
The geometry of the provinces is converted from the internal PostGIS binary to GeoJSON, using the ST_AsGeoJSON
function in line 7. The jsonbset
function in line 6 is there to avoid
the geometry ending up as a long text description in the properties of each Feature: That would be superfluous as the geometry will already
be encoded as a GeoJSON geometry.
Run the query using pgAdmin to see it in action. The result should be similar in content to figure 1 below (it will be unformatted, to get that and the syntax-highligthing you'd have to load it into an 'intelligent' editor such as Notepad++ or PyCharm).

We are now ready to work on the service script. Proceed to create a new file using the code in Listing 2, and save it as /student/<<SNUMBER>>/services/provinces.py.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
import cgi, cgitb cgitb.enable() import psycopg2 import json ############# Get Data ################ # connect to the db db = psycopg2.connect("host='gisedu.itc.utwente.nl' port='5434' \ dbname='exercises' user='exercises' password='exercises'") # Get a cursor dbcursor = db.cursor() #create and execute query SQL = """ SELECT json_build_object( 'type', 'FeatureCollection', 'features', json_agg( json_build_object( 'type', 'Feature', 'properties', jsonb_set(row_to_json(prov)::jsonb,'{geom}','0',false), 'geometry', ST_AsGeoJSON(geom)::json ) ) ) FROM (SELECT gid, admin_name, pop_admin,geom, cntry_name FROM world.prov_pop) as prov WHERE cntry_name = 'Netherlands'; """ dbcursor.execute(SQL) # using json.dumps to create valid json, using fetchone()[0] to only get the first (and only) object jsonstring = json.dumps(dbcursor.fetchone()[0]) ############# write to Response Object ################ # first we need to output http header for JSON data: print ('Content-type: application/json; charset=utf-8') print ('') # this empty print is required! print ( jsonstring ) |
The provinces.py
script first loads necessary modules: cgi
to be able to act as a
CGI application, json
for handling JSON objects, psycopg2
for connecting to Postgres
DB's. Next, it connects to the DB and runs the query. Finally the results containing the geometries are sent
back to the client in JSON format [33–35]. It is important to state the appropriate content-type
for the response, as we do in line 33, so that the content of the response can be parsed properly
by the client application.
Test the new script in your browser to verify that it works.
Those of you that payed attention, would probably think: "could this not also be achieved by using a WFS...?" And they'd be correct! By using the MapServer's WFS capability to create GeoJSON output we'd have achieved the same result. We chose to use a Python service, because this is easier to parameterize later, as you will find out in the Challenge section...
Now that we have a service in place, we can use it to overlay the boundaries of Dutch provinces on the map. Our new service does not use a standard OGC service request/response structure (such as in a WFS), but it does however generate a data stream which is fully compliant, a correct GeoJSON object. We can therefore use OpenLayers to consume the service. Update the simple-website/scripts/viewer.js by adding the highlighted code below at the end:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
dotdotdot country_provinces = new ol.layer.Vector({ source: new ol.source.Vector({ url: '../services/provinces.py', format: new ol.format.GeoJSON({ defaultDataProjection: 'EPSG:4326', projection: 'EPSG:3857' }) }), name: 'Country Provinces' }); navMap.addLayer(country_provinces); dotdotdot |
To render the provinces layer, we have defined a new layer using OpenLayers ol.layer.Vector object (in line 3-12). One of the vector object's properties is url, which we use to load the data from the service we just created. The response data is structured as GeoJSON, so we use this as the format definition for the vector data [6]. As part of the format object we define the projection in which the original data is, EPSG:4326, and then the projection to be used for display, EPSG:3857 [7-8]. OpenLayers will take care of executing the corresponding transformation. Finally we specify a display name for this layer [11], and then add it to the map [13].
Reload the page to see how the polygons of the provinces are drawn on the map. Now the provinces are shown in a default style as we have not specified any style ourselves. To fix that, we just need to create, and use, a style object. Insert the following code in the simple-website/scripts/viewer.js file
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
dotdotdot var prov_style = new ol.style.Style({ stroke: new ol.style.Stroke({ color: 'MediumPurple', width: 2 }), fill: new ol.style.Fill({ color: 'rgba(147, 112, 219, 0.2)' }) }); country_provinces = new ol.layer.Vector({ source: new ol.source.Vector({ url: '../services/provinces.py', format: new ol.format.GeoJSON({ defaultDataProjection: 'EPSG:4326', projection: 'EPSG:3857' }) }), style: prov_style, name: 'Country Provinces' }); navMap.addLayer(country_provinces); dotdotdot |
Now we should have a more attractive provinces layer. If you refresh the page the result should look like the image in the figure below.

However, the map may not be as nicely zoomed in to the Netherlands as we show here. The user can do that of
course using the OpenLayers zoom tools, but it would be nice if the system does that... This can be achieved
by using the mapView.fit()
function of OpenLayers (see line 8 in the listing below). But
because the loading and rendering of the data is an asynchronous process, we need to wait a bit until
the loading is completed. That is achieved by running the fitting function as the result of a so-called
callback promise
: the promise .once
is fulfilled once the 'change'
event
has happened (line 6):
1 2 3 4 5 6 7 8 9 10 11 12 |
dotdotdot navMap.addLayer(country_provinces); // wait until the 'change' event is fired (when all data is loaded) country_provinces.getSource().once('change', function () { // then zoom in to extent: navMap.getView().fit(country_provinces.getSource().getExtent()); }); dotdotdot |
Now we have a web site that consumes various services. But they do so only for The Netherlands (or whatever other country you have hard-coded in your Python services SQL queries)... The next and last section contains an optional challenge: to create a parameterized version where the users can make their own choice of country to show the map and bar chart for.