Mapping the population density of New Zealand with QGIS, SQLite and TileMill

There are over 4.4 million people living in New Zealand, but they’re not evenly distributed across the country. On my travels, I’ve been to some very remote areas like Doubtful Sound in Fjordland with a permanent population of one! Where do the New Zealanders live? Let’s create a population density map of New Zealand.

Doubtful Sound

According to Lonely Planet, 63% of New Zealanders live on the North Island, 20% on the South Island, 10% in Australia, 5% in the rest of the world, and 2% are travelling! The area of New Zealand is 268,021 km2 (Norway has 385,252 km2), which means there are about 16.5 New Zealanders per km2. Norwegians have a bit more space, we're only 15.5 persons per km2.

Auckland

I want to map the population density as a choropleth map, using darker colors for higher densities. But which units or geographical areas should I use? The regions (territorial authorities) I mapped in my last blog post are too big to show the population density of the country. Instead, I decided to use area units, which are non-administrative areas that are in between meshblocks (the smallest geographic unit for which statistical data is collected by Statistics New Zealand) and territorial authorities. In urban areas, area units normally contain a population of 3,000–5,000.

The census is the official count of how many people and dwellings there are in New Zealand. It takes a snapshot of the people in New Zealand and the places where we live. The last census was in 2006. There was supposed to be a new census in 2011, but it could not be completed due the national state of emergency after the Christchurch earthquake. A new census is scheduled for 2013.

I therefore decided to download the area units used for 2006 census. These can be found on Koordinates.com. I used QGIS to removed ocean, lakes, inlets and tidal bay areas where no or very few people live – which gives me the shape of New Zealand the people recognize.

The yellow area shows the extent of the original dataset, while the red color shows the area units after removing osean, lakes, inlets and tidal bay areas.

I’ve also used the field calculator in QGIS to calculate the size (in square kilometers) of each area unit polygon.

The dataset is in the New Zealand Transverse Mercator 2000 (NZTM2000) projection, which is using a metric coordinate system. The formula to calculate the area in km2 is $area / 1,000,000. There are one million square meters in a square kilometer.

The largest area unit is Fjordland with 8287 km2. The smallest area unit is
Motuopae Island in Tauranga, with less than 0.03 km2.

We can save the dataset as a SQLite database by right-clicking on the layer name in QGIS and select “Save as…”.

Using SQLite will allow us to add the population data to the same file.

Now that we have the area units, we need some population data covering the same geographical areas. Luckily, getting statistical data in New Zealand is as easy as getting geospatial data. I used the beta release of NZ.Stat to download estimated population data for 2006-2011.

It was not possible to download the dataset as a CSV file when I tried, so I downloaded the Excel version. I then created a CSV file with this format (first line is the column names):

no;name;pop2006;pop2007;pop2008;pop2009;pop2010;pop2011
500100;Awanui;370;360;360;360;360;360
500202;Karikari Peninsula-Maungataniwha;4350;4320;4300;4320;4370;4350
500203;Taipa Bay-Mangonui;1610;1590;1570;1560;1550;1540
500204;Herekino;2020;2000;1980;1960;1960;1950
500205;Ahipara;1160;1160;1180;1170;1170;1180
500206;North Cape;520;500;490;480;460;450

I used SQLite Database Browser to import the population dataset to the file containing the area units (File -> Import -> Table from CSV file).

Finally, we can use TileMill to create population density map. You can create a layer from a SQLite database like this:

I'm using this SQL query to join the area units and the population data:


SELECT OGC_FID, au_no, au_name, geometry, area, pop2011, (pop2011/area) AS popdens FROM nz_area_units_2006_census AS area JOIN population ON au_no = no

The query is also calculating the population per km2 (pop2011/area), which will be our mapping variable. Also note that you need to specify the SRS projection string for the dataset (NZTM2000):

+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

I've used a 9-class yellow-orange-brown color scheme from ColorBrewer to represent increasing population density.

This is my CartoCSS:


The map looks like this:

Recommended reading: Working with TileMill and SQLite (MapBox)

I'm going to create an interactive version with a map legend in my next blog post. Stay tuned!

Mapping the regions of New Zealand with MapShaper and Leaflet

There are 16 regions of New Zealand, and during my trip I’ll be visiting 15 of them (I’m saving the Taranaki region for my next visit). This blog post shows how you can download and simplify region boundaries, and add them to a map as an interactive layer.

The region is the top tier of local government in New Zealand. Eleven are governed by an elected regional council, while five are governed by territorial authorities (the second tier of local government) which also perform the functions of a regional council. The current regions came about in 1989, when they replaced more than 700 boards which had been formed in the preceding century.

The geographic extents of the regions are largely based on drainage basins, following drainage divides such as the Southern Alps. 

The dataset I’m using is from Statistics New Zealand, the national statistical office. Koordinates.com has made life easier by collecting geographic data from Statistics New Zealand on their portal. I’ve downloaded a polygon dataset from the last census in New Zealand (2006).

The dataset includes more than 60,000 coordinates used to draw the boundaries for each region. As I’m going to display the vector data directly in the web browser, I should simplify the dataset by reducing the number of coordinates while preserving the shape of the regions. Leaflet has built-in support for vector simplification (to make the map more responsive), but you shouldn’t rely on this alone as the data still has to be loaded and processed by Leaflet.
  
You can do vector simplification with QGIS (Vector -> Geometry Tools -> Simplify geometries), but there is a better tool, MapShaper (about), which gives you more options and preserves topology (shared borders). I loaded the shapefile into MapShaper and set the simplification level to 50% using the Special Visvalingam method. This will reduce the number of coordinates by 90%, but you still have to zoom in very far to see the difference.

Next, I’m converting the simplified shapefile to GeoJSON with ogr2ogr:

ogr2ogr -t_srs EPSG:4326 -f "GeoJSON" -lco COORDINATE_PRECISION=3 nz-regional-councils-simplified.json nz-regional-councils-simplified.shp

By restricting the number of coordinates (COORDINATE_PRECISION=3) to the map scale or number of zoom levels, you can further reduce the size of the dataset.

I’m assigning the GeoJSON object to a JavaScript variable (regions), and this is the code required to create a Leaflet layer with custom style, region highlighting and popups:

The map looks like this:

Fullscreen map

Mapping New Zealand: Creating a road map

I’ve travelled around New Zealand by bus, car, boat, and by foot. In a previous blog post, I created a “Where I’ve Been Map” using markers for each of the places I’ve stayed. Let’s add some lines showing the route between the markers.

I've destroyed my National Geographic map while travelling around the New Zealand. How can I add the roads travelled to my digital map?

Travelling around New Zealand can be quite an experience. This is a well known sight.

I’ll concentrate on the roads in this post, and deal with boat and foot tracks later. As I didn’t carry a GPS when travelling by bus or car, I had to find the road lines from a different source. I’m using road data from LINZ for my map, but I’m sure you can achieve the same with OpenStreetMap data.

First, I created an empty shapefile in QGIS for my road data. Then I marked and copied the roads I’ve travelled from the LINZ shapefile, and pasted the road segments into my own shapefile. I also had to cut road lines to remove parts where I’ve not been. Lastly, I merged the road segments together. QGIS is a great tool for tasks like this.

Editing road lines in QGIS.

I used ogr2ogr to convert the shapefile into GeoJSON, which is the vector format supported by Leaflet:

ogr2ogr -t_srs EPSG:4326 -f "GeoJSON" -lco COORDINATE_PRECISION=3 nz-tour-road.json nz-tour-road.shp

I then assigned the the GeoJSON object to a variable:

var roads = {
  "type": "FeatureCollection",
  "features": [...]
};

To create a GeoJSON layer with custom styles in Leaflet only requires a few lines of code:

var roads = L.geoJson(roads, {
  style: {
    color: '#333',
    weight: 1.5,
    opacity: 1
  }
});

As I wanted to combine the road lines with the place markers I decided to use the Leaflet LayerGroup:

This is my new “Where I’ve Been” map:

Fullscreen map

Kingston road along Lake Wakatipu.

Towards Glenorchy and the Southern Alps.

Mapping New Zealand: Clustering DOC Huts with Leaflet

In New Zealand, long distance walking or hiking for at least one overnight stay is known as tramping. There is a great network of over 950 huts throughout New Zealand operated by the Department of Conservation (DOC). I’ve just stayed in four of the huts while tramping the Abel Tasman Coast Track, and you need a lifetime to reach them all. It’s much quicker to map them.

DOC huts in the Abel Tasman National Park.

Awaroa Hut

I was fortunate to visit the geospatial unit at the DOC office in Wellington. They have a lot of ineresting conservation projects going on, and it seemed to be a great place to work in New Zealand, especially when you can combine digital work with practical field work.

You can download DOC's geospatial data for free on their data portal. Unfortunately, the datasets are not available as shapefiles, so it might require some wizardry to extract the data from KML. Luckily, the Koordinates guys have done this already, and you can download the shapefile from their website. The version I’m using was updated 10 September 2012.

The dataset contains 967 huts, and adding them all to my Leaflet map will make it too cluttered to be useful. Instead, I wanted to try the great animated marker clustering plugin, created by Dave Leaver here in New Zealand.

I start off by converting the shapefile to GeoJSON with ogr2ogr:

ogr2ogr -t_srs EPSG:4326 -f "GeoJSON" -lco COORDINATE_PRECISION=4 nz-doc-huts.json doc-huts.shp

I could use the GeoJSON directly, but I converted it to a more compact format to save some bandwith:

var data = [
  [-43.3767,170.5685,"Scone Hut",1],
  [-42.9494,171.7047,"Sudden Valley Biv",0],
  [-41.5237,172.5604,"Granity Pass Hut",1],
  [-44.9261,168.2144,"Steele Creek Hut",0],
  [-44.6221,168.4481,"Earnslaw Hut",0],
  [-45.3852,167.6192,"Luxmore Hut",4],
  ...
];

The last number corresponds to the hut categories:

0. Basic Hut/Bivvy
1. Standard Hut
2. Serviced Hut
3. Serviced-Alpine Hut
4. Great Walk Hut

This is the code you need to create clustered huts with custom markers:

The JavaScript map function, which converts one array into another, is not supported in Internet Explorer 8 and lower. If you need to support prehistoric browsers, use a for-loop instead. The icons I’ve used (which could have been more descriptive) are from Map Icons Collection.

The interactive map looks like this:

Fullscreen map

I especially like the animated feature of this clustering plugin (try zooming in and out), as it makes it more obvious for the map user what a cluster really is.

Whariwharangi Hut, built as a farmhouse in about 1896.

Mapping New Zealand: Where are the hot and cold springs?

New Zealand has a large number of cold and hot springs, where water flows to the surface of the earth from underground. Where can you find these springs? This blog post will show you can visualise a spring point dataset on a Leaflet map. 

You can find spring point dataset at LINZ Data Service, which contains 147 significant springs either by size or location. I downloaded the dataset as a shapefile, and converted it to GeoJSON using ogr2ogr:

ogr2ogr -f "GeoJSON" -lco COORDINATE_PRECISION=4 nz-spring-points.json nz-spring-points-topo-150.shp 

The resulting GeoJSON FeatureCollection can be assigned to a JavaScript variable:

var springs = {"type":"FeatureCollection","features":[{...

... and visualised with Leaflet's GeoJSON layer and L.circleMarker:

Here is the result:


Fullscreen map

Te Waikoropupū Springs in Golden Bay, discharging 14,000 litres of water per second, - the largest freshwater springs in New Zealand.

The horizontal visibility of the water in the springs has been measured at an average of 63 metres, second only to sub-glacial water in the Antarctic.

How to control your Leaflet map with URL parameters

I needed a way to control the appearance for my Leaflet map of New Zealand, especially when embedding a map with an iframe. Leaflet have no native support for URL parameters for, like the Permalink control for OpenLayers (there is a plugin + Leaflet-hash). This blog post will show you how you can add support for URL parameters for map position (latitude/longitude), zoom level and active layers.

First, I'm defining my layers in an JavaScript object:

var layers = {
  'Seafloor': L.tileLayer('tiles/nz-seafloor/{z}/{x}/{y}.png'),
  'Topographic': L.tileLayer('tiles/nz-topo/{z}/{x}/{y}.png'),
  'Tour': L.cartoDB('http://thematicmapping.cartodb.com/api/v2/sql?q=SELECT * FROM new_zealand_tour')
};

You can use the same object to populate the Leaflet layers control. Next, I'm extracting the URL parameters with a regular expression (from Papermashup):

var params = {};
window.location.href.replace(/[?&]+([^=&]+)=([^&]*)/gi, function(m, key, value) {
  params[key] = value;
});

All URL parameters are now stored as key/value pairs in the params object. We need to create an array of active map layers:

if (params.layers) {
  var activeLayers = params.layers.split(',').map(function(item) {    
    return layers[item];
  });
}

The array map function will create a new array with a reference to the visible layers on the map. Please note that the map function is not supported in Internet Explorer < 9 (you can easily rewrite the function with a for loop to support this browser).

We can now create a map with our URL parameters:

var map = new L.Map('map', {
  center: [params.lat || -41.2728, params.lng || 173.2996],
  minZoom: 4,
  maxZoom: 12,  
  zoom: params.zoom || 6,
  layers: activeLayers || [layers.Seafloor, layers.Topographic]
});

I'm using the Logical OR (||) operator to add default values. Lastly, I'm adding my map layers to the layers control, allowing users to toggle the visibility:

L.control.layers({}, layers).addTo(map);

This is the full example:

We can now control the map med URL parameters like (the code in this example is changed to support base layers and overlays):

http://earthatlas.info/nz/?lat=-42.4&lng=168&zoom=8



http://earthatlas.info/nz/?lat=-39.296&lng=174.065&zoom=11&layers=Springs


Easy peasy!

Creating a "Where I've Been" map with Leaflet and CartoDB

I wanted to create a "Where I've Been" map of my trip around New Zealand. The map should be easy to update through a simple web interface. Previously, I've used Google Fusion Tables for tasks like this. Now, I want to try CartoDB, which is Google Fusion Tables on steroids.

I used the CartoDB free plan to create a simple table containg the name and location of the places I've been:

You can easily query your data with the HTTP-based SQL API:

http://thematicmapping.cartodb.com/api/v2/sql?q=SELECT * FROM nz_tour

This request will return all the data in the default CartoDB JSON-format. Add format=geojson, if you want the data in the GeoJSON format supported by Leaflet:

http://thematicmapping.cartodb.com/api/v2/sql?q=SELECT * FROM nz_tour&format=geojson

The SQL API also supports JSONP, allowing you to query your data directly from JavaScript in the web browser. Just add a callback parameter:

http://thematicmapping.cartodb.com/api/v2/sql?q=SELECT * FROM nz_tour&format=geojson&callback=myFunction

Leaflet don't have native support for JSONP or CartoDB, but there are a few plugins around (1, 2). It's also very easy to extend Leaflet with the functionality you want, so I decided to create a new CartoDB layer class by extending the L.GeoJSON class:

This class rewrites the initialize function and adds a JSONP function (you can also use a library like jQuery). This is how you can create a new layer with this class:

L.cartoDB('http://thematicmapping.cartodb.com/api/v2/sql?q=SELECT * FROM nz_tour')

This will create a layer with the default Leaflet markers:

If you want custom markers and popups, all the functionality available for the GeoJSON layer can also be used for the CartoDB layer:

The interactive map looks like this:

Fullscreen map

You can add a marker to the map by adding a new row to your CartoDB table. Dead simple!