thematic mapping blog

Norway will open its topographic datasets to the public!

2013-04-30T21:14:00.001+01:00

Today, 30th April 2013, is a milestone of the mapping history of Norway. Together with two governmental ministers, the Norwegian Mapping Authority announced they will open its topographic datasets to the public, free of charge. Norway will follow progressive countries like Denmark, Finland and New Zealand, considering geographic data as a public good.

The topographic datasets at 1:50,000 scale will be freely available later this year, together with address, road and cadastre data. A database of 950,000 place names was released to the public one month ago. Hopefully, digital elevation data will be free as well, but it's not yet stated.

Cheers!

Cross-country skiing in Norway

2013-04-28T17:47:00.001+01:00

It's been a great winter in Norway. Pleny of snow, lots of sunshine and good skiing conditions, except from my last trip where we caught in a storm and luckily found our way back to the hut. Outdoors was what got me into mapping in the first place. A hike is much more enjoyable with a good trekking map, and on my last trip I really experienced how crucial a GPS can be.

As I spent my free time skiing there's been no time for blogging, so I decied to post a few photos instead. Cross-country skiing in Norway is highly recommended!

Mt. Gaustatoppen, where you can see one sixth of Norways mainland from the top.

A view of Rondane, the oldest national park in Norway.

Mountain pasture used for the grazing of sheep and cattle in the western part of Norway.

By the Hardangerfjord, the third largest fjord in the world.

Skiing in Dovre National Park.

PS! I'll start a new blog series on mapping with D3.js in June. Stay tuned!

How to minify GeoJSON files?

2012-11-18T21:48:00.000Z

You can't do web mapping these days without knowing your GeoJSON. It's the vector format of choice among popular mapping libraries like Leaflet, D3.js and Polymaps. Size matters on the web, especially if you want to distribute complex geometries, like the world's countries. The challenge is even bigger if you want to target mobile users - or support web browsers with poor vector handling (IE < 9). This blog post will show you how to minify your GeoJSON files before sending them over the wire.

The first thing you should do is to generalize your vectors so they don't contain more detail than you need. In a previous blog post, I was able to remove 90% of the coordinates without loosing to much detail for map scale I wanted to use. This will of course have a great effect on the file size.

Today, I'm going to use country borders from the Natural Earth dataset. These datasets are already generalized for different scales (1:10m, 1:50m, and 1:110 million), so I'll use them as they are. The 1:110m (small scale) and 1:50m (medium scale) shapefiles will cover the needs for the thematic world maps I plan to make:

The 110m and 50m country polygons shown in QGIS.

Let's open the datasets in QGIS. If you look at the attribute table you'll see that each dataset contains 63 attributes, which makes them very versatile. For your web maps, you probably need just a few of the attributes, and you should remove the ones you don't need. I'm keeping the country name and the ISO 3166-1 country codes (alpha-2, alpha-3, and numeric), which can be used to link country geometries to statistical data.

Only keep the attributes you need.

Next, we can convert the shapefiles to GeoJSON with ogr2ogr:

ogr2ogr -f "GeoJSON" -lco COORDINATE_PRECISION=1 ne_110m_admin_0_countries.json ne_110m_admin_0_countries.shp

ogr2ogr -f "GeoJSON" -lco COORDINATE_PRECISION=2 ne_50m_admin_0_countries.json ne_50m_admin_0_countries.shp

The important thing is that I'm only keeping one decimal (coordinate precision) for the 110m dataset, and two decimals for the 50m dataset, which is sufficient for my map scales. This will reduce the size of the GeoJSON files by more than half. The size of the 110m GeoJSON is now 207 kB and the 50m version is 1,897 kB. But we can do better.

The files contains a lot of whitespace, which is waste of space. I planned to use Sublime Text to remove the whitespace, but it were not able to handle the 50m GeoJSON file, so I switched to Notepad++. I used these regular expressions:

Find: "([^a-z.]) "
Replace: "$1"
This will remove all whitespace which is not succeeding a letter or a dot, which are present in country names.

Find: "\n,"
Replace: ","
Remove line breaks (keeping some for readability).

Find: "\.0([,\]])"
Replace: "$1"
Remove trailing zeros.

This will reduce the file size of the 110m GeoJSON from 207 to 156 kB, without loosing any data quality. More than 400k of whitespace characters was removed from the 50m GeoJSON file, reducing the file size from 1,897 to 1,481 kB.

If your web server is supporting gzipping on-the-fly, the 110m GeoJSON will end up being 45 kB and the 50m version will be 430 kB. Not bad!

And if this is too much work, you can always download the final GeoJSON files on thematicmapping.org.

NB! Mike Bostock’s TopoJSON would allow us to compress the GeoJSON even more, while preserving topology (shared borders between countries) - but we would need to use a map client supporting the format. Looks promising!

Mapping New Zealand - a summary

2012-11-04T19:18:00.004Z

I've had a fantastic two months study trip to New Zealand. Unfortunately, I had to go back to Norway this week to fill up my bank account - just when the summer was arriving in New Zealand. I'm going to miss the beautiful country with its great people.

New Zealand is the perfect country to map, as an isolated country surrounded by a vast ocean, and because of all the free data available. I hope my Mapping New Zealand blog series has been useful for others as well:

It's a lot of exciting things happening on the New Zealand mapping scene, and the Kiwis are very welcoming people. I want to thank all the nice people I met on my journey, who was very willing to share their knowledge and experiences:

Organisers and participants at GeoCart'2012
The Koordinates.com guys
PlanetInAction (Paul van Dinther), you should check out his great Google Earth based helicopter simulator.
Organsiers and participants at the Mountain Cartography Workshop
The geospatial team at Landcare Research in Palmerston North.
The mapsmiths at Geographics - the best cartographers in New Zealand!
The geospatial team at Department of Conservation in Wellington.
Brent Wood and Kevin Mackay at NIWA.
NewTopo (Geoff Aitken).

I've been a real map nerd in New Zealand, but I had lots of time to explore the country too. Here are a few of my photos:

New Leaflet plugin to handle multiple TileMill layers

2012-11-04T18:12:00.003Z

My setup for the population density map of New Zealand made it easy to create new choropleth maps with New Zealand census data. This blog post explains how you can use Leaflet to switch between multiple interactive layers created with TileMill.

I wanted to create a map of the social geography of New Zealand, using the Index of Deprivation from the Department of Public Health, University of Otago. I downloaded a Excel sheet containing data for the 2006 census area units, which I also used for my population density map. I simply added the data to the same SQLite database, and created the map using the same techniques described in two previous blog posts (1, 2).

The Index of Deprivation is constructed from nine Census 2006 variables, and provides a summary deprivation score from 1 to 10. A score of 1 is allocated to the least deprived 10 percent of areas, and 10 is allocated to the most deprived 10 percent of areas. You'll find more information about the index in the Atlas of Socioeconomic Deprivation in New Zealand.

Wax allows you to easily add an interactive TileMill map to Leaflet or other mapping libraries. Adding more than one interactive layer is not that straightforward, so I wrote a Leaflet plugin:

This plugin allows you to switch between various layers (interactive or not) and it will automatically load and display map legends, and remove the elements when they're no longer needed. It's easy to use the plugin:

Include a wax property when you create a tile layer, containing a TileJSON object or an URL to a TileJSON file.

var population = L.tileLayer('tiles/nz-popden/{z}/{x}/{y}.png', {
attribution: 'Statistics New Zealand',
wax: 'tiles/nz-popden.tilejson'
});

After you've created the map you simply add:

L.wax(map);

That's it! :-)

You can try to switch between the various basemaps using the layers control below:

Fullscreen map

Creating map labels with TileMill

2012-11-04T16:40:00.002Z

It's one obvious thing missing from my maps of New Zealand: map labels. This blog post will show you how to create a transparant layer with map labels with TileMill, which can be added as an map overlay in Leaflet or other mapping libraries.

To create map labels, you need a point dataset containg at least a position and the label text. As we'll see later, information about type of place, importance etc. will also be useful. For my New Zealand maps, I'm using a dataset from LINZ Data Service: NZ Geographic Names (Topo, 1:500k)

I started by using QGIS to convert the shapefile into SQLite database (by right-clicking on the layer name in QGIS and select “Save as…”). This allows us run SQL queries against the data in TileMill. The dataset includes three attributes or columns, - name, size and a code describing the type of place:

I'm opening the SQLite file in TileMill:

I'm using this SQL query to load the data from the SQLite database:

SELECT * FROM nz_labels ORDER BY size DESC

This will sort the labels after size, in descending order. I'm assuming that larger size means more important labels, and this will instruct TileMill to start with the most important labels (there is not enough space to show all labels on all zoom levels). Remember also to include the SRS projection string for the data source (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

CartoCSS includes a provides several ways to style map labels or text, and I haven't got the time to explore all the possibilities. This is the my CartoCSS for zoom level 5-7:

I start by defining the fonts I'm going to use. It's a convention in cartography to label water features (lakes, sea, rivers) in serif italic faces, but small serif fonts on screen can be hard to read.. New high resolution screens, like Apple's Retina display, will probably make serif fonts more popular for web maps. You can use TypeBrewer to test various font combinations.

I'm starting with the places marked as "METR" (probably metropolitan areas, although using this term in New Zealand sounds a bit strange...). I'm aligning the labels (text-horizontal-alignment) to place them over the ocean, so they're not obscuring the map. A halo is added around the labels (text-halo-fill / text-halo-radius) to make the text visible on top of various basemaps.

Auckland, Wellington, Christchurch and Dunedin are the only labels shown on zoom level 5 (in the Web Mercator projection). I'm then adding towns on zoom level 6, and populated places on zoom level 7.

Labels shown on zoom level 6.

The minimum distance between each label (text-min-distance) is 30 pixels to avoid label collisions. This means that many towns or populated places will not be drawn because it's not enough space. More labels will be revealed as you zoom in. Although I sorted the features after the size attribute, it's probably not enough to make a good selection of labels for each zoom level. I'm sorry if I left out your town! :-)

I'm continuing like this until zoom level 12, gradually adding more labels for features like suburbs, lakes, mountains etc.

The map labels were exported as a separate map using the MBTiles format, and added to my Leaflet map as a map overlay. This is the result:

Fullscreen map

Exploring the MapBox stack: MBTiles, TileJSON, UTFGrids and Wax

2012-11-04T12:51:00.001Z

In my last blog post, we created a population density map of New Zealand using QGIS, SQLite and TileMill. Today, we’re going to publish this map to the web using various MapBox inventions. I'll also show you how to publish an interactive TileMill map on your own web server using some PHP and JavaScript wizardry.

I love MapBox. The team behind this platform has created a series of new specifications, allowing us to create fast, good looking and interactive maps. The downside is the limited support for other map projections than Web Mercator.

TileMill allows you to add legends and tooltips to your maps. I’ve added a legend to my population density map with a HTML snippet describing the map and the color scale.

The tooltip shows when the user hovers over or clicks on the map. It allows us to show dynamic content - additional data, images, charts - for each map feature. I want to show the name, total population, area and population density for each feature:

The data fields for the layer are wrapped in curly Mustache tags. These tags will be replaced by data when you interact with the map. You can use the full Mustache template language.

The easy way to publish this map is to upload it to MapBox Hosting, and use the embed code provided. If you want to publish your map on your own web server, this is an alternative route:

To export an interactive map from TileMill, you need to use the MBTiles format. This is an innovative SQLite-based format specification capable of storing millions of map tiles in a single file. The format is also supported by various 3rd-party applications, and I'm sure we'll see a greater adoption in the future.

Within the MBTiles file, the map legend, the tooltip template and information about map extent, zoom levels etc. is stored in a format named TileJSON. This is also an open specification, providing a consistent way of describing a map, making it easier to load and display a map the way it’s meant to be seen. The TileJSON for my map looks like this:

If you add interactivity to your map (tooltips), your MBTiles file will also include the most impressing part of the MapBox specifications: UTFgrids. This JSON-format allows us to add thousands of interactive points or polygons through interactivity data grids, and it will even work in older browsers with limited support for vector data.

So how do we turn our MBTiles file into an interactive map? Previously, I've used MBUtil to extract the contents from MBTiles into a directory structure. But by doing this, we loose the benefits of the MBTiles format, like storing a map in a single file and dealing with redundant images. What we need is a script on our web server that will extract content from our MBTiles file on demand. I decided to try a PHP script from infostreams (this is probably not the most scaleable solution). The script supports the full MBTiles specification, including TileJSON and UTFGrids. Installation is simple: just put the .php file and the .htaccess file in the same directory as your .mbtiles files. The .htaccess file includes a rule that rewrites requested URLs on the fly, so the map data is available un URLs like:

TileJSON: http://earthatlas.info/nz/tiles/nz-popden.tilejson
Map tile: http://earthatlas.info/nz/tiles/nz-popden/9/504/320.png
UTFGrid: http://earthatlas.info/nz/tiles/nz-popden/9/504/320.json

So when we have our backend sorted, how can we recreate our interactive map with Leaflet or other JavaScript mapping libraries? This is way the MapBox team created Wax, which is a client implementation of the MBTiles interaction specification. You just include the wax script together with your mapping library of choice, and then you can add interactivity with a few lines of code:

I've also done some extra JavaScript coding to allow switching between various interactive map layers. I'll save that for a later blog post.

The Leaflet map looks like this (there seems to be an issue with the latest Wax distribution and Google Chrome):

Fullscreen map

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

2012-10-26T11:58:00.002+01:00

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

2012-10-22T07:10:00.003+01:00

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

2012-10-20T10:53:00.001+01:00

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

2012-10-20T09:15:00.001+01:00

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 huts = [
[-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?

2012-10-08T13:20:00.002+01:00

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

2012-10-08T09:32:00.001+01:00

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

2012-10-08T02:09:00.003+01:00

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!

Mapping New Zealand: Exporting tiles from TileMill

2012-09-30T12:45:00.001+01:00

I wanted to publish a tiled version of my New Zealand maps. I first tried to upload my maps directly to MapBox.com from TileMill, but their free plan (50 MB upload storage) was not enough for all my zoom levels. This blog post describes an alternative, and a bit more complicated, way to publish your maps from TileMill.

TileMill uses a clever format called MBTiles. This is an open specification from MapBox capabale of storing millions of tiles in a single SQLite database. You can download your maps as MBTiles from the TileMill interface. This was easy to do for my topographic map, and a bit more tricky for my seafloor map which crosses the antimeridian. TileMill don't allow you to select an area for export which crosses the antimeridian, without making a full roundrip around the Earth:

Almost 45,000 tiles are needed to cover this big area for zoom level 0 to 9, while only 4700 tiles are needed to cover the area with data. Although TileMill has support for reduntant tiles, it's not the result I want. So I downloaded the map as two different MBTiles, one for each side of the antimeridian.

There are ways you can serve tiles directly to the web from your MBTiles, but I haven't had time to check how to do this with my shared hosting provider, Bluehost. Instead, I decided to extract all the images from the MBTiles and serve them directly from image folders. You can use MBUtil, also from MapBox, to achieve this:

mb-util nz-seafloor-1.mbtiles tiles/nz-seafloor-1
mb-util nz-seafloor-2.mbtiles tiles/nz-seafloor-2

The tiles are stored in z and x subdirectories such that they have a filesystem path like 0/0/0.png - a format that is supported by most web mapping clients. After the tile export, I just merged the two seafloor folders into one.

To create an interactive map is Leaflet is very easy. Just upload the tile folders to your web server, and write a few lines of JavaScript:

The map looks like this:

Fullscreen map

The map is centered at the "Geographic Center of New Zealand", which is here in Nelson.

The geographical "centre of New Zealand" allegedly lies on a hilltop near the city centre of Nelson.

We're heading off for the Abel Tasman Coast Track tomorrow, and there won't be any blog updates the next week :-)

Mapping New Zealand: Doing magic with TileMill

2012-09-28T16:26:00.000+01:00

Finally, I got some time to play with TileMill. I've added some textures to my shaded relief map, and I'm quite happy with the result. Here are a few samples:

Tongariro National Park

Marlborough Sounds

Westland National Park

Milford Sound

All the data I've added are from NZ Topo 500 from LINZ Data Service, and I've tried to find natural looking colors for my map stylesheet:

Creating a shaded relief map of New Zealand

2012-09-28T09:54:00.001+01:00

After creating a shaded relief map of New Zealand's seafloor, and it's time to move onshore. I'm going to use the same slope-enhanced hillshade technique. It was a bit tricky, but I think I've found a workflow that works.

While we got a 250 meter resolution Digital Elevation Model (DEM) of New Zealand's seafloor, you can download a 25 meter resolution DEM for the North and South Island from Landcare Research. The free datasets are available from the LRIS portal, where you can also read about the accuracy of the data. I got an invitation to visit Landcare Research when I was passing through Palmerston North. They do a lot of open source stuff, and they have some exciting projects going on!

Also, if you have a budget, you should definately consider buying the 8 meter resolution DEM from the best cartographers (or mapsmiths which they like to call themselves) in New Zealand - Geographix. They also provide some free samples for you to play with. You will be amazed by the level of detail in their well-crafted maps! They were kindly offering me an office space while I was visiting Wellington, and it was great to be part of such an inspiring environment - and what a view they have from their office on the edge of the Botanic Garden:

Wellington seen from the Geographx office and studio in the old Dominion Observatory.

Let's move on with our elevation model from Landcare Research. I first tried to merge the two DEMs for the North and South Island, but I ran into memory problems at a later stage because of the high resoltion. The commands below are covering the North Island, but the same set of techniques applies to the South Island DEM. I downloaded the GeoTIFF version of the datasets.

The datasets are in the New Zealand Transverse Mercator 2000 (NZTM2000) projection. As this projection is not supported by the map studio I want to use, TileMill, we need to reproject it into Web Mercator using gdalwarp:

gdalwarp -t_srs EPSG:3785 -tr 76.4370282745361 76.4370282745361 -r bilinear lris-nzdem-north-island-25-metre-GTiff/1-0001-0001.tif nz-dem-north-11.tif

At the same time I'm downscaling (-tr) the dataset so it has the same resolution as zoom level 11 in the widely used Web Mercator tiling scheme. Where do I get this 76.4370282745361 number come from?

So to cover the world at zoom level 11 you need almost 4.2 million tiles, and it becomes 4 times more for every zoom level you add to the party. Luckily we're only mapping a small part of the world, New Zealand. Web Mercator consider the world to be a square, 40,075,016.68 meters wide and tall (the length of the Equator). With a tile size of 256 pixels you get the a maximum resolution of 40,075,016.68 / 256 = 156,543.03390625. You get the map scale by dividing the resolution by 0.00028 (pixel size). The rest is simple math. The scale is only correct at the Equator. For 41° south latitude, running through the Cook Strait between the North and South Island, the map scale for zoom level 11 should be around 1:360,000 if my calculations are right (multiply by 1.325). Read more about the map projection and tiling scheme here.

Let's create a color relief from the DEM. I've created a typical color scheme which progresses from greens for lower elevations up through yellows/browns, and on to grays and white at the highest elevations:

gdaldem color-relief -co compress=lzw nz-dem-north-11.tif color-relief.txt nz-dem-north-11-relief.tif

Next, I want to get rid of the green ocean by only keeping the areas within the dataset I created in my last blog post (I've splitted the shapefile in two, one for each island):

gdalwarp -co compress=lzw -dstalpha -cutline shape/nz-north-webmerc.shp nz-dem-north-11-relief.tif nz-dem-north-11-relief-cut.tif

To create a nice looking hillshade, I'm going to apply the same technique as I did for the seafloor map (where you get all the details). First we create a "normal" hillshade:

gdaldem hillshade -alt 60 nz-dem-north-11.tif nz-dem-north-11-hillshade.tif

Then we create a slopeshade:

gdaldem slope -co compress=lzw nz-dem-north-11.tif nz-dem-north-11-slope.tif

gdaldem color-relief -co compress=lzw nz-dem-north-11-slope.tif slope-ramp.txt nz-dem-north-11-slopeshade.tif

Then we use Mapnik XML to combine the hill- and slopeshading:

nik2img.py nz-dem-north-hillslopeshade.xml nz-dem-north-11-hillslopeshade.tif -d 10497 14188

(for the South Island DEM the dimension -d is: 12637 14895)

We also to some cutting to get rid of the ocean:

gdalwarp -co compress=lzw -cutline shape/nz-north-webmerc.shp nz-dem-north-11-hillslopeshade.tif nz-dem-north-11-hillslopeshade-cut.tif

We can combine the color relief and the enhanced hillshade in TileMill (I'm keeping them separate so we can apply the hillshade to different textures):

We can also add New Zealand's seafloor to the map:

Mapping New Zealand: Editing and merging shapefiles with QGIS

2012-09-26T10:00:00.000+01:00

The data you want to map is not always in the format you need. I was looking for a dataset containing polygons for all land areas (islands) in New Zealand (1:500,000 scale), but I couldn’t find it on the LINZ Data Portal. So I tried to create it myself from the data I found:

Update Oct. 8th 2012: LINZ has added new polygon datasets for all land areas in New Zealand (see comments below), so you don't have to do the job described in this blog post. The techniques for editing and merging shapefiles are still valid.

NZ coastlines (Topo, 1:500k)
This dataset (nz-coastlines-topo-1500k.shp) contains coastlines for North Island, South Island and Stewart Island, - but not for offshore islands (most notably the Chatham Islands) and the smaller islands around the coast. The coastlines are represented as lines, while I need the data as island polygons. This is how I converted the lines into polygons with QGIS:

Vector -> Geometry Tools -> Lines to polygons

I named the new dataset "nz-main-island-polygons.shp".

(I wasn't able to do the same with the 1:250,000 version of the coastlines dataset, as the coastlines are divided into multiple line segments that didn’t merge in the correct order. Please let me know if anyone are able to convert this dataset into polygons! At 1:50,000 scale, LINZ is providing a dataset containing polygons for all islands.)

I’ve not found data for the offshore islands at 1:500,000 scale, but the smaller islands around the coast are available:

NZ island polygons (Topo, 1:500k)
This dataset (nz-island-polygons-topo-1.shp) also contains inland islands (islands in lakes), which I don’t want. I made a copy of this dataset named "nz-coast-island-polygons.shp" and removed the inland islands manually with "Select features by freehand" in QGIS.

I now want to merge the two shapefiles into one. This can easily be achieved in QGIS:

Vector -> Data Management Tools -> Merge shapefiles to one

We now got a polygon shapefile (nz-land-polygons.shp) containing all land areas of New Zealand, except the missing offshore islands.

Lastly, I also need this dataset in the Web Mercator projection (for my next blog post). I’m using ogr2ogr:

ogr2ogr -s_srs EPSG:2193 -t_srs EPSG:3785 nz-land-polygons-webmerc.shp nz-land-polygons.shp

(Written on a rainy day at the Hopewell lodge in Marlborough Sounds. A higly recommended place to releax a few days!)

Mapping New Zealand: Dealing with the antimeridian in TileMill

2012-09-26T09:36:00.001+01:00

I want to play with my new seafloor map in TileMill, which is a great map design studio. It’s built on top of Mapnik, but instead of XML you can style your maps with CartoCSS, which is much more user friendly. The biggest limitation of TileMill is that it only supports the Web Mercator projection (even though the underlying rendering engine, Mapnik, handles multiple map projections very well). Hopefully, TileMill will add support for other map projections and custom tiling schemes in the future.

Web Mercator is clearly not the ideal map projection to use for countries far north, like Norway (use UTM 33N), or countries not so far south, like New Zealand (use NZTM2000). Another problem with Web Mercator, and other projections centered at the prime meridian, is that you often get troubles crossing the antimeridian on the opposite side. This is unfortunate for New Zealand, where the animeridian is just east of the North Island, leaving the Chatham Islands on the other side. When I tried to reproject my seafloor map into Web Mercator, I got a dataset that made a roundtrip around the Earth. This blog post tells you how to avoid this trip by dividing the source dataset into two parts.

Gisborne, in northeastern New Zealand, proudly claims to be the first city in the world to greet the sun each day.

You’ll see that this seafloor map is different from the one I created earlier, as I’ve masked out the islands with a white color. I did this with Mapnik XML using this LINZ dataset. I want to add land data as separate layers in TileMill.

The seafloor map of New Zealand is 12,000 x 15,200 pixels and it crosses the antimeridian (180° longitude). This is not a problem with the original map projection, as it’s centered at 100° east longitude, and not the prime meridian (0° longitude). What happens if we reproject the dataset into Web Mercator

gdalwarp -s_srs "+proj=merc +lon_0=100 +lat_ts=-41 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -t_srs EPSG:3785 -co compress=lzw -r bilinear nz-seafloor.tif nz-seafloor-webmerc.tif

(To use Mercator 41 with gdalwarp you need to specify the full Proj4-definition.)

The result is a very wide image spanning the globe (and I’m not able to get the same resolution as the original image):

Instead of this wide image, I decided to split the image in two at the antimeridian with gdal_translate:

gdal_translate -srcwin 0 0 7723 15200 -co compress=lzw nz-seafloor.tif nz-seafloor-1.tif

gdal_translate -srcwin 7724 0 4276 15200 -co compress=lzw nz-seafloor.tif nz-seafloor-2.tif

This will create two new map images, one on each side of the antimeridian. There is one column of pixels, which was unfortunate enough to cross the antimerdian, that is now missing.

Next, I’m reprojecting both images into Web Mercator using gdalwarp:

gdalwarp -s_srs "+proj=merc +lon_0=100 +lat_ts=-41 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -t_srs EPSG:3785 -co compress=lzw -r bilinear nz-seafloor-1.tif nz-seafloor-1-webmerc.tif

gdalwarp -s_srs "+proj=merc +lon_0=100 +lat_ts=-41 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -t_srs EPSG:3785 -co compress=lzw -r bilinear nz-seafloor-2.tif nz-seafloor-2-webmerc.tif

Because of the missing pixel, the images don’t overlap entirely at the antimeridian. I’m therefore stretching the images a tiny bit (which is not recognizable at this scale) to make it seamless:

gdal_translate -a_ullr 17482842.051 -7842209.647 20037508.34 -2796127.081 -co compress=lzw nz-seafloor-1-webmerc.tif nz-seafloor-1-webmerc-fix.tif

gdal_translate -a_ullr -20037508.34 -7842313.166 -18622711.547 -2796127.081 -co compress=lzw nz-seafloor-2-webmerc.tif nz-seafloor-2-webmerc-fix.tif

We now have to map images that we can load into TileMill, which will stitch them together again.

Creating a custom map tiling scheme for New Zealand’s seafloor

2012-09-22T13:09:00.000+01:00

First of all I want to thank Brent Wood at National Insitute of Water & Atmospheric Research (NIWA) for inviting me to their office after writing my first blog post about seafloor mapping in New Zealand. Luckily I happened to be in Wellington where they are located, the best city in New Zealand (see photos below)! I was also fortunate to meet Kevin Mackay, the Marine Database Manager who’s responsible for the bathymetry dataset I’ve been using to create this map. You’ll be surprised about the amount of work and data needed to create a bathymetry dataset for such a large and hidden area of our planet. This mapping has been beneficial to New Zealand, as they have been able to claim rights to 1.7 million square kilometers of extended continental seabed, an area about six times New Zealand’s total land area.

In the last blog post we created a seafloor map of New Zealand using shading relief. The map image we created was very large, 12000 x 15200 pixels, with a resolution of 250 meters per pixel. To create an interactive web map we need to cut the image into tiles for a number of zoom levels. As I’m using a custom map projection (Mercator 41), I also have to define a custom tiling scheme. This is my suggestion:

The red area shows the extent of the dataset. The yellow squares show the tiles for zoom level 0 to 2. The biggest square covers an area of 4,096,000 x 4,096,000 meters. Data nerds will know that this number is easily divided by 256, the de facto size of map tiles in pixels. The tile area is defined by these bounds (left, bottom, right, top): 4200000, -6096000, 8296000, -2000000. This table shows the resolution and scale for all zoom levels:

Zoom level 6 is maching the resolution of our dataset, 250 meters per pixel. We now have the information we need to define the tiling scheme in TileCache:

[nz-seafloor]

type=Mapnik

src=EPSG:3994

bbox=4200000,-6096000,8296000,-2000000

maxResolution=16000

mapfile=/home/turban/NZ/Bathymetry/nz-seafloor.xml

I’m not going into detail here, read more about how you can define custom tiling schemes with TileCache, Mapnik and Leaflet in this blog post. I’m using tilecache_seed.py to render the map tiles:

tilecache_seed.py -f nz-seafloor 0 7

Next, we can use Leaflet to load the tiles into an interactive map:

This is what the map looks like:

Show fullscreen

----

A few photos from the best captial city in the world, Wellington:

Creating a seafloor map using shaded relief

2012-09-20T05:54:00.000+01:00

I'll use the bathymetry data from NIWA to create a seafloor map of New Zealand using shaded relief techniques. I used the same techniques when I created the terrain map of Jotunheimen in Norway, so I won't go into detail here.

Let’s start by creating a color relief of the seafloor. I’ve used a blue color gradient similar to the one used by NIWA, progressing from light shallows to dark depths.

Source: NIWA

gdaldem color-relief -co compress=lzw nz-bathymetry.tif depth-gradient.txt nz-bathymetry-relief.tif

I want to enhance the view by using a couple of hillshade techniques. Let’s create a plain hillshade with gdaldem:

gdaldem hillshade -alt 60 -az 135 -co compress=lzw nz-bathymetry.tif nz-bathymetry-hillshade.tif

Notice that the light source (-az) is not from the default top left position (azimuth = 315°), but in a bottom right position (azimuth = 135°). The reason for this is that bathymetry data are inverted compared to elevation data, with higher depths having higher values.

The problem with this hillshade is that it’s all grey, and this grey color will blend with our color relief.
To enhance the hillshade we can apply another technique. This technique is described in the document “Creating Slope-Enhanced Shaded-Relief Using Global Mapper” by Kent D. Brown at Utah Geological Survey. The technique clears out the gray from low-angle slopes so the map colors are not "muddied" by shading. The document describes how to do this in Global Mapper, but we’ll try to do the same the hard way with GDAL and Mapnik.

As described in my previous blog post, creating slope shading is a two-step process. First we need to create a raster from the bathymetry data where each pixel contains an angle, varying from 0 to 90 degrees as the seafloor ranges from horizontal to vertical:

gdaldem slope -co compress=lzw nz-bathymetry.tif nz-bathymetry-slope.tif

Secondly, we use this raster to slope shading by assigning a color to each pixel, using the same technique as for color relief. What we’ll do is to turn all slopes of 0° white, and any slope above 5° black. Anything in between will be in various shades of gray. Adding these two lines to the configuration file will do the trick:

We can then create slope shading with this command:

gdaldem color-relief -co compress=lzw nz-bathymetry-slope.tif slope-ramp.txt nz-bathymetry-slopeshade.tif

Next we can use Mapnik to combine the hill- and slope shade:

nik2img.py nz-bathymetry-hillslopeshade.xml nz-bathymetry-hillslopeshade.tif -d 12000 15200

Lastly, we can combine this enhanced hillshade with the color relief, together with a few datasets from LINZ Data Service:

nik2img.py nz-seafloor.xml nz-seafloor.tif -d 12000 15200 --projected-extent 4800000 -5900000 7800000 -2100000

In my next blog post, I'll turn this image into an interactive web map.

Creating a Seafloor Map of New Zealand

2012-09-19T12:04:00.003+01:00

New Zealand is the perfect country to map, not only because of all the free data available - it also has a very interesting topography, above and below the sea surface. Actually, around 95 percent of New Zealand’s territory is under water. So I wanted to start off my mapping tour by creating a seafloor map of New Zealand.

On my way to Whakaari / White Island, New Zealand’s only active marine volcano. Why is it there and what's hiding below the surface?

You can get bathymetry data (underwater depth) from the National Institute of Water and Atmospheric Research (NIWA). The dataset covers New Zealand's Exclusive Economic Zone and is free to use for non-commercial purposes (terms and conditions). This page describes how the data was collected. I downloaded the “NZ 250m ESRI ASCII grid 530Mb” version of the dataset.

The zip containes two files, RegionalBathymetry2008.asc and RegionalBathymetry2008.prj. The first file contains the actual data in an ASCII grid format. The last file contains information about the projection used, and if you open this file in a text editor it will show you this information:

Projection    MERCATOR
Datum         WGS84
Spheroid      WGS84
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0
Parameters    
 100  0  0.0 /* longitude of central meridian
 -41  0  0.0 /* latitude of true scale
0.0 /* false easting (meters)
0.0 /* false northing (meters)

What does it tell us?

The dataset is in a Mercator projection. This is clearly not the best map projection to use for New Zealand, but the "official" projection, New Zealand Transverse Mercator 2000 (NZTM2000), is only applicable for the main New Zealand island group.
It’s using the WGS84 spheroid.
The coordinate unit is meters.
The central meridian is not the prime meridian (0° longitude), but 100° east longitude. This is useful, as we avoid problems with the antimeridian (180° longitude), which is passing through the area. I’ll explain these problems in a later blog post.
The scale of a map using this projection will only be true at 41° south latitude, running through the Cook Strait between the North and South Islands.

The code for this projection in ESPG:3994 and the Proj4 definition (we'll use this one later) is "+proj=merc +lon_0=100 +lat_ts=-41 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs".

You can use the gdalinfo command to obtain more information about the dataset:

gdalinfo -mm RegionalBathymetry2008.asc

Driver: AAIGrid/Arc/Info ASCII Grid
Files: RegionalBathymetry2008.asc
       RegionalBathymetry2008.prj
Size is 12115, 15391
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["latitude_of_origin",100],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",-41],
    PARAMETER["false_northing",0],
    UNIT["METERS",1]]
Origin = (4795705.551502100200000,-2067845.457881099500000)
Pixel Size = (250.000000000000000,-250.000000000000000)
Corner Coordinates:
Upper Left  ( 4795705.552,-2067845.458)
Lower Left  ( 4795705.552,-5915595.458)
Upper Right ( 7824455.552,-2067845.458)
Lower Right ( 7824455.552,-5915595.458)
Center      ( 6310080.552,-3991720.458)
Band 1 Block=12115x1 Type=Float32, ColorInterp=Undefined
    Computed Min/Max=-3367.677,10069.260
  NoData Value=-9999

In addition to the projection facts, what does it tell us?

The image or grid size is 12,115 x 15,391 pixels or cells.
Each pixel or grid cell covers an area of 250 x 250 meters.
The extent of the area is defined by four coordinate pairs. The size of the area is 3,028,750 x 3,847,750 meters.
The image is single band and the debth data is stored as 32-bit floating-point numbers. The maximum depth is 10,669 meters. Interestingly, the minimum depth is -3,368 meters. This tells us that the dataset also contains onshore elevation data, and the lowest number is actually 3.368 meters above sea level. I’m not sure why this is almost 400 meters below the Mount Cook, the highest mountain in New Zealand (3,754 m).

Let’s add the information we’ve got on Google Earth:

The red area shows the extent of the dataset. The yellow lines show the equator and the central meridian (100° east longitude). The pin is located at the origin of the coordinate reference system (0, 0). The black line is at 41° south latitude, where we have true scale (250 meters per pixel). The resolution will be higher further south, and lower as we move north.

To make the calculations easier to follow, I’m going to crop the dataset a bit with gdal_translate:

gdal_translate -projwin 4800000 -2100000 7800000 -5900000 -a_ullr 4800000 -2100000 7800000 -5900000 -co compress=lzw RegionalBathymetry2008.asc nz-bathymetry.tif

This command will create a dataset that has a size of 12,000 x 15,000 pixels within an area of 3,000,000 x 3,800,000 meters (3,000 x 3,800 km). The resolution is still 250 meters per pixel.

In the nextblog post we’ll start playing with the data.

----

A few more photos from White Island:

Towards the crater.

Between July and August 2012 White Island showed signs of increased activity with lake and gas levels rising from inside the crater. On 5 August 2012 a minor eruption occurred, rising the alert level to 2 of 5.

The main crater. Where's gollum?!?

Walking on the moon. Virtually no vegetation survives the harsh acidic environment.

Numerous failed attempts were made to mine sulphur from White Island. Before the days of antibiotics, sulphur was used in medicines as an antibacterial agent, in the making of match heads, and for sterilising wine corks.

Getting free data in New Zealand

2012-09-11T04:12:00.000+01:00

So far on my travels through New Zealand I’ve been a bit frustrated with the the slow and costly WiFi services. Access to publicly funded data about the country is a totally different story. It’s fast, easy, and free! Countries who lock up their publicly funded data (including my own country, Norway) should definitely look to New Zealand.

So where can you get free data about New Zealand? As I’ve just been in the country for two weeks, I’m not the best person to answer this question. But I’ve done some research and the natives where happy to tell me when I asked this question on my blog a few weeks ago.

Koordinates
This great Auckland-based company has created a portal which allows you to search, sort, bookmark and share geographic data sets. Download in the the format, projection and area you need, connect with web services using industry standard services or developer friendly APIs, or explore data within the browser. It’s a success story in how geodata should be distributed.

LINZ Data Service
Land Information New Zealand (LINZ) is the national provider and administrator of a wide range of geographical information. Here you can download hydrographic and topographic data (up to 1:50,000 scale), in the same easy to use Koordinates-fueled interface.

LRIS Portal
The Land Resource Information Systems (LRIS) Portal is hosted by Landcare Research, an Environmental research organisation specialising in sustainable management of land resources. Here you can find a 25 meter Digital Elevation Model (DEM) of the North and South Island, - and a lot of data about the environment in New Zealand. Some of the data from Lancare Research is also available as a web service (WMS, WMS-C).

The National Land Resource Centre (NLRC)
A science center coordinated by Landcare Research in collaboration with other science organisations, government and businesses. They’ve created an extensive list of datasets provided by different organisations related to land resources.

geodata.govt.nz
A catalogue of publicly-funded geospatial data, and part of New Zealand's developing national spatial data infrastructure (NSDI).

data.govt.nz
A directory of publicly-available, non-personal New Zealand government held datasets. This site does not host data, instead it links to datasets held on other government websites.

DigitalNZ
An inititative led by the National Library of New Zealand helping people to use digital material from libraries, museums, government departments, publicly funded organisations, the private sector, and community groups.

Other resources:

New Zealand Bathymetry (NIWA)
GNS Science
GeoNet - real time updates of earthquakes in NZ

Have I missed something? Please add a comment below and I’ll try to keep this post updated.

---

A few photos from my travels:

At the Tongariro Alpine Crossing together with a group from the great Mountain Cartography Workshop.

The Tongario Alpine Crossing, rated as the best on-day trek available in New Zealand. Part of the track is still closed after the recent eruption at Mount Tongariro.

Mount Ngauruhoe or Mount Doom in Lord of the Rings:

Geothermal activity at Whakarewarewa.

Boiling mud pool created by acid gases and steam that cause the decomposition of minerals to form a clay. The temperature is about 90-95°C.

Pohutu geyser erups once ore twice per hour and can reach heights of up to 30 meters.

(Written at The Bean Café and Roastery in Whakatane, who has the best WiFi offer in New Zealand!)

Mapping through New Zealand

2012-08-11T23:25:00.000+01:00

This autumn I’ll switch hemisphere to experience the spring in New Zealand. I’ll spend two months down under to learn more about cartography and mapping, and to go tramping. Where should I go and what should I see? I need your advice!

My journey starts in Auckland and the GeoCart'2012 conference, before I’m heading up to Tongariro National Park and the 8th Mountain Cartography Workshop. Hopefully, Mount Tongario will stay awake. I’ve never been this close to an active volcano before! Then I’m free to go wherever I want. I roughly plan to stay one month on each island.

What I would like to know:

Are there other mapping events going on in New Zealand in September and October? Anything from conferences, small meetups to short-time courses.
Are there museums or institutions worth visiting for a map geek? I’ve seen that Te Papa in Wellington has an impressive light up satellite photo of New Zealand.
Where can I find geospatial datasets in the public domain? I would like to make my own map of New Zealand while travelling. I know about LINZ and NLRC.
Do you know of an empty office desk or a remote bach (with WIFI!) that I could lend for a week or two. I need to do some work on my Thematic Mapping API.
What will be the best tramps to do on the North Island in September and the South Island in October? I’m a dedicated cross-country skier, but I won’t bring my skis all the way from Norway - so keep it mainly below the snow line.
Other interesting sites and activities.

Please send me a message by adding a comment below, on Twitter, or by sending me an email: bjorn[AT]thematicmapping.org

And, of course - here is a "Where should I go-map"!

Thanks,
Bjørn

Showing GPS tracks with Leaflet

2012-08-08T21:50:00.000+01:00

Last week I had a great 5 day trek through Jotunheimen in Norway, the same area featured on this blog in a number of posts (1, 2, 3, 4, 5, 6, 7) . The trek included the Norway’s highest mountain, Galdhøpiggen (2,469 m), and the deepest valley, Utladalen.

A foggy day on Galdhøpiggen (2,469 m).

Visdalen close to Spiterstulen (1,100 m).

The view from Kyrkja (2,032 m)

Storutla river

Avdalen farm in Utladalen.

I used my Garmin GPSmap 60CSx to track the route, which allows me to save the tracklogs on the memory card. After the trek I had 5 GPX files, one for each day. I used GPSBabel to merge the files into one GPX, jotunheimen_track.gpx.

Leaflet don't have native support for rendering GPX files (there is a plugin though), so I decided to convert the GPX track to GeoJSON. This can easily be done with ogr2ogr program (part of the great GDAL library):

ogr2ogr -f GeoJSON jotunheimen_track.json jotunheimen_track.gpx tracks

I'm only using the tracks from the GPX files, not waypoints and routes. The command generates a file, jotunheimen_track.json, where the complete track is represented as GeoJSON. Leaflet has very good support for GeoJSON. Just assign it to a JavaScript variable and add the following lines of code (using Leaflet 0.4):

L.geoJson(mytrack, {
style: {
color: '#ff0000',
weight: 1,
opacity: 1
}
}).addTo(map);

You can style the track as you like. The resulting map looks like this:

Fullscreen map