created:
Last modified:
Permalink: kreidefossilien.de/1859

Create a nice looking relief-hillshading map and use it in OsmAnd

tags: hillshading, elevation, maps, qgis, webgis, osmand
Hillshading of Saxony (DGM20, GEOSN - License: <a href='https://www.govdata.de/dl-de/by-2-0'>dl-de/by-2-0</a>) in OsmAnd – © OpenStreetMap contributors (License: <a href='https://www.openstreetmap.org/copyright'>ODbL</a>)
Hillshading of Saxony (DGM20, GEOSN - License: dl-de/by-2-0) in OsmAnd – © OpenStreetMap contributors (License: ODbL)
In September 2019 the State Survey of Saxony (Germany) published a few of their datasets under a "Open-Data" license. Datasets like the satellite imagery (DOP) or the digital elevation models (DEM - German: Digitales Geländemodell DGM) are now available for free, without any charge and it seems with a proper license for even a commercial reuse. The DEM-datasets are offered in different resolutions and level of details, ranging from the DGM1, DGM2 up to DGM20. The latter one seems to be detailed enough to enhance a offline navigation app like OsmAnd+ for Smartphones, but not having troubles with file sizes of hundreds MB or even Gigabytes of data. Creating an enhanced hillshading layer with a tolerable file size was the aim of this small project.

The article describes the process of creating a custom hillshade image for the usage in the OsmAnd App based on the now public available DGM20 dataset of Saxony (Germany). This article will focus on the technical details and solving the problems I had, whilst creating it. The result is a much more detailed hillshade-layer than the SRTM, initially served by the Osmand+ hillshade-contour-plugin. The hillshade layer in this example only covers the area of the Free State of Saxony (Germany).

If you just landed here to get and use the result - it's a sqlitedb - for the OsmAnd Application scroll down to the very end of this page and download the *.zip. READ the instructions — at least the one in the readme.txt

This article addresses users who are comfortable in working with a GIS, in this case with QGIS. You should be already familiar with the "concepts" of geographic projections and certain requirements for the used file formats and applications. You should have a slight idea what some of the used QGIS/GDAL-tools are actually good for. 

Structure

  1. Get some XYZ-Data

  2. Convert XYZ to TIFF and (if needed) reproject

  3. Hillshading and Slope

  4. Combine Hillshade and the color-relief'd Slopes with Raster Calculator

  5. Conversion for OsmAnd/Geopaparazzi/TileMill/…

  6. Avoid a click adventure and use the shell (Linux/Unix)

  7. Summary & Download: used data and the processed example

Introduction

Hillshading in OsmAndUsing an offline navigation app, is a good solution if you're having no or bad internet connection in the field. This approach is quite considerate in resources, like the battery lifetime. The State survey provides a pre-rendered Service for even the high-resolution DGM as a WMS. The Smartphone-application I use doesn't provide the possibilty to process WMS-requests. A needed workaround for this would be to set up a Tileserver handling the requests, caching and serving the XYZ-tiles. It can be achieved with Tilecache.

But the pre-rendered hillshade-tiles of the GEOSN are not as good as they could be. Steep areas just drown in black, the overall areas are just too dark. Since the OsmAnd+ App doesn't provide a proper blending-mode, but just the "regular" additive-blending with respect to the transparency. The result does not look very promising. Supporting the blending mode "Multiply" would address some of the unwanted effects, like brightening/darkening of the map. The advantage of the Multiply-blending-mode over the transparency-approach will be shown. The OsmAndd+ app is available for Android and iOS and is using the openstreetmap-data. I will not explain all of the steps made with QGIS or GDAL in detail or explain the GIS-basics.

 

Software and Tools used, Data

Linux Ubuntu with certain pre-installed utilities like "cat" and "sort"

Data sources:

1. Get the XYZ-Data and modifiy it

In this example I used the DEM-20 (DGM20) dataset of Saxony, available at the GEOSN opendata website. For smaller areas high-resolution datasets like the DGM1 or DGM2 can be processed the same way. The DGM20 dataset is served as 20x20 km tiles. To merge all the data into one processable file, I used the Linux bash-tool "cat".

$ cat *.xyz > dgm20.xyz

The GDAL-driver requires a certain sorting of the XYZ-data, documented here. Cells with same Y coordinates must be placed on consecutive lines. For a same Y coordinate value, the lines in the dataset must be organized by increasing X values. If not properly sorted, loading the ASCII-Grid into QGIS is resulting in an error:

