[GIS] raster tiles from vector data with GDAL – how to avoid aliasing artefacts

Question

Prequesites

I have vector data in a Shapefile N47E007_CONTOURS.shp produced with

wget -O N47E007.zip http://e4ftl01.cr.usgs.gov/SRTM/SRTMGL1.003/2000.02.11/N47E007.SRTMGL1.hgt.zip
unzip N47E007.zip

# Convert from WGS 84 to Web Mercator and from HGT to GeoTIFF
gdalwarp -s_srs EPSG:4326 N47E007.hgt \
         -t_srs EPSG:3785 N47E007_WebMercator.tif \
         -co "COMPRESS=LZW" -of GTiff -r bilinear

# Compute contour vector data
gdal_contour -a elev -i 20 N47E007_WebMercator.tif N47E007_CONTOURS.shp

(Contour lines based on one patch (N47E007) of NASA SRTM elevation data.)

Process

From this, I generate TMS/WMTS tiles to be used by leaflet, by first rasterizing the vector data from the Shapefile with gdal_rasterize and then tiling it with gdal2tiles.py. I repeat both steps for every zoomlevel, so I can rasterize to a larger raster for higher zoomlevels, so that lines will have the same on-screen thickness on each one:

#!/bin/sh
for z in $(seq 1 15)
do
    gdal_rasterize \
        --config GDAL_CACHEMAX 1024 \
        -ts $(calc "2851 << (${z} - 12)") \
            $(calc "4220 << (${z} - 12)") \
        -ot byte -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 \
        -ot Byte -burn 0 -a_nodata 255 \
        -l N47E007_CONTOURS N47E007_CONTOURS.shp N47E007_CONTOURS_rasterized_z${z}.tif
    gdal2tiles.py \
        -r near -z${z} -a 255,255,255 \
        -e N47E007_CONTOURS_rasterized_z${z}.tif \
        N47E007_contours_tiled
done

Adapting resolution to zoomelevel

(See also my answer to Problem with contour lines thickness in bigger zoom levels.)

calc "2851 << (${z} - 12)" and calc "4220 << (${z} - 12)" compute the size of the image to rasterize to (N47E007_CONTOURS_rasterized_z${z}.tif) depending on the zoomelevel $z.

<< is the binary left-shift operator, so this translates to 2851⋅2 ^z-12 and 4220⋅2 ^z-12. This means that

• for zoomlevel 12, the intermediate raster image
  N47E007_CONTOURS_rasterized_z12.tif will be 2851×4220 pixels,
• for zoom level 13 (N47E007_CONTOURS_rasterized_z13.tif) it will be
  twice that in each direction (5702×8440) and
• for zoom level 11 (N47E007_CONTOURS_rasterized_z11.tif) half that
  of zoomlevel 12, i.e. 1425×2110.

Imprecise parametrization

The numbers 2851 and 4220 (the size for zoomelevel 12) were guessed or measured on screen or otherwise derived in an imprecise manner.

Problem

The resulting tiles show some artefacts that aren't visible in the intermediate (untiled) rasterized images:

• double-width segments of lines,
  e.g. 13/4270/5331.png vs. corresponding excerpt of N47E007_CONTOURS_rasterized_z13.tif

  vs.

• gaps in lines,
  e.g. at tile edges between

  14/8536/10665.png 14/8537/10665.png
14/8536/10664.png 14/8537/10664.png

  
  

  corresponding excerpt from N47E007_CONTOURS_rasterized_z14.tif for comparison:

If I skip -r near (so rasterization uses the default resampling mode 'average') I instead get some washed-out pixels which make the lines look unsharp.

My theory

The coordinate reference systems already match (due to the gdalwarp step before calculating the contours), so no general transformation/re-projection should be required. But because the size to rasterize to only approximates the on-screen size of the geographic area covered by the N47E007 data patch, some rescaling must still happen in gdal2tiles.py, to produce correctly placed and sized seamless tiles.

With -r near this causes some kind of macroscopic Moiré effect, because the pixel grids of the tiles and of the intermediate raster only almost match. (Additionally, gdal2tiles.py might not look accross tile boarders when looking for the nearest pixel in the source image, causing the gaps.) With average (the default for -r), the same causes effect the washed-out look because values from more than one source pixel get 'mixed' to get the value for a given tile pixel.

Thus, if I knew the correct (zoom-level-dependant) on-screen-size of the geographic area of the data patch and would use that for rasterization, the artifacts should go away.

Questions

• Is my theory correct?
• If so, how can I determine the precise -ts parameters for gdal_rasterize (or, alternatively, -tr parameters) so that no rescaling (or a reschaling without effect) will be done by gdal2tiles.py?

Answer 1

Is my theory correct?

It seems so, as I get flawless results with the method described below.

If so, how can I determine the precise -ts parameters for gdal_rasterize (or, alternatively, -tr parameters) so that no rescaling (or a reschaling without effect) will be done by gdal2tiles.py?

While gdal2tiles.py is meant to be used as a command, we can load it as a python module and use the very same functions it uses itself for the size computations. Thus, I've made me this little python script:

from __future__ import print_function
import sys
from gdal2tiles import GlobalMercator

def main():
    lat, lon = [ float(arg) for arg in sys.argv[1:3] ]
    zoom = int(sys.argv[3])
    px, py = lat_lon_to_pixels(lat, lon, zoom)
    print(px, py)

def lat_lon_to_pixels(lat, lon, zoom):
    gm = GlobalMercator()
    mx, my = gm.LatLonToMeters(lat, lon)
    return tuple( int(round(p)) for p in gm.MetersToPixels(mx, my, zoom) )

if __name__ == '__main__':
    main()

I then use the WGS 84 coordinates from gdalinfo N47E007.hgt (exctracting them with the regular expressions /^\s*Upper Right\s*\(\s*(?<east>\d+\.?\d*),\s*(?<north>\d+\.?\d*)\) and /^\s*Lower Left\s*\(\s*(?<west>\d+\.?\d*),\s*(?<south>\d+\.?\d*)\)/) to call lat_lon_to_pixels from the python script above for the north-east corner and the south-west corner and compute the differences, which I then use as argument for the -ts option of gdal_rasterize. (This only works because this HGT file has WGS 84 as coordinate reference system and WGS 84 and WebMercator are axis-aligned to each other.)