Tiling an RGB Raster

Continuing on from the previous post, this entry will describe a workflow for saving the satellite image mosaic created using Google Earth Engine to an RGB raster file. The image will then be processed to scale its values, and tiled for use as a web map layer.

Exporting the Image

After running the GEE script from the last post, prepare a layer for export by selecting bands from the rgb layer and sorting them in RGB format (i.e. bands 4, 3, 2). The resulting RGB layer can then be exported to a GeoTiff (requires Google Drive with sufficient storage available) by additionally running the code below and then clicking Run > Set Output in the Tasks panel. Note that the export may take quite a long time to run (in my case 39 mins), but once started the GEE editor can be closed. For starting in a fresh editor session, a full script is available here, which skips the visualization and other unneeded elements of the original script, and then runs the export.

// Reorder bands to correct 'True Color' RGB ordering.
var sentinel2mosaic = rgb.select('B4_median').rename('b4')
  .addBands(rgb.select('B3_median').rename('b3'))
  .addBands(rgb.select('B2_median').rename('b2'));

// Export the image, specifying the CRS, scale, and region.
Export.image.toDrive({
  image: sentinel2mosaic,
  description: 'sentinel2mosaic',
  crs: 'EPSG:4326',
  //crsTransform: affine,
  region: geom,
  scale: 10
});

Processing

One the export is complete, the file can be downloaded from Google Drive and processed. The intent is to tile the image using gdal2tiles, but if we check the data type of the file with gdalinfo as gdalinfo sentinel2mosaic.tif -mm, we can see the values are stored as a Float64 data type, whereas the gdal2tiles documentation states the following:

Inputs with non-Byte data type (i.e. Int16, UInt16,…) will be clamped to the Byte data type, causing wrong results. To avoid this it is necessary to rescale input to the Byte data type using gdal_translate utility.

This means before we can tile the image, we first need to convert the data type from Float64 to Byte, and scale the 64 bit values to fit within this greatly reduced range (0-255). As recommended above, we can use gdal_translate for this as:

gdal_translate -ot Byte -scale sentinel2mosaic.tif sentinel2mosaic-8bit.tif

Although if you then open the resulting file in QGIS, you will likely see that the image appears very dark:

To fix this, right click on the layer and select Properties > Symbology, and under Min / Max Value Settings, set Cumulative count cut, and Contrast enhancement as Stretch to MinMax, then hit Apply:

What happened here was the extreme top and bottom 2% of values were excluded from the color ramp, which were otherwise setting a too large range that allocates unexpected colors to the majority of values. This fix only applies to the layer in QGIS, however, so to fix the file itself, we too need to exclude the top and bottom 2% of values when scaling. One way to get these values can be with a simple Python script as below:

python3
import gdal
import numpy as np

ds = gdal.Open('sentinel2mosaic.tif')
for n in range(1, 4):
    band = ds.GetRasterBand(n)
    arr = band.ReadAsArray()
    sel = np.isnan(arr)
    p2 = np.percentile(arr[~sel], 2)
    p98 = np.percentile(arr[~sel], 98)
    print('{}: {} {}'.format(n, p2, p98))

Which returns the values:

1: 147.5 1341.0
2: 292.5 1184.0
3: 268.0 955.0

We can then re-run the scaling and override the default min-max values with the 2nd and 98th percentiles calculated above as:

gdal_translate -ot Byte -scale_1 147.5 1341.0 0 255 -scale_2 292.5 1184.0 0 255 -scale_3 268.0 955.0 0 255 sentinel2mosaic.tif sentinel2mosaic-8bit.tif

Loading this to QGIS should look fairly comparable to the last layer – but now without needing to tweak the symbology.

Tiling

The image can now be tiled using gdal2tiles with something as below, exporting tiles for zoom levels 2 to 12 to a directory called tiles.

gdal2tiles.py -z 2-12 sentinel2mosaic-8bit.tif tiles

A leaflet.html file is also exported, containing a web map, which with a few tweaks (e.g. updating the tiles url) can be uploaded online.

