The following is a method for converting QGIS projects into stylized Cloud Optimised Geotiffs (COG)s. The process is devised as an alternative to raster tiling and explores the use of COGs as a raster-tile-esque format along with QGIS as a styling editor replacing Tilemill.

For those unfamiliar with the COG format, there is an excellent explanation here.  

For COG uses in QGIS, users readers can look here.

Basic Method

  1. Develop QGIS project
  2. Build in zoom rules for scales
  3. “Export as Image” based on scale and zoom resolution
  4. Images converted to overview structure, with the lowest scale being the base of the overviews
  5. Combine images into COG.
  6. Upload COG to S3 and consume via webviewer or QGIS

The Method in Detail

The method relies on two processes readily available today. One is the QGIS “Export to Image” function and two, is a quirk of GDAL overviews, where an overview can read an overview of itself. 

Let’s build a hypothetical QGIS project with a hillshade.tif, multiplied into an elevation DEM.tif gradient, and a vector roads layer placed over top. With the roads layer in place, build a few varying road widths like:

  • road_width = 1 at 250k scale
  • road_width = 5 at 100k scale
  • road_width = 10 at 50k scale

With these rules in place, let’s next export a few images from capturing the style changes using the QGIS “Export to Image” function. ‘Export to Image” allows users to export images from the project based on extent, resolution, and scale. Sort of like taking snapshots of a project at different zoom levels and in different areas of interest.  When used properly, the exported image will honour the zoom rule for that scale.  So, from the example, users could export three images using the scale capturing the road width changes. Let’s say the exports are named:

  • 50k.tif
  • 100k.tif
  • 250k.tif

Great, we can export images with zoom rules honoured in the images, but what does it matter? Well, if we could “stack” these images in QGIS, have the resolution set to half of each other,  and have them automatically turn on/off as users zoom, we would see the changes happening for each scale. This is entirely possible and how overviews work in geospatial;  we now just need to hack the overview method for the purposes of viewing our images in this manner. 

Here is the basic premise: in order to get our images to act as overviews, we want to reduce the geospatial image resolution (not DPI) by half of the image before it. So, if the image at the bottom (50k.tif) is 14, then the 100k.tif is 28.  These are strange DPIs, but will make more sense in a moment. 

The resolution we are seeking for each image is based on the screen resolution, scale, and pixel width of the tile. In its most basic form, the resolution of the for a given scale can be found by multiplying the scale (in metres) by 0.00028.

  • Resolution = Scale * 0.00028
  • Where, 0.00028 is metres/pixel.  

Knowing this equation, we can begin to build a “tile matrix” composed of Zoom Level, Scale, and Resolution.   Here is good explanation of the tile matrix from the OGC.

Matrices are available for a number of projections.  The NZTM tile matrix developed by LINZ is here. Inspecting the NZTM matrix,  the resolutions we are seeking for the image scales are:

Snippet of NZTM Tile Matrix for the Web

Now we know the image resolution we need,  when we “Export to Image” from QGIS we can set the proper resolution for the images as if they were stacked into a pyramid. Let’s revise the export of the images and use the correct resolution:

  • 50k.tif (image resolution = 14)
  • 100k.tif (image resolution = 28)
  • 250k.tif (image resolution = 70)

I do this programmatically using pyQGIS leveraging QPainter.  You can see a full python script here.

Here is where the Overviews now come into play. We want the 50k.tif image to see the 100k.tif image as an overview and so on down the line.  Fortunately, we can rely on a strange quirk of GDAL to do this.  Using GDAL, or any other raster tool, convert your image stack look like so:

  • 50k.tif -> No change
  • 100k.tif -> 50k.tif.ovr
  • 250k.tif -> 50k.tif.ovr.ovr

GDAL command line example: 

gdal_translate \
	100k.tif \
        50k.tif.ovr \
        -of "GTiff" \
        -r bilinear \
        -co PROFILE=BASELINE \
        -co BIGTIFF=YES \
        -co TILED=YES

* Note: we are stripping the header from only the overview images using PROFILE=BASELINE. We don’t need the additional information and it makes a lighter file.

Now we have:

  • 50k.tif 
  • 50k.tif.ovr
  • 50k.tif.ovr.ovr

You could now load only the 50k.tif file into QGIS and see the changes from the other files happen as you zoom in. If you wanted and you only planned to use QGIS, you could leave it here, but this post is about how to bring this into a single file format (COG), ready for the web.

The last step is to bring all these files into a single file format, the COG. When creating a COG, GDAL will allow users to use external overviews we already generated instead of generating them itself.  That is what we do here.  We point to the single base file, set the output format to COG, and compress the tiled overviews.

gdal_translate \
	50k.tif \
	50k-cog.tif \
	-of COG \

In the end, we have a single COG tif file, with internal overviews tiled and acting like “zoom” scales.

From this point, the COG can be transferred to a location like S3 and accessed by Openlayers. You can view loading a three band COG in Openlayers here, to get a sense on how to access the file.  This is how I load the COG using Openlayers:

const url = "https://d3cywq4ybqu7io.cloudfront.net/cogs/as-raster-tile/50000-cog.tif"

const cogSource = new GeoTIFF({
  sources: [
  convertToRGB: true,

const cog = new TileLayer({
  crossOrigin: 'anonymous',
  source: cogSource,
  extent: extent,

* the notable bit here is that you might need to use: convertToRGB: true

You can view the website with the COG in action here.

The COG may also be directly access from S3 in QGIS without the need to download the file. In the “Add Layers”, you can access the COG as a “Raster” using the “http” connection.


There are a lot of moving parts in this process and it can get a bit confusing about what step goes where. Here is a link to the Git repo I built with a few examples. Please do contact me with questions and I will try to help however I can.

Basic COG in Openlayers:  Single Band Tif

This post covers loading a raw COG Tif and manipulating the values in JS. If you are unfamiliar with the COG Tif format, see here for an explanation. Openlayers has a few good examples on how to load COGs.  Some of this is a repeat of their examples and some goes a little more in depth. 

  • Example site is running here
  • Github repo is here

The following will cover:

  1. Loading a COG via Openlayers 6
  2. Demonstration of a pop-up to query the COG value

Data in the site:

  1. Sea Surface Temperature (SST) from JPL. Clipped to New Zealand EEZ
  2. LINZ Aerial Imagery Basemap. The link is temporary.  I highly suggest you get your own link from here.  

Few Notes:

  1. The COG and the website are built using WEB MERCATOR projection, so there are no special steps to build the site with regards to projection 
  2. I am only covering the JS in this example, see here for a full HTML/CSS/JS example, you can look here

Creating the COG

There are many tutorials out there on creating a COG.  I am not going to repeat these.  Instead, I am providing a link to a COG I created. The COG used in this example was created using a standard python build like:

creation_options = [

    format = "COG",
    creationOptions = creation_options

I put the COG on S3 and made the link public


Loading a raw COG tif

First, set the url to the COG and load it as a source:

var urls3 = 'https://d3cywq4ybqu7io.cloudfront.net/cogs/sst/clipped-eez-nztm-20200101090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02_fill_cut_warp_cog.tif'

var cogSource = new ol.source.GeoTIFF({
    normalize: false,
    sources: [
        url: urls3,
        min: 277,
        max: 300,
        nodata: -32768,

Here, we set the min/max of the value from the COG.

Second, we build the colour gradient for the COG.  This will tell the client how to colour the values of the COG:

var cogBand = ['band', 1]

var defaultColor = {
color: [
    276.9, [255, 255, 255, 0],
    277, [19, 22, 180, 1],
    284, [70, 111, 207, 1],
    289, [196, 229, 183, 1],
    294, [217, 164, 73, 1],
    300, [199, 69, 40, 1]

Here, we are using a linear interpolation to create a smooth gradient between the values.

Third, load the COG with colorization:

var cog = new ol.layer.WebGLTile({
    visible: false,
    crossOrigin: 'anonymous',
    source: cogSource,
    style: defaultColor,

Here, the visibility is set to ‘false’.  When the map initializes on the browser, I want it turned off at the start to allow users to toggle the layer on/off.

Setting up the query button

This is a standard pop-up button from Openlayers.  It is developed so users can query the SST map and see the actual value at that location. You’ll need to do some work in the HTML and CSS. See here for more detail in the code: 

First, build the pop-up in the JS at the top of the file:

// Pop up set
var container = document.getElementById('popup');
var content = document.getElementById('popup-content');
var closer = document.getElementById('popup-closer');

var overlay = new ol.Overlay({
    element: container,
    autoPan: {
      animation: {
        duration: 250,

closer.onclick = function () {
return false;

Second, set the pop-up to query the COG layer:

// Set onclick to return values from COG
map.on('singleclick', function(evt) {
    var coordinate = evt.coordinate;
    var data = cog.getData(evt.pixel);
    var celcius = data[0] - 273.15
    var codeText = "Temp in Celcius"
    content.innerHTML = "<div class='popupText'>Sea Surface Temperature: <strong>" +     
        celcius.toFixed(2) + "</strong><div class=returnVal>" + codeText + "</div></div>";

Those are the main bits. The remainder of the code in the JS is setting up a base map, building the toggle button, and preparing the query pop-up.

Of course, you can clone the repo and do what you like to make it your own special site.

As usual, this example is only one way to do this. There are many other setups and special circumstances that will negate what is shown here.

Let me know if this this helped or if you have questions. Also, if you used this, I’d love to see what you built!

Basic GRASS GIS with BASH, plus GDAL

As a follow-up to the last blog, I thought it would be helpful to demonstrate how next to break up the elevation example into individual watersheds. The reason being, in my last example I demonstrated the process using a square raster tile. Hydrological processes are not accurate when run on square tiles. It is best to run further processes, like stream order extraction, on complete watersheds.

If you are reading this and want to follow along, it is very helpful to go back and read this blog. The following builds off that last example.

I’ll break this into two parts:

1. Creating multiple watershed boundaries of different sizes with GRASS and using a basic loop in BASH for the process.
2. Clipping the original raster by the watershed boundaries using GDAL and SQL with a basic loop in BASH.

Creating Watershed Boundaries with Varying Sizes

This section will modify the script created in the last post.

Notice in the previous script, when I used r.watershed, I set the threshold value to 100000. You can use this value to set the minimum size of your watershed based on cell units. I see this as minimum number of pixels per watershed since your cell unit is the size of your pixel.

r.watershed elevation=$fillDEM threshold=$threshold accumulation=$accumulation drainage=$drainage stream=$stream basin=$basin --overwrite

Have a play setting the threshold value to different values, like 2000000 or 5000. You will get very different outputs in your watersheds vector layer. Try to determine the upper most threshold value you can input before GRASS can no longer determine a watershed.

Warning: really small threshold values will take very long to process and give very poor results, like thousands of very tiny watersheds of little to no value in analysis.

Once you have found you upper most limit for you threshold, determine a few lesser values and build them into a list. Maybe take the upper limit, half it, then half it again.

Make that into a list into a variable:

$ list=$(echo 2000000 1000000 500000)

Then loop the list to iterate through your threshold values when running r.watershed in GRASS. All we change from the last script is the threshold value to use the item from the list. Note we add the threshold value as part of the name too.

Use the standard form of a BASH loop:

for i in $list
    echo $i 

Loop with GRASS watershed commands:

list=$(echo 2000000 1000000 500000)

for i in $list
    # Run watershed operation on fill sink raster
    r.watershed elevation=$fillDEM threshold=$threshold accumulation=$accumulation drainage=$drainage stream=$stream basin=$basin --overwrite

    # Convert Basin (watershed) to vector format
    r.to.vect input=$basin output=$basinVect type=area column=bnum --overwrite

    # Export catchment to vector format
    v.out.ogr input=$basinVect output=$basinVectOut type=area format=ESRI_Shapefile --overwrite

Full script example here

Your outputs this time will be a collection of shapefile watersheds constructed by your threshold value.

Clipping your Raster by Individual Watersheds using GDAL, OGR and SQL

The watersheds layers are full of nice watersheds, however you cannot simply clip the elevation raster by the full vector file. For the clip, you need to select each individual watershed first. The following will demonstrate how to list each of those watersheds and use them as a clip file via a loop.

To start, set up a test directory for your processed raster outputs:

$ mkdir grass_test/raster_watersheds

The basic command for clipping a raster with a vector in GDAL is:

$ gdalwarp -of GTiff -cutline input_shapefile_for cut.shp input_raster_to_be_cut.tif output_raster.tif

You could create a shapefile for each watershed and clip your raster by that watershed individually using the command above, but the process would be fraught with errors and take a very long time.

It is much easier to automate the process using a combination of ogrinfo, sql and gdal.

We’ll do this with an example using JM_basinVectOut_2000000.shp you just created.


In the example, we’ll:

1. Read and list individual watersheds from the vector layer, by listing all the ‘cat’ values leveraging the ogrinfo -sql switch.
2. Using the list, we’ll clip the original raster by each individual watershed

First, list the ‘cat’ values in the vector watershed layer.

The basic ogrinfo command is:

$ ogrinfo /grass_test/JM_basinVectOut_2000000.shp

For this operation; however, we need to list the information of each row in the attributes:

$ ogrinfo -geom=NO -q -sql "SELECT cat FROM JM_basinVectOut_2000000" /grass_test/JM_basinVectOut_2000000.shp

In the above ogrinfo command lists all of the rows for the ‘cat’ column, excluding any geographic information for the file. To do this we use the ‘-sql’ switch available in ogrinfo.

Note when using the ‘-sql’ switch you need a layer name, ‘JM_basinVectOut_2000000’ from the shapefile. The layer name is just the name of the shapefile with the .shp removed.

Use grep to select the parts from the output you want and sed to clean up the outputs for a tidy list.

$ ogrinfo -geom=NO -q -sql "SELECT cat FROM JM_basinVectOut_2000000" /home/ireese/grass_test/JM_basinVectOut_2000000.shp | grep 'cat (Integer)' | sed s/'cat (Integer) =//g'

Make it a list variable:

$ watershedList=$(ogrinfo -geom=NO -q -sql "SELECT cat FROM JM_basinVectOut_2000000" /home/ireese/grass_test/JM_basinVectOut_2000000.shp | grep 'cat (Integer)' | sed s/'cat (Integer) =//')

Before we loop the process in gdal lets look quickly at why this list was created.

In gdalwarp we are going to again leverage the ‘-sql’ switch in order to select each individual watershed row and clip only by that row. We’ll be writing the a basic ‘sql’ similar to the above, but instead filtering by each watershed:

SELECT cat FROM watershed WHERE cat=$i

So, putting all the bits together, using clip and SQL, in gdalwarp we get:

$ gdalwarp -of GTiff -dstnodata -9999 -cutline $inputVector -csql "SELECT cat FROM $inputVectorLayerName where cat='$i'" -crop_to_cutline input_raster_to_be_cut.tif output_raster.tif

Note, I am using ‘-crop_to_cutline’ and ‘-dstnodata’. Basically I’m saying, crop the output raster the bounding box of the input vector and set a ‘nodata’ value of -9999 to the pixels outside the clip area.

Now we need put this all together in a script to loop the process.


# set base path

# Set raster as variable
rasterName=$( basename $raster | sed 's/.tif//g' )

#prep your input vectors
inputVectorLayerName=$(basename $inputVector | sed 's/.shp//')

#create your watersheds list
watershedList=$(ogrinfo -geom=NO -q -sql "SELECT cat FROM $inputVectorLayerName" $inputVector | grep 'cat (Integer)' | sed s/'cat (Integer) = //')

for i in $watershedList
    gdalwarp -of GTiff -dstnodata -9999 -cutline $inputVector -csql "SELECT cat FROM $inputVectorLayerName where cat='$i'" -crop_to_cutline $raster $outDirRast/${rasterName}_${i}.tif

Copy the above script, check your file paths, and save it in a text editor as:



To run the script:

$ bash grass_test/rasterClipByWatershed.sh

Final Note
In the above process, it may be necessary to run some clean up operations to remove invalid geometries in the watersheds vector layer. I did not include this, but it is helpful to resave the watersheds shapefile first by running a buffer operation with a buffer value of ‘0’. This will remove the invalid geometries. You can run this in GRASS or using ogr2ogr with PostGIS. Here is an ogr2ogr with PostGIS example:

$ ogr2ogr -f "ESRI Shapefile" output_vector.shp input_vector.shp -dialect sqlite -sql "select id, ST_buffer(Geometry,0) as geom from input_vector" -overwrite

All the scripts from the past two posts can be found here.

Let me know of this helpful. I could potentially move on to doing some more basic hydro analysis with GRASS and BASH, but these things take time and energy. It would be good to know if these sorts of blogs have value to the wider geospatial world.



But this wasn’t always the case.

GRASS GIS was, for a long time, something I dismissed as ‘too complex’ for my everyday geospatial operations. I formulated any number of excuses to work around the software and could not be convinced it had practical use in my daily work. It was ‘too hard to set-up’, ‘never worked well with QGIS’, and ‘made my scripting processes a nightmare’.

I am here to officially say I’ve been very wrong for a very long time. GRASS GIS is pretty amazing and is a wonderfully easy tool to script.

Recently I scripted a model for extracting river centrelines from high resolution elevation data. During the process I thought how useful it would have been, way back, to have a simple example for setting up an environment and scripting a BASH process using GRASS. So, I built one for myself.

The following example provides the steps for a simple catchment extraction on a piece of LINZ 8m elevation data. The example is only for demonstration of running a basic BASH/GRASS set up with a hydrology command and NOT a demonstration for how to do hydrology using GRASS. The 8m elevation data is not the best data for hydrological extraction; however, these data do provide a nicely sized dataset to use as practice and will give results. Also please note, catchment extraction on a square raster is will not give accurate results at the edges of the raster.

In this example we will:

1. Download a small piece of elevation data from the LINZ Data Service
2. Build a GRASS environment to process these data
3. Build a BASH script to process the catchments
4. Import the elevation into the GRASS environment
5. Perform some basic GRASS operations (fill and watershed)
6. Export raster format for viewing
7. Export the vector catchments to shapefile

The following assumes you are working in a Linux environment and have a basic knowledge of BASH scripting. I tested this process using Ubuntu 18.04 and GRASS 7.4. I have 16GB RAM on my machine.

If you do not all ready have it, you can install GRASS as follows:

$ sudo apt-get install grass-core

If interested, but not necessary, you can install functionality for building GRASS plug-ins:

$ sudo apt-get install grass-dev

Create yourself a directory to work in:

$ mkdir grass_test

Download a piece of elevation data from the LINZ Data Service, place in your directory:


From the Tiles Table, I downloaded the JM tile in EPSG:2193.

Within your directory you will need to build a ‘PERMANENT’ folder in order for GRASS to do its magic. This folder will be set to operate in the projection of your data. Be sure your data is in the same projection as the environment you built. The data is downloaded in EPSG:2193, New Zealand Transverse Mercator (NZTM), so we set the GRASS environment to work in this projection:

$ grass -c epsg:2193 -e grass_test/GRASS_ENV

‘-c’ will create your directory using epsg:2193 and ‘-e’ will exit once this operation is complete.

Running this command will build a PERMANENT folder in:


providing all the bits GRASS needs to operate. Your folder structure will look like:


Have a look inside the folder and see what GRASS built for itself.

The key from now on is to always run your BASH script through this environment. Running your script in this environment will allow you to perform BASH and GRASS commands at the same time.

First, let’s look at the the basic command to run a BASH script using a GRASS environment.

$ grass grass_test/GRASS_ENV/PERMANENT --exec sh grass_test/catchment.sh

The above says, launch GRASS using this environment


and execute, –exec, this script


Let’s build the BASH script. In this we will:

1. Import the elevation data. GRASS likes to work in its own data formats
2. Set the region for where the operation will be performed. GRASS needs to know where the operation is going to be performed.
3. Perform the operations, r.fill.dir and r.watershed
4. Export a raster to .tif format
5. Export the vector outputs to .shp format


# set base path

# Set raster as variable

# Set a base name for the data. This is used to demonstrate that normal
# BASH commands can be used in this process, along side GRASS
rasterName=$( basename $raster | sed 's/.tif//g' )

# Import raster data
r.in.gdal input=$raster output=$rasterName --overwrite

# Set region. IMPORTANT so GRASS knows where the data is located.
# This region is set for the duration of the following commands
g.region rast=$rasterName

# Fill sinks
r.fill.dir input=$rasterName output=$fillDEM direction=$directionDEM areas=$areasDEM --overwrite

# Export a raster for viewing
r.out.gdal input=$areasDEM output=$areaOut

# Run watershed operation on fill sink raster
r.watershed elevation=$fillDEM threshold=$threshold accumulation=$accumulation drainage=$drainage stream=$stream basin=$basin --overwrite

# Convert Basin (watershed) to vector format
r.to.vect input=$basin output=$basinVect type=area column=bnum --overwrite

# Export catchment to vector format
v.out.ogr input=$basinVect output=$basinVectOut type=area format=ESRI_Shapefile --overwrite

There you have it. A simple BASH script, set up for running some GRASS commands, running on your LINUX machine.

Putting it all together

Set up you environment, copy the above script, paste it into a text editor, save it as:


Now run

$ grass grass_test/GRASS_ENV/PERMANENT —exec sh grass_test/catchment.sh

Your output catchments shapefile should be similar to the image below



Click here for the next step in this process; how to clip your raster using the watershed layer using GDAL and loops.

Let me know if this was helpful or if you would like to see any changes.

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
The visualization is 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.”