„Ungridded dataset: At line 123, change of Y direction“

So use the bash-tool "sort" with the follwing command.

$ sort -nk 2 -nk 1 dgm20.xyz > dgm20_sorted.xyz

If there are multiple values for XY-coordinate like

…
442020 5666000 307.26
442020 5666000 307.30
…

there will be an error:

Ungridded dataset: At line 123456789, X spacing was 0.000000. Expected > 0 value.

So some look-ahead regex would be neccessary. Didn't solved it. Just take the already merged and sorted DGM20 fulldataset which I concatenated from an early version of the DGM20-datasets, when no overlapping was in them .

Considering the amount of data, this approach is the probably the most effective way, at least it doesn't take minutes or even hours to sort the data. The resulting file is about 1 Gigabyte in size (packed 87 MB here), even it is "just" a ASCII-file. That is one reason not to concatenate all of the available high resolution datasets of the DGM2 or even the DGM1. One ASCII file of the DGM1, covering an area of 2×2 km² is already 90 MB in size. You can imagine how big the handling raster-files would become. You can load the dgm20_sorted.xyz into QGIS as a raster and style it with the "hillshading"-option. But I prefer to convert the XYZ-data with the "gdal_translate" into TIFF, before using it with QGIS.

2. Convert XYZ to TIFF and (if needed) reproject

In short, create a DEM-TIFF, calculate the slope and hillshading and compress the files afterwards. The "target" projection (mbtiles/sqlitedb) is EPSG:3857.

The reprojection of the EPSG:25833-data to EPSG:3857 requires a intermediate format. It is not possible to reproject the XYZ-ASCII (EPSG:25833) into EPSG:3857 out of the box.The calculation of the slope seems to take an obscene amount off time, when compressing the DEM-data before. A resolution (meters per pixel) of 5 for the DGM20 seems to be sufficient. The reprojection to EPSG:3857 before calculating the hillshade is needed to avoid some ugly, artifact-like squares in the map.

3. Hillshading and Slope

The hillshading for an big area with many different reliefs can be annoying. Either the mountains are quite good shaded, but areas with a flat relief are not well shown. Focusing on flat areas can resulting in a drowning in black for the steeper mountain areas. Additionally a mountain will only be shaded on one side. The other side is just bright and you actually get no information on what is there. Either you add a "weighted" additionally light source to it or you are using the calculated slope, where steep areas just appear to be darker, not regarding the light source(s). If you combine it via the blending mode "multiply", you can gain relief information back.

Another approach to handle this, is adding a multi-directional shading where there is not only one light source in the NW (315°). But it can happen that the "feeling" for what is actually a depression and which areas are a elevation can be lost. If you're not anticipating it, it can be quite confusing. The same problem appears when your viewing slopeshade maps. The slope (0 to 90°) itself gives you no information about a depression or a elevation in relation to the surrounding area. It's just anticipating.

The table shows a comparison of different approaches: „simple“, enhanced Hillshading with 50 % transparency, blending mode „multiply“ and the combination of calculated slopes and the hillshades with multidirectional lightsources.

OSM without Hillshading
Just the OpenStreetMap-View
OSM without Hillshading
simple Hillshading (HS) with 50% Transparency
OSM without Hillshading
HS with blending mode "Multiply"
OSM without Hillshading
Slope with 50% Transparency
OSM without Hillshading
slope, Multiply
OSM without Hillshading
Combination of slope and HS, Multiply
OSM without Hillshading
slope+HS, Multiply, decreased transparency
OSM without Hillshading
3 light sources for the HS, calculated slopes (Multiply)

Calculating the slope can be done with the gdal-library. The original data is in EPSG:25833, UTM-Grid (all values in metres). That's why we calculate the slope FIRST and reproject the raster afterwards. The hillshading should be done AFTER the reprojection to EPSG:3857 or you get some weird and ugly artifacts (blocks and lines). It is not that bad when using the resampling method bilinear.

gdaldem slope -of GTiff -b 1 -s 1.0 -compute_edges inputdemdata-epsg25833.tif outputslope.tif

Then assign a color map with color-relief and a colorramp-file.

gdaldem color-relief -alpha outputslope.tif -co compress=deflate -co zlevel=6 color_slope.txt output-slopeshade.tif

The color_slope.txt looks like this:

