Aligning Rasters to Webmercator grid

On raster import, rasters need to be projected to webmercator first (srs=epsg:3857). Then they need to be scaled to the closest scale of the webmercator grid, preferring an higher scale over a lower scale.

• Start from the scale of the raster, example:

  input_scale=`gdalinfo HYP_HR-wm.tif | grep '^Pixel Size' | sed 's/.*(\([^,]*\).*/\1/'`

• Divide the resolution of tiles in zoom 0 of webmercator grid (z0=156543.03515625) by the input scale

  factor=`echo "scale=10; $z0/$input_scale" | bc`

• Find the closest power of 2 to the factor (on the lower bound)

  pow2=`echo "pw=l($factor)/l(2); scale=0; pw/1" | bc -l`

• Elevate the factor to the found closer power of 2:

  out_scale=`echo "scale=10; $z0/(2^$pow2)"| bc`

• To be sure the final scale is not (for appriximation reasons) higher than the official grid value, it can be reduced by a 0.01% of its value:

  out_scale=`echo "scale=10; $out_scale-($out_scale*0.0001)" | bc`

Now that the output scale has been computed, we use gdalwarp to downscale the raster:

gdalwarp -tr $out_scale -$out_scale HYP_HR-wm.tif HYP_HR-wm-al.tif

The output raster should then have a pixel scale which is 0.01% lower than the one for zoom level = $pow2. The percentage of reduction could be reduced, or even omitted, after verifying that the scale does not become slightly higher than the target (on this reguard we could maybe allow the mapnik pgraster plugin to use a tolerance when looking for appropriate overviews, to avoid picking the higher resolution level for a small difference).

Important notes

This is the formula proposed above:

out = z0 / (2^floor(log2(z0/in)))

where:

• z0 is the webmercator grid size
• in is the pixel size in meters of the input raster
• out is the output pixel size in meters.

The goal of the reprojection described above is to adjust the base raster and its overviews to well-known zoom levels for rendering, so less work is required scaling tiles.

In experiments implementing this algorithm, the results obtained for request time at different zoom levels for http://www.naturalearthdata.com/downloads/10m-raster-data/10m-shaded-relief/ (large size):

for zoom in $(seq 0 12); do
    for i in $(seq 1 40); do
	url="http://development.localhost.lan:8282/api/v1/map/3e6d7e4df386debf1fd79042e16263dc:1418139427175.6401/$zoom/0/0.png"
	curl -w "%{time_total}\n" -o /dev/null -s $url;
    done | awk '{ sum += $1; n++ } END { if (n > 0) print sum / n; }';
done

Zoom level	Adjusting scale (s)	Vanilla (s)
0	0.029	0.04385
1	0.040775	0.053775
2	0.0324	0.051875
3	0.030675	0.052975
4	0.027625	0.047375
5	0.027	0.04535
6	0.037475	0.031075
7	0.035325	0.043375
8	0.109175	0.0431
9	0.390525	0.136025
10	1.57533	0.4952
11	6.44325	1.95802
12	26.334	8.03935

In conclusion:

• This improves performance in the average case
• There's a resolution loss as this is written
• Higher zooms must be prevented in the UI, because they can overload mapnik
• To prevent the resolution loss, there are two things to do:
  • avoid chaining two calls to gdalwarp. Pixel size in the target ref system can be calculated using gdaltransform
  • use round or ceil instead of floor to adjust the power of 2 to be used in the adjust formula

There's also something worth mentioning: gdalwarp will also adjust the size of the raster if no output size is specified, to a value convenient for further transformations:

-et err_threshold: error threshold for transformation approximation (in pixel units - defaults to 0.125).