Fishnets, Honeycombs and Footballs – Better spatial models with hexagonal grids.

I love the unique mathematical properties of hexagonal tiling for spatial analysis and the inherent possibility to create beautiful bi-variate hex maps.

GEOLYTIX asked me to share some of my thoughts about hex grids. What makes them special and why are they better suited for all aspects of our work than fishnets. There will be an in-depth look at the design of hex maps with Leaflet and Openlayers in a series of future blog posts and I am planning to do a workshop at FOSS4G 2017 in Boston. In this post I want to give a brief overview of regular grid sampling and the benefits of hexagonal grids over rectangular grids.

 

Hexagons are the closest shape to a circle that can form a regular grid.

Regular tiling of an Euclidean surface (eg. a map) can be achieved with triangles, squares (fishnets), or hexagons. The latter two can be further tessellated into triangles. But why is hexagonal tessellation best for maps?

The honeycomb conjecture itself was only proven in 1999 by American mathematician Thomas C. Hales. There is however plenty of accessible material on the nature of hex structures. My favourite is this TED.ed video.

It all boils down to the high perimeter / area ratio. Structural strength. That which is good for physical models should be good for an analytical model as well. But how does the need for less bee wax relate to making better maps? A map is a Cartographic (flat) representation of a distinct part of the spherical Earth. We must therefore flatten its surface in order for features to be drawn on a flat screen.

We ain’t going back to make you use cathode ray tubes for our web maps but fewer edges mean that less material needs to be bend when we put a hexagonal grid over a curving surface.

Think about the white bits on a golden age football.

Every point on the surface of a hexagon is on average closer to the centre compared with points in a triangle or rectangle. This causes less distortion of the shape. Here is a more relevant representation which shows a hexagonal grid on an azimuthal projection which uses the location of GEOLYTIX offices in London as centre.

For the scientifically minded there is a peer reviewed article from Birch et al on ‘Rectangular and hexagonal grids used for observation, experiment and simulation in ecology’.

Turning heterogeneous areas into regular grid surfaces

The next two grid maps show European countries split into equally sized hexagonal and rectangular grids.

Apart from Luxembourg being completely ignored in the rectangular grid due a shape mismatch at the selected resolution it can be seen that the hexagonal grid matches the international boundaries better than the rectangular presentation.

The next screenshot shows the distribution of registered unemployed in communities near Vienna. Very little variation can be seen in the rural communities outside the capital.

We get a much better insight into the actual distribution of unemployed labourers if we apply the figures from the communities to a hexagonal grid (1km cell height) which has been created from the Global Human Settlement (GHS) population grid. We see that large stretches of land are forests and do not represent the populated areas in which we are interested for a social model.

Unlike fixed reporting geometries, grid cells can be sized according to the required model scale. If we zoom out of the region shown in the images above we can no longer distinguish small urban communities. A hexagonal grid with 4km cell height however clearly shows the population distribution across the extended region.

Gridded networks

An important advantage of the hexagonal grid is the unambiguous definition of nearest neighbourhood: Each hexagon has six adjacent hexagons in symmetrically equivalent positions (Birch et al. 2000). In contrast, the rectangular
grid has two different kinds of nearest neighbour: orthogonal neighbours sharing an edge and diagonal neighbours sharing only a corner.

Networks can be better represented with a hexagonal grid. The next two images show the London Tube network as a rectangular grid, sometimes called a Manhatten grid, as well as a hexagonal grid. The colours represent the different tube lines.

Being able to connect to diagonal neighbours improves the representation of the network in a hexagonal grid.

There are numerous blog posts which have pointed out the benefits of different types of grids. This post would not have been possible without the work and dedication of others. I want to point to this excellent gis.stackexchange post as well as this blog post by Matt Strimas-Mackey. Matt does some great work with R and highlights the square grid’s benefits of orthogonal simplicity and computational ease. With this in mind:

We choose to create hex tiles in this decade and do the other things, not because they are easy, but because they are hard.

OSM tileserver on a vanilla Ubuntu 16.04 droplet.