nv 255 255 255 0
0 255 255 255 0
5 255 255 255 0
90 1 1 1 255

The first column is the value for the slope angle (0° to 90°). The 2nd to 4th values are for red, green, blue (RGB) and the last value within one row is the value for the alpha channel (transparency), where 0 is fully transparent and 255 is opaque. The term "nv" stands for nodata value and will be assigned to pure white with full transparency. Same with areas where the slope does have an angle between 0 and 5 °. Such flat areas are actually not that interesting in a DGM20 - or at least not for our purpose now.

The darkest value will be 1 1 1 255, not pure black (like 0 0 0 255) since we might need this value for nodata values in the process later. Either in QGIS or with image-processing software. The output is a TIFF with grayscale band and second band - the alpha band.

Generating a hillshade from the elevation data can be done within QGIS, or by "gdaldem hillshade" (what I prefer).

gdaldem hillshade -z 1 -combined -compute_edges -alt 45 -co compress=deflate -co zlevel=6 inputdemdata-epsg3857.tif output-hillshade.tif

4. Use Raster Calculator to combine Hillshade and the color-relief'd Slopes

The result of this hillshading-process is a grayscale image with pixelvalues ranging from 1 to 254 (1 band only). Now comes the tricky part. Combining the calculated slopeshade and hillshade images. It can be done in QGIS by using the blendding mode "multiply", but using the gdal raster calculator seems to be more efficient when handling similar datasets more than once on a regular base.

gdal_calc.py -A output-slopeshade.tif --A_band=4 -B output-hillshade.tif --B_band=1 -C output-slopeshade.tif --C_band=1 --outfile=slope-hill-shades-combined.tif --calc="round_(((1-A/255.0)*B + (A/255.0)*0.5*C),0)" --NoDataValue="255" --type=Byte

We use the 2-band slopeshade raster and the 1-band hillshade image and "merge" (blend/combine/stack/whatever) it together with the raster calculator. It can be achieved with virtual rasters (vrt and gdalbuildvrt) as well, but I didn't manage to do it up to this time. The formula has been taken from here. It's a conversion from RGBA (output-slopeshade.tif band 1 "A" and it's 2. band - the alpha band "C") to it's corresponding RGB value. There will be no transparency band after the calculation. The "background" raster will be the output-hillshade.tif ("B"). It took me a while to figure out, but later it is important to set the output file-type to "Byte".

Take into account that we calculate with float values, not just integer values. (255.00 vs. 255). It seems odd to calculate within the RGB (0…255) range with float values, but the result wouldn't be the same if not using "255.0" (or "255.") float-values. Round the resulting pixel values with "round_(A,0) back to integer range. Otherwise gdal2mbtiles will not work and show just a noisy raster.

At this time OsmAnd (it's plugin) doesn't support blending modes multiply or overlay (a combination of multiply and screen blending), so the layer will be set semi-transparent. That's resulting in a pretty bright map, where city names become lighter than usual. I don't like that. That why we do some last calculation, where we create some kind of alpha mask on the combined hillshade-slopeshade raster with gdaldem color-relief and a modified color-ramp (combinealpha.txt)

0 0 0 0 255
32 0 0 0 240
64 0 0 0 220
96 0 0 0 180
128 0 0 0 150
250 0 0 0 0
255 255 255 255 0

And using this gdaldem-command:

gdaldem color-relief slope-hill-shades-combined.tif -co compress=deflate -co zlevel=6 combinealpha.txt -alpha dgm20-hillshade-final.tif

And processing finished. The result looks like this in QGIS. The PNG is transparent, so no multiply-mode necessary.

transparent Hillshading

5. Convert to MBTILES/SQLITEDB for OsmAnd/Geopaparazzi/TileMill/…

Converting it to MBTILES within QGIS or gdal2mbtiles or gdal2mbtiles.py and the processed raster is almost ready for OsmAnd Just convert it to sqlitedb with mbtiles2osmand.py and name it

Hillshade Germany Saxony DGM20.sqlitedb

and put it into the OsmAnd-folder

… Android/data/net.osmand.plus/files/tiles/

Seems to be undocumented, but the first part of the filename should always be "Hillshade", so the OsmAnd-plugin will load it by default.

6. Avoid a click adventure and use the shell — Sorry, Linux/Unix only

Here some bash-code for all the stuff we've done after sorting the data.
Copy and paste. Just put the color_slope.txt, combinealpha.txt and the mbtiles2osmand.py in the same folder where the dgm20.xyz (sorted) is located. Install gdal2mbtiles or modify section 11.

clear &&
inputfile=dgm20.xyz &&
outputfilename=dgm20 &&
resolution="5 5" &&
sourcesrs=EPSG:25833 &&
targetsrs=EPSG:3857 &&
rampfile=color_slope.txt &&
alphamask=combinealpha.txt &&
minzoom=6 &&
maxzoom=14 &&
png8colors=32 &&

echo -e "Create nice looking hillshading map with GDAL and convert them to MBTILES. The SQLITE-Files can be used in Apps like OsmAnd+ (SQLITEDDB) or Geopaparazzi (MBTILES - atm jpg only) \n
Inputfile: $inputfile \n
Outputfiles contain \"$outputfilename\" \n
X,Y Resolution: $resolution (units of the SRS) \n
Inputfile-SRS: $sourcesrs \n
Target SRS: $targetsrs \n
colordata (slope): $rampfile \n
colordata (alpha-masking): $alphamask \n
zoomlevels (mbtiles): $minzoom"-"$maxzoom \n\n

Files with \"temp\" contain uncompressed data and will be deleted to save storage. Using the uncompressed files for the calculation of the slope, shades and other raster operations is faster than using the (lossless) compressed data.\n $alphamask"," $rampfile"," mbtiles2osmand.py needs to be in the same folder as \"$inputfile\". gdal2mbtiles needed aswell \n

Abort now to modify the parameters … \n" &&
sleep 1 &&

echo -e "\n1. Convert XYZ to TIFF … " &&
gdal_translate -tr $resolution -r bilinear -of GTIFF -a_nodata -9999 -a_srs $sourcesrs $inputfile $outputfilename"_temp.tif" &&

echo -e "\n2. Calculate the slope using the uncompressed DEM-TIFF, or it would take hours … " &&
gdaldem slope -of GTiff -b 1 -s 1.0 -compute_edges -co bigtiff=if_safer $outputfilename"_temp.tif" $outputfilename"_slopetemp.tif" &&

echo -e "\n3. Reproject the DEM-TIFF and SLOPE-TIFF from $sourcesrs to ${targetsrs} . You might need to add the -srcnodata and -dstnodata parameter in the gdalwarp command…" &&
gdalwarp --config GDAL_CACHEMAX 1024 -wm 500 --config NUM_THREADS ALL_CPUS -s_srs $sourcesrs -t_srs $targetsrs -tr $resolution -r bilinear $outputfilename"_slopetemp.tif" $outputfilename-$targetsrs"_slopetemp.tif" &&
gdalwarp --config GDAL_CACHEMAX 1024 -wm 500 --config NUM_THREADS ALL_CPUS -s_srs $sourcesrs -t_srs $targetsrs -tr $resolution -srcnodata -9999 -dstnodata -9999 -r bilinear $outputfilename"_temp.tif" $outputfilename-$targetsrs"_temp.tif" &&

echo -e "\n4. Compress DEM-TIFF and SLOPE-TIFF in target-projection… " &&
gdal_translate --config GDAL_CACHEMAX 1024 --config NUM_THREADS ALL_CPUS -co compress=deflate -co zlevel=9 -co tiled=yes -co blockxsize=2048 -co blockysize=2048 -co tfw=no -co interleave=pixel -co predictor=3 -co bigtiff=if_safer -co nbits=32 -ot float32 -tr $resolution -r bilinear -of GTIFF -a_srs $targetsrs $outputfilename-$targetsrs"_slopetemp.tif" $outputfilename-$targetsrs"_slope.tif" &&
gdal_translate --config GDAL_CACHEMAX 1024 --config NUM_THREADS ALL_CPUS -co compress=deflate -co zlevel=9 -co tiled=yes -co blockxsize=2048 -co blockysize=2048 -co tfw=no -co interleave=pixel -co predictor=3 -co bigtiff=if_safer -co nbits=32 -ot float32 -tr $resolution -r bilinear -of GTIFF -a_srs $targetsrs $outputfilename-$targetsrs"_temp.tif" $outputfilename-$targetsrs".tif" &&

echo -e "\n5. Create Pyramid-files (ovr) for fast handling in GIS … " &&
gdaladdo -ro --config COMPRESS_OVERVIEW DEFLATE $outputfilename-$targetsrs"_slope.tif" 2 4 8 16 && gdaladdo -ro --config COMPRESS_OVERVIEW DEFLATE $outputfilename-$targetsrs".tif" 2 4 8 16 &&

echo -e "\n6. Create a hillshade and slopeshade (see txt-files) …" &&
gdaldem color-relief -alpha -co bigtiff=if_safer $outputfilename-$targetsrs"_slopetemp.tif" $rampfile $outputfilename-$targetsrs"_slopeshade_temp.tif" &&
gdaldem hillshade -z 3 -combined -compute_edges -alpha -alt 45 $outputfilename-$targetsrs"_temp.tif" $outputfilename-$targetsrs"_shade_temp.tif" &&

echo -e"\n7. Compress the calculated files …" &&
gdal_translate --config GDAL_CACHEMAX 1024 --config NUM_THREADS ALL_CPUS -co compress=deflate -co zlevel=6 -co tiled=yes -co blockxsize=2048 -co blockysize=2048 -co tfw=no -co interleave=pixel -co predictor=2 -co bigtiff=if_safer -r bilinear -of GTIFF -a_srs $targetsrs $outputfilename-$targetsrs"_slopeshade_temp.tif" $outputfilename-$targetsrs"_slopeshade.tif" &&
gdal_translate --config GDAL_CACHEMAX 1024 --config NUM_THREADS ALL_CPUS -co compress=deflate -co zlevel=6 -co tiled=yes -co blockxsize=2048 -co blockysize=2048 -co tfw=no -co interleave=pixel -co predictor=2 -co bigtiff=if_safer -r bilinear -of GTIFF -a_srs $targetsrs $outputfilename-$targetsrs"_shade_temp.tif" $outputfilename-$targetsrs"_shade.tif" &&

echo -e "\n8. Calculate a combination of Slope and the casual Hillshade (RGBA to RGB conversion with a 0.5 factor) … "&&
gdal_calc.py -A $outputfilename-$targetsrs"_slopeshade_temp.tif" --A_band=4 -B $outputfilename-$targetsrs"_shade_temp.tif" --B_band=1 -C $outputfilename-$targetsrs"_slopeshade_temp.tif" --C_band=1 --outfile=$outputfilename"-combined_temp.tif" --calc="round_(((1-A/255.0)*B + (A/255.0)*C*0.5),0)" --NoDataValue="255" --type=Byte &&

echo -e "\n9. \"Thinning\" - Create a alpha-mask from combinealpha.txt (make the brightest and flatest areas transparent)" && gdaldem color-relief $outputfilename"-combined_temp.tif" -co compress=deflate -co zlevel=6 $alphamask -alpha $outputfilename"-hillshade-final.tif" &&

echo -e "\n10. Delete the uncompressed \"temp\"-files … manually or remove the echo command before rm \n" &&
echo rm $outputfilename-$targetsrs"_temp.tif" $outputfilename"_temp.tif" $outputfilename"_slopetemp.tif" $outputfilename-$targetsrs"_slopetemp.tif" $outputfilename-$targetsrs"_shade_temp.tif" $outputfilename-$targetsrs"_slopeshade_temp.tif" $outputfilename"-combined_temp.tif" &&

echo -e "\n11. Create MBTILES for configured zoomlevels. Grayscale png8 with $png8colors colors for quality-filesize improvement. Very big TIFFs may need to be uncompressed for processing in gdal2mbtiles. If diskspace is short temporarily set: \$ export TMPDIR=/home/<username>/tmp" &&
gdal2mbtiles --min-resolution $minzoom --max-resolution $maxzoom --no-fill-borders --png8 $png8colors --name kreidefossilien.de --resampling bilinear --layer-type overlay $outputfilename"-hillshade-final.tif" &&

echo -e "\n12. Convert the mbtiles to sqlitedb … " &&
python3 mbtiles2osmand.py $outputfilename"-hillshade-final.mbtiles" $outputfilename"-hillshade-final.sqlitedb" &&

echo -e "\n\n… Done. \n Certain intermediate files has not been deleted by this script. Keep the one you need and delete the rest by yourself (copy the output at section 10) ."

Some resources on hillshading, raster calculation and the other stuff.

 

7. Summary & Download: used data and the processed example(s)

The following 4 pages may interest you as well:

Comments (0)






Allowed tags: Add a new comment: