raster tiles

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.

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.