I want to write down a few simple steps which helped me to setup an OpenStreetMap tileserver in under 2 hours. There are quite a few good tutorials out there which helped me to get things rollin. Without this great tutorial by Xiao Guoan I wouldn’t have been able to have this running yet. I tried to get the switch2osm tutorial for Ubuntu 14.04 working on Ubuntu 15.04 a few months back and failed terribly. The experience of compiling Mapnik from source is harrowing. Don’t waste your time and use the 16.04 packages as described in this article.

The tutorial was completed on a vanilla Ubuntu 16.04 droplet from DigitalOcean. The spec for the server is 2GB RAM with 40GB SSD. At $20 per month this is very good value from DigitalOcean.

We love OpenStreetMap data at GEOLYTIX and hope to roll out our own styles in the future.

Let’s get started…

Open a shell as root user and create a new user first. I call this user ‘osm’ and add the user to the sudo group. I add my public key to the authorized_keys file in order to open the shell for this user in the future.

My preferred shell editor is nano but you can use vi, vim or whichever editor you are most comfortable with.

adduser osm
usermod -aG sudo osm
su osm
sudo mkdir ~/.ssh
sudo nano ~/.ssh/authorized_keys

Paste your key into the authorized_keys file. This will allow you to open a shell as the OSM user in the future.

My first sudos are an apt update followed by an upgrade and the installation of a few dependencies which we will need later on.

sudo apt update sudo apt upgrade
sudo apt install curl unzip autoconf gdal-bin libtool

I also install fail2ban and setup a UFW firewall blocking all but http, https and ssh.

Installing Postgis

Next install PostgreSQL, PostGIS and a couple of tools which we need to load OSM data to our server.

sudo apt install postgresql postgresql-contrib postgis postgresql-9.5-postgis-2.2

Let’s switch to the user ‘postgres’; Open psql and change the password for this user. \q quits psql.

sudo -u postgres -i
postgres$ psql
\password
\q

Create a Postgresql user ‘osm’ and a database for the OSM data load. I call the database ‘england’. Add the postgis and hstore extensions to the newly created database.

postgres$ createuser osm
postgres$ createdb -E UTF8 -O osm england
postgres$ psql -c "CREATE EXTENSION hstore;" -d england
postgres$ psql -c "CREATE EXTENSION postgis;" -d england
postgres$ exit

Getting the data

I get OSM data dumps from Geofabrik. I suggest to create a directory to store the raw data.

mkdir ~/osm_data
wget http://download.geofabrik.de/europe/great-britain/england-latest.osm.pbf -O ~/osm_data/england.pbf

Getting the data into PostGIS

Install osm2pgsql and allocate a swap file.

sudo apt install osm2pgsql
sudo fallocate -l 2G /swapfile
sudo chmod 600 /swapfile
sudo mkswap /swapfile
sudo swapon /swapfile
sudo nano /etc/ssh/ssh_config

Check these settings in the ssh_config file.

Host *

    ClientAliveInterval 30 TCP

    KeepAlive yes

    ClientAliveCountMax 99999

Then restart the sshd service.

sudo service sshd restart

The load into PostGIS itself will require a style. I create a new folder for OSM styles and get Andy‘s default OSM carto style.

mkdir ~/osm_styles
cd ~/osm_styles
git clone https://github.com/gravitystorm/openstreetmap-carto ~/osm_styles/openstreetmap-carto

Now load the data into PostGIS using osm2pgsql.

osm2pgsql --slim -d england -C 3600 --hstore -S ~/osm_styles/openstreetmap-carto/openstreetmap-carto.style ~/osm_data/england.pbf

Install Apache2, Mapnik, Carto

Apache2 is the tile server’s web server.

Mapnik is the rendering engine at the core of making OSM tiles.

CartoCSS is a syntax language for OSM styles.

sudo apt install node-carto libmapnik3.0 libmapnik-dev mapnik-utils python-mapnik apache2 apache2-dev

Build and install mod_tile

mod_tile is a plugin for Apache2 which ties the web server and mapnik render engine together. We need to make and install mod_tile from source.

git clone https://github.com/openstreetmap/mod_tile.git
cd mod_tile/
./autogen.sh
./configure
make
sudo make install
sudo make install-mod_tile

Build the style

We also need to build the OSM style. Style require some world shapefiles. OpenStreetMap-Carto styles comes with a script to get and index these shapefiles.

