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.


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

And here is a look at the same bathy hillshade with coloring

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:

Aerial Imagery
Then merge the models together

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."

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.


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.


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/

I produced a public repository for some of the scripting work. This repository is not specific to the above project but does contain some of the base processing I did on the Wellington elevation data: https://github.com/IReese/wellyvation/tree/master/utils
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:


# 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.

# Create supporting variables

# 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
# Variables, dEm and dSm, are created for naming purposes and moving data to
# the correct directories 

	# 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
		# 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
			namecompare=$( echo $j  | sed "s#.*$capvar##" | sed 's#_.*##' )
			echo $namecompare
			if [ $i = $namecompare ]
				echo $j	>> ${!tempvar_list}/$i.txt	

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

	for k in $listsnames
		# 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 

	# 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


View the FULL BASH SCRIPT HERE via Git Gist.

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.”