Creating a shaded relief map of New Zealand

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: