QGIS TO STYLIZED COGs

Image | | xycarto

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:

ScaleResolutionZoom
250000707
100000288
50000149
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 \
	-co COMPRESS=JPEG

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: [
	{
  	url:url,
	},
  ],
  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.

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

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

Image | | xycarto

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 = [
    "BIGTIFF=YES",
    "BLOCKSIZE=512",
    "RESAMPLING=BILINEAR",
    "COMPRESS=DEFLATE",
    "NUM_THREADS=ALL_CPUS"
]

gdal.Translate(    
    "clipped-eez-nztm-20200101090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02_fill_cut_warp_cog.tif",
    “warp_anti.tif”,
    format = "COG",
    callback=gdal.TermProgress_nocb,
    creationOptions = creation_options
)

I put the COG on S3 and made the link public

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

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: [
    'interpolate',
    ['linear'],
    cogBand,
    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 () {
overlay.setPosition(undefined);
closer.blur();
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);
    console.log(data[0])
    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>";
    overlay.setPosition(coordinate);
  })

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!

Warping Across Anti-meridian: NZTM to Web mercator

Image | | xycarto

TL;DR

Python example here.

NZTM to Web Mercator

There are a number of times when I am building a web map when I need to warp a GTiff across the anti-meridian. New Zealand sits at the edge of this line, with a portion of it on the other side, The Chatham Islands. Doing a direct gdalwarp from NZTM (2193) to Web Mercator (3857) will “split” the file with the portion across the anti-meridian being placed to the other side of the projection.

Doing a straight warp from NZTM to Web Mercator:

gdalwarp -s_srs EPSG:2193 -t_srs:3857 in_tif.tif webmer_out_tif.tif

produces the following:

A number of people have devised methods to get around this, but I wanted to post one way I found that could be useful, albeit verbose. This is in by no means the best method; however, I have found it useful as a solution when I have a request to build a visualization using Web Mercator.

The short description is this:

The goal is to get the extents in Web Mercator for your GTiff that cross the anti-meridian. These extents can be passed to gdal.Warp as the output extents in your warp command

A Github repo with an example script can be found here. The code is intentionally verbose to help keep track of what is happening along the way.

The example uses Python and a number of available libraries, including:

The real star of the method is fiona. Fiona has a method in the command that allows reprojection across the anti-meridian: antimeridian_cutting=True. When using this, fiona will output a multi-polygon in two pieces similar to the image above. Using the individual extents from each polygon in the multi-polygon, the extents can be leveraged to reconstruct a polygon across the anti-meridian. The reconstructed polygon extents can then be used as the output extents of you gdalwarp operation.

A little more in depth:

  • Find extent of input file
  • Reproject this extent to WGS 84 using fiona with antimeridian_cutting=True. This creates a multi-polygon
  • Find extent of each multi-poly
  • Shift left polygon over anti-meridian, next to right polygon
  • Build extent of polygons combined
  • Reproject WGS polygon to Web Mercator
  • Shift maxx of Web Mercator polygon across the anti-meridian
  • Find extent of new Webmercator polygon
  • Use new extent as output bounds in gdal.Warp process

Again, here is link to the repo with an example.

The example in the repo uses data from National Institute of Water and Atmospheric Research (NIWA) in New Zealand. You can get a download of the example data via Koordinates here.

Normally, I would try to build the code in the blog post, but the example too large. Also, I usually try to break down all the import highlights in the code; however, there is just too much to cover in the method. Hopefully, this will get you a little further along.

Static Vector Tiles II: Openlayers with Custom Projection

Image | | xycarto

Building off the work from my previous post on vector tiles, I wanted to develop a second process for use in Openlayers. There are many tutorials out there demonstrating the use of vector tiles with predefined styles and tutorials explaining how to build a fully vector tile map online. My aim instead is to demonstrate the processes to build:

  1. A vector tile set up overlaying a raster tile cache
  2. How to integrate a custom style.json
  3. How to do all of this in a custom projection using Openlayers

As always, there is more than one method to do any of this. What is shown below can definitely be improved upon. Please take this and make it better.

Github Respository is here.

Website example is here.

This example uses Openlayers6. The map is served in the NZTM projection (EPSG:2193). No effort has been made in ordering the labels. The labels display as the Openlayers decluttering orders them. Included in the JS code is a basic pop-up window demonstrating how to get information from the vector tile.

The project is built as a static set up. The vector tile cache is built directly into the website. THIS IS NOT OPTIMAL, but does demonstrate the principle. Ideally, you would have a location like AWS S3, to serve your tile cache from.

In order to use a custom projection, you will need to build an XYZ tile cache. MBTiles do not handle projections other than Web Mercator (EPSG:3857).

Basic steps

  1. Download or reproject the shapefile in NZTM
  2. Upload shapefile to PostgreSQL database with PostGIS extensions
  3. Tile PostgreSQL table into NZTM (EPSG:2193) XYZ tile cache using TRex
  4. Construct Openlayers6 JS for tile consuption

Sample Data

https://data.linz.govt.nz/layer/50280-nz-geographic-names-topo-150k/

Note:

The Geographic Names layer is clipped and filtered for this example. I clipped only to the Wellington Region and filtered the data only to use:

desc_code = BAY, METR, LOC, POP, TOWN, SBRB

Upload to PostgreSQL

shp2pgsql -s 2193 /data/wellyRegion_townBay.shp public.wellyRegion_townBay_nztm | psql -h localhost -d <yourDatabaseHere> -U <youUserNameHere>

TRex Tiling

TRex will create an XYZ tile cache in the projection of your choosing. You will need to know the resolutions and bounding box of your projection in order to make this work. I was fortunate to have this information at hand thanks to a great tutorial from LINZ.

TRex uses a config file for tiling. The config used in this example is here

The command used to run TREX:

t_rex generate --progress true --maxzoom=14 --minzoom=0 --extent=174.627603,-41.613839,176.259896,-40.737190  --config /configpsql_points.toml

TRex will generate gzip pfb’s. If you prefer to unzip them:

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

Openlayers JS

The Openlayers for this is version 6.  <script> tags needed are:

<link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.5.0/css/ol.css" type="text/css">

<script src="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.5.0/build/ol.js"></script>

//cdnjs.cloudflare.com/ajax/libs/proj4js/2.3.15/proj4.js

https://unpkg.com/ol-mapbox-style@6.3.2/dist/olms.js

For the full JS example

NZTM Construct in Openlayers

Building the projection for Openlayers

// set NZTM projection extent so OL can determine zoom level 0 extents.
// Define NZTM projection
proj4.defs("EPSG:2193","+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");

// Register projection with OpenLayers
ol.proj.proj4.register(proj4);

// Create new OpenLayers projection
var proj2193 = new ol.proj.Projection({
	code: 'EPSG:2193',
	units: 'm',
	extent: [827933.23, 3729820.29, 3195373.59, 7039943.58]
});

// NZTM tile matrix origin, resolution and matrixId definitions.
var origin = [-1000000, 10000000];
var resolutions = [ 8960, 4480, 2240, 1120, 560, 280, 140, 70, 28, 14, 7, 2.8, 1.4, 0.7, 0.28, 0.14, 0.07 ];
var matrixIds = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16];

Applying to Raster and Vector Tiles

Raster Tiles

Another great tutorial from LINZ regarding the set up of an XYZ for raster tiles.

// Tile Services Map
var urlTemplate =
  "https://tiles.maps.linz.io/nz_colour_basemap/NZTM/{z}/{x}/{y}.png";

// Set raster layer
var layer = new ol.layer.Tile({
  source: new ol.source.XYZ({
    url: urlTemplate,
    projection: proj2193,
    attributions: ['<a href="http://data.linz.govt.nz">Data from LINZ. CC BY 4.0</a>'],
    tileGrid: new ol.tilegrid.TileGrid({
      origin: origin,
      resolutions: resolutions,
      matrixIds: matrixIds,
      extent: [827933.23, 3729820.29, 3195373.59, 7039943.58]
    })
  })
});
Vector Tiles

Set up the vector layer to use custom projection:

// Set vector layer
var placeSource = new ol.source.VectorTile({
  cacheSize: 0,
  overlaps: true,
  tilePixelRatio: 1, // oversampling when > 1
  tileGrid: new ol.tilegrid.TileGrid({ 
    origin: [-1000000, 10000000],
    maxZoom: 16,
    tileSize: 4096,
    extent: [827933.23, 3729820.29, 3195373.59, 7039943.58],
    resolutions: resolutions,
  }),
  extent: [827933.23, 3729820.29, 3195373.59, 7039943.58],
  format: new ol.format.MVT(),
  projection: ol.proj.get('EPSG:2193'),
  url: 'https://xycarto.github.io/static.vector.tiles.openlayers.nztm/tiles/wellyRegion_townBay_nztm/{z}/{x}/{y}.pbf'
});

var vectorMap = new ol.layer.VectorTile({
  declutter: true,
  source: placeSource,
  renderMode: 'vector',
  zIndex: 10
  })

Styling

For the style file example. 

  1. The method in this example is loading the vector tile and overaying it on a raster tile cache. In order to accomplish this, a vector tile cache must be loaded first to the map, THEN the rules from the style JOSN are applied using:
fetch('./styleText.json').then(function(response) {
  response.json().then(function(glStyle) {
    olms.applyStyle(vectorMap, glStyle, 'wellyRegion_townBay_wgs');
  });
});
  1. The above uses olms.applyStyle. To access this function you will need to add the scipt tag to your HTML:
https://unpkg.com/ol-mapbox-style@6.3.2/dist/olms.js

Notes

  1. Not fully complete. Working example only
  2. Runs slow in Safari
  3. This is in no way the only way to do this. My hope is someone takes this and make it better.
  4. Still needs label work for use on mobile devices.

RGB Elevation Creation for 3D Online Mapping (Terrain RGB)

Image | | xycarto
Mt. Ruapehu, Ngauruhoe, and Tongariro

The following is about how to build an elevation file to work with procedural-gl.js. I just really wanted to build my own elevation dataset and thought it would be helpful to share how I did it.

Procedural-gl.js has changed the game in attainable 3D online mapping. When I came across this project, I was very excited and wanted to get to mapping ASAP.  There is so much potential in just a few lines of JS code; just add a base map, an elevation file and boom, you have a smooth functioning 3D map, mobile ready. I cannot thank Felix Palmer enough for putting this out there. You can check out the source code and some really great maps he developed here.

I’m not one to leave well enough alone. I wanted to develop a method to use New Zealand specific elevations with my goal being to eventually incorporate New Zealand’s LiDAR elevations. However, before I get to building elevation models for LiDAR, I wanted to test building an elevation model with the LINZ 8m DEM. This dataset covers the entire nation and was light weight enough to test on the ol’ home Linux box.

Method

As mentioned above, to run procedural-gl.js in its most basic form, you’ll need a tile cache for your base map and a tile cache for your elevation file.  The base map I have covered, but I was unfamiliar for what was needed for the elevation tile cache. Fortunately, I came across SyncPoints tutorial on how to do this.

I am not going to rewrite their entire process.  They did a really good job explaining the what, why, and how.  Instead, I will layout what I did for the NZ specific data. Much of this process can be further refined and this is really just the proof of concept.

Two quick notes before I get started:

  1. For now, this process is only done in Web Mercator.  I always want to push the limits and do my online work in NZTM, but for now, I am going to work with what I have.
  2. I am not going to talk about building a raster tile cache for the base map.  I was fortunate to have a web mercator tile cache for NZ available.

We want to create an elevation file, in PNG rgb format, rendered into an XYZ tile cache directory, with 512x512px tiles.

The basic steps are this:

  1. Remove all errant noData pixels from a geotiff.  
  2. Convert your elevation geotiff to rgb (single band to three band), in Byte format.
  3. Render your rgb geotiff into xyz raster tiles. Use a directory format.
  4. Consume through procedual-gl.js

List of Tools

Data

The data I used for this was the full LINZ 8m DEM elevation. Make sure your data is in Web Mercator projection (EPSG:3857).  See some earlier blogs about how to retroject across all the tiles using BASH and gdalwarp.

LINZ 8m DEM

Process

1. With the VRT, we can gather all the elevation tiles under one file name.  

gdalbuildvrt elevation.vrt *.tif

2. Remove nodata

Using the VRT, we can preform the noData conversion across all the elevation tiles as if it were a mosaic.  Better yet, the output from this command will produce a mosaic out the other side. GDAL is magic.

gdalwarp  -t_srs EPSG:3857 -dstnodata None  -co TILED=YES  -co COMPRESS=DEFLATE  -co BIGTIFF=YES  elevation.vrt elevation_noData_mosaic.tif

3. RGB-ify the elevation mosaic.

rio rgbify -b -10000 -i 0.1 /elevation_noData_mosaic.tif elevation_noData_mosaic_rgb.tif

4. Create your XYZ tile cache in EPSG:3857

GDAL2Tiles is a quick way to render xyz tile caches if you have single layer and are working in Web Mercator.  Be sure to use 3.1 or greater to get the functionality of xyz caching and tile size control.  You need these both.  In the end, I used the GDAL Docker developed by perrygeo, with modification, --with-python, and got access to GDAL 3.2.  It was a lifesaver, however, I did need to modify the Dockerfile to add GDALs Python bindings and rebuild the Docker.

gdal2tiles.py --s_srs=EPSG:3857 --zoom=0-16 --xyz --tilesize=512 --processes=7 elevation_noData_mosaic_rgb.tif /data/tile-cache

Localhost Testing

Once I built the tile cache, I built a quick site on my localhost to test. I embedded the tile cache into the site (place the tile cache in the same directory as your web files), then I ran it under a local server.  The advantage of embedding your tile cache for testing is that it allows the local server to serve out the raster tiles as well. There is no need to set up a tile server. The local server does this for you.

I switched it up this time used the npm http-server to avoid some cors issues I was encountering. 

http-server . —cors 8000

Code and Website

You can view the code for the HTML and JS here:

https://github.com/xycarto/3D_Aotearoa/

The running website can be found here:

https://xycarto.github.io/3D_Aotearoa/

It works great on mobile devices, but will happily display on your desktop.