cd ~/osm_styles/openstreetmap-carto
./get-shapefiles.sh
nano project.mml

We need to make some edits in the project.mml file before we can build the style.

Replace ‘gis’ with database name which you chose earlier e.g. ‘england‘.

In nano this is easily done with the find and replace function (crtl + \). There should be around 70 instances of the database name in the project.mml.

We can build the style now using the carto tool.

carto project.mml > style.xml

Configure renderd

Renderd is the rendering demon for OSM tiles.

sudo nano /usr/local/etc/renderd.conf

In the renderd configuration check the host definition, the location of the OSM style, and the mapnik plugin.

[default]

XML=/home/osm/osm_styles/openstreetmap-carto/style.xml

HOST=localhost

[mapnik]

plugins_dir=/usr/lib/mapnik/3.0/input/

Copy the sample rendered.init file into the SystemV init directory in order to start the service when the system boots. Make the file executable with chmod. And open nano to do some config edits.

sudo cp ~/mod_tile/debian/renderd.init /etc/init.d/renderd
sudo chmod a+x /etc/init.d/renderd
sudo nano /etc/init.d/renderd

Check following configs in the renderd startup.

DAEMON=/usr/local/bin/$NAME

DAEMON_ARGS="-c /usr/local/etc/renderd.conf"

RUNASUSER=osm

Make sure that mod_tile is good to go for renderd before we start the service.

sudo mkdir -p /var/lib/mod_tile
sudo chown osm:osm /var/lib/mod_tile
sudo systemctl daemon-reload
sudo service renderd start
sudo service renderd enabled

Let’s edit the mod_tile.load file for Apache2.

sudo nano /etc/apache2/mods-available/mod_tile.load

In this file we tell Apache2 to load the load the mod_tile module.

LoadModule tile_module /usr/lib/apache2/modules/mod_tile.so

Enable the module and open the default site config for Apache2.

sudo ln -s /etc/apache2/mods-available/mod_tile.load /etc/apache2/mods-enabled/
sudo nano /etc/apache2/sites-enabled/000-default.conf

In this file we set the actual web service and cross origin access.

<VirtualHost *:80>
  ErrorLog ${APACHE_LOG_DIR}/error.log
  LoadTileConfigFile /usr/local/etc/renderd.conf
  ModTileRenderdSocketName /var/run/renderd/renderd.sock
  ModTileRequestTimeout 0
  ModTileMissingRequestTimeout 30
  Header set Access-Control-Allow-Origin "*"
</VirtualHost>

Now we need to restart Apache2 and renderd.

sudo service apache2 restart
sudo service renderd restart

Done! Give yourself a pad on the back and open tile zero in a browser. Here is my OSM world image tile. Note… I tend to use https for pretty much everything and so should you. I am also bad with IP addresses and give all my servers a name and URL redirect.

IMHO the easiest way to check whether you tile server plays well with a web map framework is to setup a jsFiddle or codePen, etc.

It doesn’t matter whether you use Openlayers or Leaflet. I like Leaflet because it is so simple and small. Plus, there is finally a proper Leaflet release.

Here is a quick jsFiddle for the tile server which I just setup.

All good? Well there is one last thing which I wasn’t able to workout. For me renderd does not start when I reboot the server. I have tried all of the suggestion which I could find on the web but so far no cigar.

Here is my gis.stackexchange post in this regard. It would be awesome if you can help.

WFS-T with OpenLayers 3.16

This post was originally published for OL 3.5 in June 2015. There is no good reason not to use the latest version of OpenLayers and I am updating the code examples today for OL 3.16.

I am currently using GeoServer 2.8. The data store being a PostGIS 2.1 database.

For my test service I use a very simple table with an ID and geometry only. The geometry is defined as geometry, no type nor projection. It is important that the geometry field is called geometry. Otherwise inserts may create records with empty geometry fields. A constraint must be set on the ID or GeoServer will not be able to insert records into the table.

CREATE TABLE wfs_geom
(
  id bigint NOT NULL,
  geometry geometry,
  CONSTRAINT wfs_geom_pkey PRIMARY KEY (id)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE wfs_geom
  OWNER TO geoserver;

CREATE INDEX sidx_wfs_geom
  ON wfs_geom
  USING gist
  (geometry);

At the core of the OL javascript snippets is the ol.format.WFS.writeTransaction function which takes 4 input parameter. The first 3 parameter define whether the data should be inserted, updated, or deleted from the data source. The 4th parameter takes the form of ol.format.GML and passes on information about the feature type, namespace, and projection of the data.

The writeTransaction node must be serialised with an XMLSerializer to be used in a WFS-T post.

The three use cases (insert/update/delete) and the AJAX call are shown in the following code snippet.

var formatWFS = new ol.format.WFS();

var formatGML = new ol.format.GML({
    featureNS: 'https://gsx.geolytix.net/geoserver/geolytix_wfs',
    featureType: 'wfs_geom',
    srsName: 'EPSG:3857'
});

var xs = new XMLSerializer();

var transactWFS = function (mode, f) {
    var node;
    switch (mode) {
        case 'insert':
            node = formatWFS.writeTransaction([f], null, null, formatGML);
            break;
        case 'update':
            node = formatWFS.writeTransaction(null, [f], null, formatGML);
            break;
        case 'delete':
            node = formatWFS.writeTransaction(null, null, [f], formatGML);
            break;
    }
    var payload = xs.serializeToString(node);
    $.ajax('https://gsx.geolytix.net/geoserver/geolytix_wfs/ows', {
        service: 'WFS',
        type: 'POST',
        dataType: 'xml',
        processData: false,
        contentType: 'text/xml',
        data: payload
    }).done(function() {
        sourceWFS.clear();
    });
};

Inserts are called from the drawend event of the OL interaction. The .clear() function on the WFS source reloads the source after each transaction. This ensure that feature IDs are correct for new features and that deleted features are removed from the view.

The ID of a modified feature is stored and the update transaction is posted once the feature is unselected. In order to successfully post an update transaction the boundedBy property must be stripped from the feature properties. A clone of the feature is used to achieve this.

map.addInteraction(interactionSelect);
interaction = new ol.interaction.Modify({
    features: interactionSelect.getFeatures()
});
map.addInteraction(interaction);
map.addInteraction(interactionSnap);
dirty = {};
interactionSelect.getFeatures().on('add', function (e) {
    e.element.on('change', function (e) {
        dirty[e.target.getId()] = true;
    });
});
interactionSelect.getFeatures().on('remove', function (e) {
    var f = e.element;
    if (dirty[f.getId()]) {
        delete dirty[f.getId()];
        var featureProperties = f.getProperties();
        delete featureProperties.boundedBy;
        var clone = new ol.Feature(featureProperties);
        clone.setId(f.getId());
        transactWFS('update', clone);
    }
});

A complete working example is hosted in this jsFiddle.

A previous working example for OL 3.5 is still hosted on Google App Engine.

The complete code example for the OL 3.5 version is on GitHub.

I plan to update my Google App Engine projects and GitHub repositories after I complete my switch to CodeEnvy.

OpenLayers 3.16 – Three flavours of Geoserver WFS as ol.layer.Vector

Today I want to recap on the different ways of adding a vector layer from a Geoserver WFS. I used OL 3.5 when I published this post in May 2015. Since then I have switched to OL 3.16 and decided to update this post rather than to write a new post. There is really no reason not to be working with the latest version of OpenLayers. If you are familiar with WFS vector sources in OpenLayers 3.x or if you used my previous code snippets then please review the update notes. There have been a few significant changes between OL 3.5 and OL 3.16.

The examples are based on the vector-wfs example from the OpenLayers page.

We create a vector layer with a WFS source like so:

var layerWFS = new ol.layer.Vector({
  source: new ol.source.Vector({
    loader: function(extent) {
      $.ajax('http://demo.opengeo.org/geoserver/wfs', {
        type: 'GET',
        data: {
          service: 'WFS',
          version: '1.1.0',
          request: 'GetFeature',
          typename: 'water_areas',
          srsname: 'EPSG:3857',
          bbox: extent.join(',') + ',EPSG:3857'
        }
      }).done(function(response) {
        layerWFS
        .getSource()
        .addFeatures(new ol.format.WFS()
          .readFeatures(response));
      });
    },
    strategy: ol.loadingstrategy.bbox,
    projection: 'EPSG:3857'
  })
});

We use a bbox loading strategy and use the extent in the loader function for the WFS source.

By default the WFS request will return XML from Geoserver. Inside the .done() function we define WFS format object which will be used to read the features from the WFS response and add these to the source.

Here is a working JSFiddle for this code.

As an alternative to the WFS format we can directly request GeoJSON from GeoServer. GeoJSON will be more compact and easier to parse.

var layerWFS = new ol.layer.Vector({
    source: new ol.source.Vector({
        loader: function (extent) {
            $.ajax('http://demo.opengeo.org/geoserver/wfs', {
                type: 'GET',
                data: {
                    service: 'WFS',
                    version: '1.1.0',
                    request: 'GetFeature',
                    typename: 'water_areas',
                    srsname: 'EPSG:3857',
                    outputFormat: 'application/json',
                    bbox: extent.join(',') + ',EPSG:3857'
                }
            }).done(function (response) {
                layerWFS
                .getSource()
                .addFeatures(new ol.format.GeoJSON()
                  .readFeatures(response));
            });
        },
        strategy: ol.loadingstrategy.bbox,
        projection: 'EPSG:3857'
    })
});

We define the outputFormat as ‘application/json’ inside the AJAX call.

In order to read the features from the GeoJSON response we use a GeoJSON format object inside the .done() function.

Here is a working JSFiddle.

An alternative to the JSON format is the JSONP format which allows us to receive data from Geoserver without changing the CORS access on the webserver. JSONP is by default disabled and must be enabled in Geoserver.

The AJAX call must configured as shown below for the callback with padded JSON.

var layerWFS = new ol.layer.Vector({
    source: new ol.source.Vector({
        loader: function (extent) {
            $.ajax('http://demo.opengeo.org/geoserver/wfs', {
                type: 'GET',
                data: {
                    service: 'WFS',
                    version: '1.1.0',
                    request: 'GetFeature',
                    typename: 'water_areas',
                    srsname: 'EPSG:3857',
                    outputFormat: 'text/javascript',
                    bbox: extent.join(',') + ',EPSG:3857'
                },
                dataType: 'jsonp',
                jsonpCallback:'callback:loadFeatures',
                jsonp: 'format_options'
            })
        },
        strategy: ol.loadingstrategy.bbox,
        projection: 'EPSG:3857'
    })
});

window.loadFeatures = function (response) {
    layerWFS
    .getSource()
    .addFeatures(new ol.format.GeoJSON().readFeatures(response));
};

We set the output format in the AJAX call to be ‘text/javascript’. We don’t use a .done() function on the AJAX call but define dataType, jsonp, and jsonCallback parameter. The response from the WFS call will be sent to the jsonCallback function. Inside the load features function we use again a GeoJSON format object to read features from the response and add these to the source of the WFS layer.

Here is a working JSFiddle.

If you have problems loading data I recommend to use the Chrome developer tools (F12). A common error is the use of a http source on a https page.

wfs_https_error

Filter for XHR on the network tab of the debug tools to inspect the details of your WFS calls.

wfs_chrome_debug

Please drop me a comment if this helped you, if you have questions or if you can recommend further improvements to these examples.

Next time… WFS-T.

 

 

AJAX mail service on Google App Engine

Setting up a mail service in the Google App Engine is simple enough following the tutorial on the Google Cloud Platform. I create a seperate class for the mail service like this:

import webapp2
from google.appengine.api import mail

class MailService(webapp2.RequestHandler):
  def post(self):
    mailFrom = self.request.get('mailFrom')
    mailMsg = self.request.get('mailMsg')
    if not mail.is_email_valid(mailFrom):
      self.response.write('wrong mail address')
      return
    message=mail.EmailMessage(
      sender='dennisargeomatica@gmail.com',
      subject='Mail from website')
    message.to='dbauszus@gmail.com'
    message.body=(
      'Mail from: %s \n'
      'Message: %s \n'
      %(mailFrom,mailMsg))
    message.send()
    message=mail.EmailMessage(
      sender='dennisargeomatica@gmail.com',
      subject='Your mail to Argeomatica')
    message.to=mailFrom
    message.body=(
      'Thank you, we have received your mail. \n'
      'Message: %s \n'
      %(mailMsg))
    message.send()
    self.response.write('mail sent')

app = webapp2.WSGIApplication(
  [('/mailService', MailService)],
  debug=True)

Coming from an ASP.net background I am used to work with update panels which allow me to do use a submit input control without the need of rendering the complete page after the submit.

The jinja2 framework, although incredibly fast, does not allow for a similar functionality. Instead we will use an AJAX call from the client to the GAE mail service.

Creating an AJAX call is dead simple with jQuery.ajax(). I call my GAE mail service from javaScript like so:

function mailService() {
  mailFrom = $('#tbFrom').val();
  mailMsg = $('#tbMail').val();
  $.ajax({
    type : 'POST',
    url : 'mailService',
    data : {
      mailFrom:mailFrom,
      mailMsg:mailMsg
      },
    success: function(response) {
      $('#mailStatus').html(response);
      }
    })
}

The function will grab the email address and message body from the DOM input elements and POST the data to the mail service. The success event handler will then receive confirmation of the process from the service and display the message on the page.

The complete site is hosted here.

The complete code is on GitHub.

Geoserver + PostGIS on OpenShift

I wanted to have my own Geoserver development environment for free (as in speech) when I found this great blog post by Steve Citron-Pousty:

Build your own Google Maps (and more) with Geoserver on OpenShift.

I am fairly familiar with Geoserver and wasn’t interested in the sample data provided on the GitHub page for this tutorial. In my setup I use the Geoserver 2.7.0 web archive which I downloaded from the Geoserver project page.

I signed up for a free account with OpenShift which gives us 3 free gears to use. These gears will hibernate after a day or so but it is possible to upgrade the account to bronze (no monthly charge) which means that the 3 free gears won’t hibernate. Upgrade to bronze requires you to submit your credit card details and is not yet available in all regions.

I run the setup from a Windows 8 machine with Powershell. Pre-requisites are Ruby, GIT, and RHC.

To create the application on Openshift use following command:

$ rhc app create geoserver -s tomcat7 postgresql-9

geoserver being the name of the new application and tomcat7 and postgresql-9 the cartridges which are used on the Openshift platform.

A local git clone will be created. From this clone remove the pom.xml since we are not interested in a Maven build.

From the git root open a secure shell into the application:

$ rhc ssh

In the application shell create a directory for the geoserver application data.

mkdir $OPENSHIFT_DATA_DIR/geoserver_data

From the application shell open a pSQL shell.

psql

In the pSQL shell enable the PostGIS extension for the default database.

CREATE EXTENSION postgis;

Exit the pSQL shell with CTRL+D and exit the application shell with ‘exit’.

Rename the Geoserver web archive to ROOT.war and move the file into the webapps directory of the GIT clone.

In .openshift\action_hooks create the file pre_start_jbossews-2.0 with the content line: ‘export CATALINA_OPTS=-DGEOSERVER_DATA_DIR=${OPENSHIFT_DATA_DIR}geoserver_data’

Geoserver would be useless if you cannot access the services from other client applications. In order to enable CORS for Tomcat edit the .openshift\config\web.xml and add following filter:

<filter>
  <filter-name>CorsFilter</filter-name>
  <filter-class>org.apache.catalina.filters.CorsFilter</filter-class>
</filter>
<filter-mapping>
  <filter-name>CorsFilter</filter-name>
  <url-pattern>/*</url-pattern>
</filter-mapping>

Now add, commit and push these changes to OpenShift.

$ git add -A
$ git commit -am "prep work geoserver"
$ git push

The setup may take a few minutes since the Geoserver web archive is rather large (~60mb).

Once the application has fully started you can login to the admin console at http://geoserver-{your account}.rhcloud.com/web/

Change the admin password first, thereafter you can start to setup your workspaces, stores, styles, and layers.

I use QGIS to upload my data into the PostGIS database which we created. In order to access the database from my local machine I need to enable port forwarding like so.

$ rhc port-forward

Hello world!

Trying to do a million things at the same time there is a risk you will forget that one command which makes your process run as intended. Neurodegeneration is a serious threat to productivity. The primary intention of this blog is a personal future reference to selected bits and bobs of work which I am currently doing. If this may be of use to others, even better. Anyways, please let me know in a comment if a post is helpful to you or whether you need help with an related issue.