SRTM Elevation Map

Elevation map of Ontario, lighter shades are higher elevations.

Digital Elevation Models (DEMs) are maps of the bare topographic surface of the Earth, excluding features like trees and buildings etc, and usually collected by satellites. DEMs can be used to map topographic elevation, relief, slope, and aspect, and are useful sources for a range of purposes such as flood risk modelling, contour mapping, correcting aerial and satellite imagery, and land cover classification.

A well known source of this data is the USGS/NASA Shuttle Radar Topography Mission (SRTM), which provides tiled DEM data for the entire globe within latitudes 60° N and 56° S, at a resolution of 1 arc-second (~30m).

In this post, I describe some methods to source and process SRTM imagery. Before you can download the data though, you’ll first need to register for a free Earthdata account.

PEC Sample DEM

As an example, we can first create an elevation map for a smaller region, in this case Prince Edward County, Ontario (the peninsular in north eastern lake Ontario). For downloading a handful of tiles, this useful tool created by Derek Watkins can be used to select tiles from an interactive map: http://dwtkns.com/srtm30m/.

Download and unzip the four tiles that cover Prince Edward County, N44W078 (pictured below) and the three others east and south of it: N44W077; N43W077; and N43W078, entering your Earthdata username and password when prompted. Also, download Canadian census divisions data from Statistics Canada (Cartographic Boundary File), which we can use to clip the SRTM data.

SRTM download by tile

With the four tiles downloaded and unzipped to .hgt files, we can merge them using GDAL, building a VRT file, a virtual dataset useful for intermediate processing steps. The below example uses Docker and the osgeo/gdal image, but can be modified for other configurations.

docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdalbuildvrt /mnt/data/pec-srtm.vrt /mnt/data/N44W077.hgt /mnt/data/N44W078.hgt /mnt/data/N43W077.hgt /mnt/data/N43W078.hgt

We can then clip the VRT file to the census divisions layer with gdalwarp as below, selecting the ‘Prince Edward’ feature as a cutline. This creates a clipped DEM GeoTiff like that displayed at the head of this post.

docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdalwarp -cutline /mnt/data/lcd_000b16a_e/lcd_000b16a_e.shp -cwhere "CDNAME = 'Prince Edward'" -crop_to_cutline /mnt/data/pec-srtm.vrt /mnt/data/pec-dem.tif

Some alternative calculations and visualizations can be created with the gdaldem tool as below, for example hillshade (shaded relief to represent 3d terrain), slope, and aspect (the direction of the slope).

docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdaldem hillshade -s 111120 /mnt/data/pec-dem.tif /mnt/data/pec-hillshade.tif

docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdaldem slope -s 111120 /mnt/data/pec-dem.tif /mnt/data/pec-slope.tif

docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdaldem aspect /mnt/data/pec-dem.tif /mnt/data/pec-aspect.tif

Entire Ontario DEM

Needless to say, it would be tedious and time consuming to manually download individual tiles to cover the entirety of Ontario (or beyond) in this way. So to scale up to a larger area we can write a script to automate the process using e.g. Python, like shown in the notebook below (available for download here).

This script uses a layer of the SRTM grid extent and a Provinces/territories boundaries layer, again from Statistics Canada, overlaying the two to determine which tiles to download. Note that the entire selection consists of 170 files, totalling ~750mb, which make take some time to download.

With the tiles downloaded we can again use GDAL to merge and crop the tiles in a similar way to the PEC sample e.g.

docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest sh -c 'gdalbuildvrt /mnt/data/on-srtm.vrt /mnt/data/*.hgt'

docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdalwarp -cutline /mnt/data/lpr_000b16a_e/lpr_000b16a_e.shp -cwhere "PRNAME = 'Ontario'" -crop_to_cutline -co COMPRESS=LZW -co BIGTIFF=YES -co TILED=YES /mnt/data/on-srtm.vrt /mnt/data/on-dem.tif

Addendum: Web Map

I you want to use the results for a web map, web tiles can be created using gdal2tiles.py, after first converting the data type via gdal_translate. The below commands generate web map tiles at zoom levels 4-10, along with the HTML content for rendering a web map, such as leaflet.

docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdal_translate -of VRT -ot Byte -scale /mnt/data/on-dem.tif /mnt/data/temp.vrt

docker run --rm -it -v $PWD:/mnt/data osgeo/gdal:ubuntu-small-latest gdal2tiles.py -e -x -t "Ontario DEM" --zoom 4-10 --processes 2 /mnt/data/temp.vrt /mnt/data/web-tiles

All roads vector map: part 1

Welcome to my blog, all about mapping Ontario, Canada (and beyond) with free and open source GIS tools and data. I am starting out with a three part tutorial, which will demonstrate creating the ON roads web map from the maps page, using a range of GIS tools, including GDAL/OGR, tippecanoe, tileserverGL, maputnik, and mapbox-gl.js. As installing these libraries may not be straight forward on all operating systems, and out of scope for this blog, I suggest using Docker and running dedicated containers for each task, where applicable. QGIS is also recommended, as a general purpose method for reviewing spatial data.

We will first need some GIS data to work with. As a case study, I am using the Ontario Road Network (ORN), which may be found at the Ontario Geohub. To follow along, first download and unzip the data, selecting the ESRI Shapefile format. Optionally preview the data with QGIS, if installed.

To run the vector tile generation tool tippecanoe, we next need to transform our shapefile to the desired GeoJSON format. To do this we can use the ogr2ogr spatial data conversion tool, available as a stand alone application or with a docker image like osgeo/gdal:alpine-small-latest, as with the example below.

This command reads the road network layer from the current directory ($PWD, mounted as /mapping-on), and filters the output attributes to only retain the geometry column using SQL (-dialect & -sql), to reduce the file size. The transformed and filtered output is written to a file called road.geojson. Omit the code before ogr2ogr and adjust the file paths if installed stand alone.

docker run --rm -v $PWD:/mapping-on osgeo/gdal:alpine-small-latest ogr2ogr -dialect sqlite -sql "select geometry from Ontario_Road_Network_ORN_Road_Net_Element" /mapping-on/road.geojson /mapping-on/Ontario_Road_Network_ORN_Road_Net_Element/Ontario_Road_Network_ORN_Road_Net_Element.shp

With the roads data now transformed to GeoJSON format, we can use tippecanoe to generate vector tiles for mapping. In brief, vector tiles are a compressed and efficient method of delivering spatial data to a browser for client side styling and rendering, in contrast to the more “traditional” server side raster (image) rendering, which can be resource intensive. Tippecanoe is one of many tools that can generate vector tiles from other data sources.

The simplest way to get started with tippecanoe is again via docker, as below. This command generates an mbtiles database of vector tiles from zoom levels 4 (low) to 12 (high), and a low-zoom detail level of 10 (from a default of 12, lower is less detail). Tippecanoe will simplify, aggregate and filter the data at each zoom level to reach an optimized tile size for mapping, while aiming to retain the shape, density, and distribution of the original data.

docker run --rm -v $PWD:/mapping-on ssastbury/tippecanoe:latest -f -D10 -Z4 -z12 -o /mapping-on/road.mbtiles /mapping-on/road.geojson

With the vector tiles generated, we can now think about serving them for styling. A free and open source way to do this is with TileserverGL, a user friendly server for vector tiles. Running the klokantech/tileserver-gl docker image with default parameters as below, will search for a .mbtiles file in the current directory, and serve the tiles on port 8080 of your machine.

 docker run --rm -it -v $PWD:/data -p 8080:80 klokantech/tileserver-gl 

Point your browser to localhost:8080 and you should see tileserver up and running. Click the Inspect button to generate a basic slippy map with an interactive preview of your tiles. The next job is to apply our own custom styling to the data, which is the topic of my next post.