Static Vector Tiles for Small Projects

Contours Wellington

Update 04/05/2021: This method has been improved and found here. https://xycarto.com/2021/04/30/static-vector-tiles-ii-openlayers-with-custom-projection/

The following covers a method for visualizing a complex vector tile dataset using a static server.  This particular blog is more for beginners and reduces the technical terms as much as possible. I wouldn’t recommend this method for large complex vector tile data sets, but instead for quick one-off sites that are looking to incorporate a small number of vector tiles where the GeoJSON might be too large/complex to handle.

If you just want to skip to the map and code, visit the Github page created for this blog.

Here is the big reveal at the beginning. If you want to serve a basic vector tile locally or with a static server:

  1. Create an tile cache with uncompressed .pbf files
  2. Place that tile cache on something like GitPages, e.g. embed the tile cache directly along with your website or embed the cache on your localhost along with the other site files
  3. In Leaflet, access the tile cache using ‘L.vectorGrid.protobuf()’

If you’d like to read more in depth:

I work on a Ubuntu OS; however, the methods below will just as easily work with iOS and Windows environments. The parts for building the site, e.g. the HTML, JS, CSS and GitPages, are OS agnostic.

If you want to follow the steps outlined you will need to have installed:

  • GDAL: I am only including this since I used GDAL to extract the contours from the original elevation file.
  • OGR2OGR: Easy tool for converting your shapefile to GeoJSON.
  • Tippecanoe: Excellent tool for generating your tile caches in Web Mercator.
  • T-Rex: Necessary for building your tile cache in a custom projection like NZTM. I tried a number of methods and T-Rex seemed to work best. Sadly, as of Nov 2020, ogr2ogr cannot currently create an MVT in NZTM projection, but it can do projections where the resolutions are halved at each zoom level.
  • Leaflet: Your friendly JS library for building your map in a web environment.
  • SimpleHTTPServer: if you’d like to do localhost testing, you can run your site with this.

You will also need a Github account and a basic knowledge of building a webpage .

Here are the two stacks depending on your needs:

  • Web Mercator: GDAL, OGR2OGR, Tippecanoe, Leaflet, GitPages
  • NZTM: GDAL, OGR2OGR, PostGIS, T-Rex, Leaflet, GitPages

The basic premise is this: with your vector file, you need to construct a tile cache with uncompressed PBF files. The tile cache is the same as any XYZ tile cache you would use for raster tiles, except you are filling it with .pbf files. That tile cache will need to reside in a location accessible by your Leaflet application. In this example we are using GitPages as our server and embedding our tiles directly along with the website we built. Technically, Gitpages is acting as the server. If you are testing and serving on localhost, just embed your tile cache with your web files.

I am going to keep the code light for this blog and instead layout the steps. You can find full code examples in my Github repository here.

Building Your Own Contours

Download Elevation in NZTM
https://data.linz.govt.nz/layer/53621-wellington-lidar-1m-dem-2013/

If you download the entire dataset, you can create the contour lines from the VRT. This saves a lot of time by not creating a mosaic.

Build VRT

gdalbuildvrt dem.vrt *.tif

Contour (@ 50m intervals)

gdal_contour -a elev -i 50 dem.vrt wellyDEMContour.shp

Web Mercator Tile Cache

If you are building a Web Mercator site, you can use Tippecanoe to render the tile cache. You need to create a GeoJSON of your shapefile first

Shape to GeoJSON

ogr2ogr -f GeoJSON -a_srs EPSG:2193 -t_srs EPSG:3857 wellyDEMContour.json wellyDEMContour.shp

Tile JSON

tippecanoe --no-tile-compression --projection=EPSG:3857 --minimum-zoom=12 --maximum-zoom=16 --output-to-directory "static.vector.tiles/contoursWebmer" wellyDEMContour.json

NZTM Tile Cache

If you are building an NZTM site you will need to use T-Rex to generate the NZTM tile cache for the vector tiles. T-Rex likes best if you can give it a PostGIS table to work from. I’d also recommend simplifying complex contour data sets. Your tile generation will be much faster.

Upload to PostgreSQL (you can upload the original NZTM shapefile)

shp2pgsql -s 2193 wellyDEMContour_NZTM.shp public.contournztm | psql -h localhost -d dbName -U userName

T-Rex Config

See full config here

Note the tile size you are setting. You will need this later for your Leaflet application

[grid.user]
width = 4096
height = 4096
extent = { minx = -1000000, miny = 3087000, maxx = 3327000, maxy = 10000000 }
srid = 2193
units = "m"
resolutions = [8960.0,4480.0,2240.0,1120.0,560.0,280.0,140.0,70.0,28.0,14.0,7.0,2.8,1.4,0.7,0.28,0.14,0.07]
origin = "TopLeft"

T-Rex Tile Cache

t_rex generate --progress true --maxzoom=12 --minzoom=9 --extent=160.6,-55.95,-171.2,-25.88 --config /static.vector.tiles/trexConfig/configpsql_contour.toml

Decompress PBF Files

find . -type f | xargs -n1 -P 1 -t -I % gzip -d -r -S .pbf %
find . -type f | xargs -n1 -P 1 -t -I % % %.pbf

Important bit for Leaflet

See the GitPages site for how to set up the Leaflet JS in NZTM


// Access the tile cache
var vector = L.vectorGrid.protobuf(vectorURL, styles);

Have fun!

Open Source Stack for Raster Tiling in Custom Projections

This post is a high level look at the recent stack I built for a raster tiling set up. I am working out some kinks in my online and network delivery of cartographic products, so I thought it was time to set up a raster tiling service to access XYZ and WTMS services from my raster tile caches. I’ll be adding maps and zooms levels in the future, so check back now and again. Antarctica is on it’s way soon!

nz_from_basemap_service

Basic Demo Service using NZTM projection is here: https://xycarto.github.io/

See below for WMTS links

Raster tiling is not the only method, but it is still a viable choice for delivering nice looking maps online, serving across networks, and designing with raster data. I am particularly enamored with the quality of the visual outputs. For me, it is akin to the difference between music in vinyl and digital formats. In addition, the process is well documented and fairly straight forward. By virtue of having been around for a while, raster tiling has a wealth of information and standards to work with, delivery from S3 is a robust process, and there is nice integration with QGIS, Leaflet and Openlayers.

I break the stack in to three areas: analysis, rendering, and delivery

Analysis
QGIS: Sketching, QC, and general geospatial work.

GDAL: Processing raster data. Configuring your rasters in an optimal format from the beginning will greatly improve your rendering speeds. I recommend creating a good set of overviews and gathering everything into a virtual raster tile (VRT).

Postgres/PostGIS: Handling your vector data. Pulling all your data from a database significantly improves rendering speeds. Don’t forget to index!

Rendering
Tilemill/Mapnik XML: Yes, I still design using CartoCSS when working with raster data. I love the simplicity of the language. Tilemill is easy enough to containerize these days too. Tilemill exports into the Mapnik XML format, essential for my process further down the line. Here is how to hack Tilemill to work in a custom projection.

Mapnik: Support for using Mapnik XML

Mapnik with Python Bindings: Necessary for using Mapnik XML documents in MapProxy

MapProxy: MapProxy is a map server and tile renderer . It is easy to build on your machine, though I recommend using a container like Docker. Specifically, I use a hack provided by PalmerJ at Github to increase rendering speeds through multi-threading.

Delivery
Amazon S3: Simple Storage Service. Amazon is pretty cheap, free in many cases, and a good place for storing your tile cache. You get an easily accessed URL for your tiles and a home for your WMTS GetCapabilities document.

WMTS: For me, the real power in a base map service is the WMTS, so, below are two links to the WTMS service for you to set up in QGIS if you’d like to have a play. Here is a quick tutorial about how to set up WMTS if you are unfamiliar.

https://s3-ap-southeast-2.amazonaws.com/basemaps.temp/nz_colour_basemap/WMTSCapabilities.nz_colour_basemap.xml
https://s3-ap-southeast-2.amazonaws.com/basemaps.temp/nz_topo_basemap/WMTSCapabilities.nz_topo_basemap.xml

XYZ: Building a web map? If your tile cache is in S3, in a TMS structure, and public you should be able to access it via simple XYZ request like so:

https://{s3-your-region-here}/{your_bucket}/{project_name}/{projection}/{z}/{x}/{y}.png

Leaflet: Leaflet will handle all the XYZ requests to the server and allow for custom projections. Have a look here for the basic HTML, CSS and JS set up.

Experimenting with Hydrological Analysis using TauDEM

Blue Rivers Ordered

Over the past few years, I’ve played around with developing ordered rivers networks for different projects. I am not an expert in hydrology, but I can get close for cartographic purposes. I am an expert; however, in asking for help from those who know best and I rely on a lot of very smart people to guide me on my journey.

Recently, I decided to put together a visualization of ordered rivers for New Zealand. I came across a very nice data set offered through the Ministry for the Environment via the Koordinates website and thought I’d like to put it to use.

The rivers vis project made me wonder if I could build this base dataset myself using some of the recently released elevation data sets via the LINZ Data Service. The short answer to my question is “sorta”. Doing it open source is not an issue, but building an accurate ordered river centerline network is another story. This is a task I cannot take on as a solo project right now, but I could do a little experimentation. Below, I’ll offer some of methods and things I learned along the way.

Tools and Data

The method I tested used TauDEM and a 1m DEM raster accessed from the LINZ Data Service. I down sampled the DEM to 2m and 5m resolutions and used small areas for testing. Finding and open source tool was easy. I sorted through a few available methods and finally landed on “Terrain Analysis Using Digital Elevation Models” (TauDEM). There are additional methods through GRASS and SAGA GIS. I chose TauDEM because I never used it before.

Method Tested

To my knowledge, there is no open source tool where a person can put in a DEM and get a networked rivers centerline vector out the other side. It requires a number of steps to achieve your goal.

The basic run down to process the DEM is to:

  1. Fill sinks
  2. Determine flow directions
  3. Determine watersheds
  4. Determine flow accumulation
  5. Stream classification
  6. Export to vector

TauDEM does require a few extra steps to complete the process, but these steps are explained in the documentation of the tool. It was more about keeping all my variables in the right places and using them at the right time. I recommend using the variable names TauDEM provides.

BASH script here

Click the arrow to the left to view the full BASH script below:

#!bin/bash

#Rough sketch for building river centerlines. Rasters have been clipped prior

BASEPATH=/dir/path/to/base

raster_list=$( find $BASEPATH -name "*.tif" )

taudem_outputs=/dir/path/to/outputs

reso=resolution_number

for i in $raster_list
do

	INPUT_RASTER=$i

	file_name=$( basename $i )

	strip_input_extension=$( echo $file_name | sed 's/.tif//' )

	reso_name=$taudem_outputs/${strip_input_extension}_${reso}res

	gdal_translate -tr $reso $reso -of GTiff $i $reso_name.tif

	fel=${reso_name}_fel.tif
	p=${reso_name}_p.tif
	sd8=${reso_name}_sd8.tif
	ad8=${reso_name}_ad8.tif
	ang=${reso_name}_ang.tif
	slp=${reso_name}_slp.tif
	sca=${reso_name}_sca.tif
	sa=${reso_name}_sa.tif
	ssa=${reso_name}_ssa.tif
	src=${reso_name}_src.tif

	ord=${reso_name}_strahlerorder.tif 
	tree=${reso_name}_tree.dat
	coord=${reso_name}_coord.dat
	net=${reso_name}_network.shp
	w=${reso_name}_watershed.tif 

	processed_input_file=$reso_name.tif

	#TauDEM Commands
	mpiexec -n 8 pitremove -z $processed_input_file -fel $fel

	mpiexec -n 8 d8flowdir -fel $fel -p $p -sd8 $sd8 

	mpiexec -n 8 aread8 -p $p -ad8 $ad8 -nc

	mpiexec -n 8 dinfflowdir -fel $fel -ang $ang -slp $slp

	mpiexec -n 8 areadinf -ang $ang -sca $sca -nc

	mpiexec -n 8 slopearea -slp $slp -sca $sca -sa $sa

	mpiexec -n 8 d8flowpathextremeup -p $p -sa $sa -ssa $ssa -nc

	mpiexec -n 8 threshold -ssa $ssa -src $src

	mpiexec -n 8 streamnet -fel $fel -p $p -ad8 $ad8 -src $src -ord $ord -tree $tree -coord $coord -net $net -w $w

done

The script is a rough sketch, but does get results.

Challenges in the Process

One major challenge for this project was the size of the input DEM and my computers available RAM. I work primarily off a laptop. It’s a good machine but no match for a proper server set up with some spacious RAM. My laptop struggled with the large hi-resolution DEMs, so I needed to down-sample the images and choose a smaller test area to get it to work.

Clip the tiff with gdal_translate -projwin and down sample with -tr

gdal_translate -tr xres yres -projwin ulx uly lrx lry input.tif output.tif


The second challenge came up because I used a bounding box to clip my test regions. I recommend not doing this and instead clip your regions using a watershed boundary. Having square shapes for your test regions will give you very inaccurate and unhelpful results. For example, major channels in your DEM will be cut at the edges of your raster. You will not get accurate results.

Clipping a raster using a shapefile, like a watershed boundary, can be achieved using gdalwarp.

gdalwarp –cutline input.shp input.tif output.tif


Results

I ran my process and QCed the results against Aerial Imagery and a hillshade I developed from the DEM. The first run gave me good enough results to know I have a lot of work to do, but I did manage to develop a process I was happy with. The tool did a great job, but the accuracy of the DEM was a little more challenging. It’s a start. I captured a good number of river channels despite my incorrect usage of a square DEM, learned a lot about how DEM resolution affects outputs, and gained knowledge around how to spot troublesome artifacts.

Well Defined ChannelsImg 1: River capture in well defined channel.

From this experiment, there are a few ideas I’d like to explore further:

1. Accuracy of the DEM. The particular DEM I worked with had a number of ‘dams’ in the flows. Notably, bridges, culverts, vegetation artifacts, and other general errors that caused water to flow in interesting directions. When working with a data set like this, I am curious how manage these artifacts.

Road issueImg 1: River diversion at road.

Artifact issueImg 1: River diversion at culvert or bridge.

2. How to go beyond borders. This analysis can be broken down by watershed, but it will be necessary to link the outflows of those watersheds to the next for accurate results.

Edge issueImg 1: Flow not captured at edge.

3. As DEMs are released with better resolution, there is a need for scaled up computing power. The process needs a large amount of RAM. What is the best computational set up for capturing the largest area?

4. Did I do this correctly? I perform this task about once every two years and usually weekends when the surf is flat and the garden is weeded, so I am not an expert. There is a lot more research to be done to determine if I am using the tools to the best of their abilities.

Wellington Elevations: Interpolating the Bathymetry

It is important to note something from the very beginning. The interpolated bathymetry developed in this project does not reflect the actual bathymetry of the Wellington Harbour. It is my best guess based on the tools I had and the data I worked with. Furthermore, this interpolation is NOT the official product of any institution. It is an interpolation created by me only for the purposes of visualization.

welly_harbour-colour-and-aerial_FULLVIEW

Part of the goal when visualizing the Wellington landscape was to incorporate a better idea about what may be happening below the surface of the harbor. Various bathymetric scans in the past have gathered much of the information and institutions like NIWA have done the work visualizing that data. As for myself, I did not have access to those bathymetries; however, I did have a sounding point data set to work with, so I set about interpolating those points.

The data set, in CSV format, was over a million points; too dense for a single interpolation. I worked out a basic plan for the interpolation based on splitting the points into a grid, interpolate the smaller bits, then reassemble the grid tiles into a uniform bathymetry.

Conversion from CSV to shp
Using the open option (-oo) switch, OGR will convert CSV to shp seamlessly

ogr2ogr -s_srs EPSG:4167 -t_srs EPSG:4167 -oo X_POSSIBLE_NAMES=$xname* -oo Y_POSSIBLE_NAMES=$yname*  -f "ESRI Shapefile" $outputshapepath/$basenme.shp $i

Gridding the Shapefile
With the shapefile in place, I next needed to break it into smaller pieces for interpolation. For now, I create the grid by hand in QGIS using the ‘Create Grid’ function. This is found under Vector>Reasearch Tools>Create Grid. Determining a grid size that works best for the interpolation is a bit of trial and error. You want the largest size your interpolation can manage without crashing. Using the grid tool from QGIS in very convenient, in that it creates an attribute table of the xmin, xmax, ymin, ymax corrodinates for each tile in the grid. These attributes become very helpful during the interpolation process.

Interpolating the Points
I switched things up in the interpolation methods this time and tried out SAGA GIS. I have been looking for a while now for a fast and efficient method of interpolation that I could easily build into a scripted process. SAGA seemed like a good tool for this. The only drawback, I had a very hard time finding examples online about how to use this tool. My work around to was to test the tool in QGIS first. I noticed when the command would run, QGIS saved the last command in a log file. I found that log, copied out the command line function, and began to build my SAGA command for my script from there.

Here is look at the command I used:


saga_cmd grid_spline "Multilevel B-Spline Interpolation" -TARGET_DEFINITION 0 -SHAPES "$inputpoints" -FIELD "depth" -METHOD 0 -EPSILON 0.0001 -TARGET_USER_XMIN $xmin -TARGET_USER_XMAX $xmax -TARGET_USER_YMIN $ymin -TARGET_USER_YMAX $ymax -TARGET_USER_SIZE $reso -TARGET_USER_FITS 0 -TARGET_OUT_GRID "$rasteroutput/sdat/spline_${i}"


I tested a number of methods and landed on ‘grid_spline’ as producing the best results for the project. It was useful because it did a smooth interpolation across the large ‘nodata’ spaces.

Once the initial interpolation was complete, I needed to convert the output to GeoTIFF since SAGA exports in an .sdat format. Easy enough since GDAL_TRANSLATE recognizes the .sdat format. I then did my standard prepping and formatting for visualization:


gdal_translate "$iupput_sdat/IDW_${i}.sdat" "$output_tif/IDW_${i}.tif"
gdaldem hillshade -multidirectional -compute_edges "$output_tif/IDW_${i}.tif" "$ouput_hs/IDW_${i}.tif"
gdaladdo -ro "$output_tif/IDW_${i}.tif" 2 4 8 16 32 64 128
gdaladdo -ro "$ouput_hs/IDW_${i}.tif"2 4 8 16 32 64 128


Here is look at the interpolated harbour bathymetry, hillshaded, with Wellington 1m DEM hillshade added over top
welly_harbour_bw_all

And here is a look at the same bathy hillshade with coloring
welly_harbour_bw-and-aerial

Visualizing the Bathymetry
With the bathymetry, complete it was simply a matter of building it into the existing visualization I built for the Wellington Region. Learn more about the project here. The visualization was four steps:

Hillshade
addedbathy_bathyonlypng
Color
addedbathy_bathyonly_withcolor
Aerial Imagery
addedbathy_bathyonly_withcoloraerial
Then merge the models together
addedbathy_final

Easy as, eh? Let me know what you think!

Note: All imagery was produced during my time at Land Information New Zealand. Imagery licensing can be found here:
“Source: Land Information New Zealand (LINZ) and licensed by LINZ for re-use under the Creative Commons Attribution 4.0 International licence.”

Building the Wellington Model with 1m DEM and DSM

As interest in LiDAR derived elevation increases, so grows the interest in the capabilities. LiDAR derived elevation data has been great for my visualization game and in helping me communicate the story out about what LiDAR can do. It all starts with a picture to get the imagination going.

wellyvation

The Wellington model derived for this project is part of an ongoing project to help increase the exposure of the Wellington 1m DEM/DSM elevation data derived from LiDAR. Step one for me is getting a working model built in QGIS, capturing still images, and increasing interest in the data.

I’ve talked about the processing of the elevation data for Wellington visualizations in the past, so for this post I’m only focusing on the blending of the data sets in building the model. This project is a good model since it encompasses a number of subtle techniques to get the model to stand out. This post is one of a two part series; the second post discusses the techniques used to derive and visualize the bathymetry for the surrounding harbor.

Let’s start with the base, Aerial Imagery.
wellyhabour_aerialonly

Blended with a hillshade
wellyhabour_aerial_withHS

DSM added for texture and context
wellyhabour_aerial_withDSMHS

Slope added to define some edges
wellyhabour_aerial_withDSMDEMSLOPEHS

Some darker shading added to the bathymetry to frame the elevation data
wellyhabour_aerial_withDSMDEMSLOPEHS_darkenframe

And finally some added bathymetry to lighten the edges at the shoreline enhancing the frame a bit more.
wellyhabour_aerial_withDSMDEMSLOPEHS_edgeframe

In the end there is some post-processing in Photoshop to lighten up the image. Honestly, this could have been done in QGIS, but I was being lazy. For the images produced, there was no need to retain the georeferencing, and when that is the case, I rely on Photoshop for color and light balancing.

The greatest difficultly in this project so far has been trying to create a universal model for the data set. I’m finding that as I visualize different regions using this model, I need to adjust the hillshading quite significantly to draw out different features. Take a look at the images here. It is the same model, but with the noticeably different gradients used in the hillshades. The techniques used for the images in this post worked well for the urban region shown, but fall apart as you move further out into the more mountainous regions. Much of the blending is too harsh and turns the mountains into a black muddled mess. I am almost there, but like any project, it takes a good bit of subtle tweaking of the blending to get a universal image to work.

The entire base mapping work is completed in QGIS. The elevation data was processed using GDAL and the bathymetric interpolations were produced SAGA GIS. There are no color palettes for this project. The aerial imagery does all the work in that department.

Base data can be found here:
DEM: https://data.linz.govt.nz/layer/53621-wellington-lidar-1m-dem-2013/
DSM: https://data.linz.govt.nz/layer/53592-wellington-lidar-1m-dsm-2013/
Aerial Imagery: https://data.linz.govt.nz/layer/51870-wellington-03m-rural-aerial-photos-2012-2013/

The next post covers the development of the bathymetry for the surrounding harbor. Thanks for having a look and let me know what you think.

Note: All imagery was produced during my time at Land Information New Zealand. Imagery licensing can be found here:
“Source: Land Information New Zealand (LINZ) and licensed by LINZ for re-use under the Creative Commons Attribution 4.0 International licence.”

The Rejects

Sometimes there is simply not enough room for all the ideas. Sometimes you want all the images to make it to the final round.

wairarapa

In a recent project to promote some of our elevation data, I was asked to present a number of ideas for a 2000mm x 900mm wall hanging. The piece was to act as a conversation starter and demonstrate some of the finer details elevation from LiDAR possesses.

In the end, the image above was the chosen candidate. Below are the drafts I initially presented for review. You can see the difference in treatment from the original ideas to the final product. Personally, I really enjoyed the images developed for the draft series, I liked the silvery undertones, and I thought it was a shame to merely let these images sit on my hard drive.
Below, you’ll find a brief description about a few challenges faced in the image development.

near_lake_ferry
nice_farm
masterton_region
random
draft_wairarapa

Artifacts and Finer Details
The hardest part of this job was drawing out the finer details of the chosen location. There was a strong interest in showing the ancient river bed; however, without a good bit of tweaking in the hillshades, the image is quite flat. After some trial and error, I found I could get a good contrast by limiting the hillshade values range to 170-190. That’s it, but the readability of the project really hinged on the simple tweak. It really made the details stand out.
That said, the gain in detail also revealed a significant artifact in the data. If you go back up and have a closer look, you will find diagonal depressions running across the images in equal intervals. These are lines from where the LiDAR scans overlap. I haven’t quite had the time to figure out how to remove these from the original data source, so for now I leave them in as conversational piece around improving LiDAR capture practices.
As usual, all map layout work was completed on QGIS, with the bulk of the data processing done using GDAL. The ‘Reject’ images for this post are direct exports from QGIS, with no manipulation apart from some down-sampling and cropping in Photoshop.

Base data can be found here:
DEM: https://data.linz.govt.nz/layer/53621-wellington-lidar-1m-dem-2013/
DSM: https://data.linz.govt.nz/layer/53592-wellington-lidar-1m-dsm-2013/
Aerial Imagery: https://data.linz.govt.nz/layer/51870-wellington-03m-rural-aerial-photos-2012-2013/

Hope you like and thanks for checking in!

Note: All imagery was produced during my time at Land Information New Zealand. Imagery licensing can be found here:
“Source: Land Information New Zealand (LINZ) and licensed by LINZ for re-use under the Creative Commons Attribution 4.0 International licence.”

Processing and Visualizing Auckland 1m DEM/DSM Elevation Data

About two years ago, I took on a cartographic project visualizing the Auckland 1m DEM and DSM found publicly via the LINZ Data Service (LDS) here: DEM, DSM. The goal at the time was to develop a base map for the extraction of high resolution images for use in various static media. It was a good piece of work with some fun challenges in scripting and gradient development. Included herein are notes about processing the data with QGIS and BASH, building the gradients, and blending the base maps using QGIS.
two_volcanoes_forWebImg 1: Conference media developed for International Cartography Conference (ICC) 2018, Washington DC.
Processing the Data
The original data download was 7GB per data set (DEM and DSM). Each data set contained 6423 individual files at 1.4MB.  This set up was pretty hard to work with in QGIS; so, initially I processed each data set for ease of viewing.  This processing included grouping the data into larger files matching the LINZ Topo 50 Map Grid Sheets and running a few processes like GDALADDO (overviews), GDALDEM (hillshades), and GDALBUILDVRT (virtual mosaic).

The basic idea of formatting the data is as follows:

gdaldem hillshade -multidirectional -compute_edges input.tif output.tif
gdaladdo -ro input.tif 2 4 8 16 32 64 128
gdalbuildvrt outputvrt.vrt *.tif
Click the arrow to the left to view the full BASH script below:

#!bin/bash

# The purpose of this script is to process the Auckland 1m DEM and DSM elevation
# data into more manageable pieces for easier viewing in QGIS...  
#... The original
# elevation tile downloads from LDS contain 6423 individual tiles. The
# downloaded elevation tiles are reworked into tiffs the same size as the NZ
# LINZ Topo50 Map Sheets 
#(https://data.linz.govt.nz/layer/50295-nz-linz-map-sheets-topo-150k/).  In
# this case, the original data contains an identifier, like 'AZ31', within 
# the tile name that associates it with Topo50 Map Sheets.  This script 
# extracts that identifier, makes a list of the files containing the identifier
# name, makes a vrt of the items in the list, creates hillshades from that vrt,
# then formats for quicker viewing in QGIS.

# All data is downloaded in EPSG:2193 and in GeoTiff format
# Auckland DEM here: https://data.linz.govt.nz/layer/53405-auckland-lidar-1m-dem-2013/
# Auckland DSM here: https://data.linz.govt.nz/layer/53406-auckland-lidar-1m-dsm-2013/

# Place ZIPPED files in directory of choice

# Set root directory for project. PLACE YOUR OWN DIRECTORY HERE.
BASEDIR=[PLACE/YOUR/OWN/BASE/DIRECTORY/FILEPATH/HERE]


# Create supporting variables
dSm_dir=$BASEDIR/dSm_elevation
dEm_dir=$BASEDIR/dEm_elevation
dSm_list_dir=$BASEDIR/lists/dSmlist
dEm_list_dir=$BASEDIR/lists/dEmlist


# Create file structure
mkdir $BASEDIR/lists
mkdir $dSm_dir
mkdir $dEm_dir
mkdir $dSm_list_dir
mkdir $dEm_list_dir

# Extract data
unzip $BASEDIR/lds-auckland-lidar-1m-dsm-2013-GTiff.zip -d $dSm_dir
unzip $BASEDIR/lds-auckland-lidar-1m-dem-2013-GTiff.zip -d $dEm_dir

# Delete zipped files
# rm -rf $BASEDIR/lds-auckland-lidar-1m-dsm-2013-GTiff.zip
# rm -rf $BASEDIR/lds-auckland-lidar-1m-dem-2013-GTiff.zip

# Loop to process both DEM and DSM data
demdsm="dEm dSm"
for opt in $demdsm
do
# Variables, dEm and dSm, are created for naming purposes and moving data to
# the correct directories 
tempvar=""$opt"_dir"
tempvar_list=""$opt"_list_dir"
capvar="${opt^^}_"

	# Identify associated Topo50 map sheet name.  Make it as a list held as a variable
	unique=$( find ${!tempvar} -name "*.tif" | sed "s#.*$capvar##" | sed 's#_.*##' | sort | uniq )

	# from the 'unique' variable, create a list of files with similar Topo50 idenifier
	for i in $unique
	do
		# List all available tiffs in directory 
		namelist=$( find ${!tempvar} -name "*.tif" -maxdepth 1 )
		# Compare unique name to identifier in available tiffs name.  If 
		# there is a match between the unique name and identifier in the 
		# tiff name, the name is recorded in a list.
		for j in $namelist
		do
			namecompare=$( echo $j  | sed "s#.*$capvar##" | sed 's#_.*##' )
			echo $namecompare
			if [ $i = $namecompare ]
			then
				echo $j	>> ${!tempvar_list}/$i.txt	
			fi
		done
	done

	# Create list of available .txt file 
	listsnames=$( find ${!tempvar_list} -name "*.txt" )

	for k in $listsnames
	do
		# list contents of .txt file into variable
		formerge=$( cat $k )
		# prepare file name to use as vrt name
		filename=$( basename $k | sed 's#.txt##' )
		#echo $filename
		#echo $formerge
		# Build VRT of elevation files in same size as Topo50 grid 
		gdalbuildvrt ${!tempvar}/$filename.vrt $formerge 
	done

	# Change directory to 'Merged DEMs'
	cd ${!tempvar}

	# Make directory to store hillshade files
	mkdir hs

	# Clean out overviews
	find -name "*.vrt" | xargs -P 4 -n4 -t -I % gdaladdo % -clean

	# Create hillshade from VRTs
	find -name "*.vrt"  | xargs -P 4 -n4 -t -I % gdaldem hillshade -multidirectional -compute_edges % hs/%.tif

	# Create external overviews of VRTs
	find -name "*.vrt" | xargs -P 4 -n4 -t -I % gdaladdo -ro % 2 4 8 16 32 64 128

	# Create vrt of elevation VRTs
	gdalbuildvrt $opt.vrt *.vrt

	# change directory to hillshade directory
	cd ${!tempvar}/hs

	rename s#.vrt## *.tif

	# Clean out old overviews
	find -name "*.tif" | xargs -P 4 -n4 -t -I % gdaladdo % -clean

	# Create external overviews of HS tiffs
	find -name "*.tif" | xargs -P 4 -n4 -t -I % gdaladdo -ro % 2 4 8 16 32 64 128

	# Create vrt of Hillshade tiffs
	gdalbuildvrt "$opt"_hs.vrt *.tif

done

Building the Gradients
Getting a natural transition through the land and sea was difficult. I was presented with two challenges; 1. building a colour gradient for elevations spanning mountain tops to undersea and 2. determining which intertidal feature to model. 

Studying Aerial Imagery from the region, I determined three zones I’d develop gradients for:

  1. Bathymetric
  2. Intertidal
  3. Terrestrial

With the gradient zones in place, I developed a few additional rules to keep the project linked visually.

  1. The colours for each gradient would be linked in tone, but distinct from each other.
  2. The bathymetry colour would frame the intertidal and terrestrial data.
  3. The deep sea blue would anchor the colour pallet.

Finally, I needed an intertidal model. Since, there are a number of different intertidal zones represented around Auckland: estuaries, sandy beaches and rocky shores and I am using only one DEM, I needed to choose which intertidal zone to represent in the image.  One colour gradient does not fit all. I had to make a choice.  I decided to use the mud flats and estuaries as the primary intertidal model.  The DEM included the channels in the mudflats and I really liked the shapes they made.  They also covered wide areas and were a prominent feature.elevation_cross_sectionImg 2: Elevation model focusing on marshy tidal zones

With the intertidal model determined and the zones set, I developed the colour gradients.

Elevation Value Colour Zone
-1.0 Bathymetric
-0.75 Bathymetric
0.5 Intertidal
1.0 Intertidal
1.7 Intertidal
1.8 Terrestrial
25 Terrestrial
100 Terrestrial
500 Terrestrial

Note: The elevation values used do not necessarily correlate with the actual elevations where these zones transition. They are a best estimation based on samplings from aerial imagery.

Blending the Layers
For the final step, I blended the layers in QGIS.  I needed the hillshades I developed from the DEM and DSM, plus the original DEM elevation. That’s it.  Three layers for the whole thing. The DEM elevation carried all the colour work, the DSM hillshades gave the detail, and the DEM hillshade added some weight to the shaded areas. Here is the order and blending for the project in QGIS:

  1. DSM Hillshade: multiply, brightness 50%, black in hillshade set to #333333
  2. DEM Hillshade: multiply, brightness 50%
  3. DEM: contrast 10%

Overall, the base image proved to be a success and the script has been useful across a number of projects. The images have ended up in a good bit of internal and conference media and I have seen steady use for almost two years now. For me, the image is getting tired and I’d eventually love to redevelop the gradients; but, I am happy to have a chance to write about it and get some more external exposure. In the future, I am looking to develop this map a bit further and present it as a web map as well. Time will tell whether this happens or not. I think it would take a directive from an outside source.

I’d be keen to hear your comments below or get in touch if you are interested in learning more.

rando_forWebImg 3: Promotional media

Note: All imagery was produced during my time at Land Information New Zealand. Imagery licensing can be found here:
“Source: Land Information New Zealand (LINZ) and licensed by LINZ for re-use under the Creative Commons Attribution 4.0 International licence.”