Adding support for Swissalti3d elevation (opendata)

Open Source Routing Engine
1

Hi,

Swissalti3d is a precise elevation data (0.5m or 2m step) available freely as opendata.
I have started investigating how to add support for it in Graphhopper. I would be glad to have feedback and maintainer guidance on how to do that in a great way.

If there is interest for it by the maintainers, my plan is to have it merged into graphhopper mainstream.
Otherwise it would stay in a fork, available for anyone interested in it.

Here are some insights / details following my primary investigation.

The datasource: opendata Swissalti3d

Area: Switzerland

Cloud Optimized Geotiff (COG), LZW compression
0.5 m COG: ~ 26 MB / tile, 770 GB / full coverage
2 m COG: ~ 1 MB / tile, 44 GB / full coverage

Projection: LV95 (EPSG:2056)

Downloading all tiles:

Select entire dataset + download csv with the URLS for all the tiles

Example of tiles from the CSV (these are 2m steps):
swissalti3d_2019_2501-1120_2_2056_5728.tif

Gdal output:
Driver: GTiff/GeoTIFF
Files: swissalti3d_2019_2501-1120_2_2056_5728.tif
swissalti3d_2019_2501-1120_2_2056_5728.tif.aux.xml
Size is 500, 500
Coordinate System is:
PROJCRS[“CH1903+ / LV95”,
…
USAGE[
SCOPE[“Cadastre, engineering survey, topographic mapping (large and medium scale).”],
AREA[“Liechtenstein; Switzerland.”],
BBOX[45.82,5.96,47.81,10.49]],
ID[“EPSG”,2056]]
Data axis to CRS axis mapping: 1,2
Origin = (2501000.000000000000000,1121000.000000000000000)
Pixel Size = (2.000000000000000,-2.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
LAYOUT=COG
Corner Coordinates:
Upper Left ( 2501000.000, 1121000.000) ( 6d 9’21.66"E, 46d14’ 3.78"N)
Lower Left ( 2501000.000, 1120000.000) ( 6d 9’22.42"E, 46d13’31.40"N)
Upper Right ( 2502000.000, 1121000.000) ( 6d10’ 8.32"E, 46d14’ 4.31"N)
Lower Right ( 2502000.000, 1120000.000) ( 6d10’ 9.08"E, 46d13’31.93"N)
Center ( 2501500.000, 1120500.000) ( 6d 9’45.37"E, 46d13’47.86"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
Min=372.139 Max=372.141
Minimum=372.139, Maximum=372.141, Mean=372.140, StdDev=0.000
NoData Value=-9999
Overviews: 250x250

Steps:

  • figure out the tiling naming conventions of Swissalti3d
  • create a Swissalti3dElevationProvider on the model of CGIARProvider

Pain points?

I am open to your feedback.
If someone want to directly participate, also let me know.

Guillaume

2

create a Swissalti3dElevationProvider on the model of CGIARProvider

A better approach might be first converting the Swissalti tiles to a format that GraphHopper understands already. Actually, I think I did this before using gdal (or QGIS? or both?), I’ll try to remember, but I think it worked quite well once I had figured out the right command.

3

Thanks @easbar for your early feedback.
Indeed, your approach would minimize the changes to do in Graphhopper (eventually, there would be no change a all!). If I publish the data pipeline tool, others could use it to prepare there data.

I will continue with my initial approach today and see what I can achieve until Monday.

My main issue currently is that the TiffDecoder used in CGIARProvider does not support the tiff from Swisstopo: the sample format is not supported.
org/apache/xmlgraphics/image/codec/tiff/TIFFImage.java

            case 32:
              if (sampleFormat[0] == 3) {
                isValidDataFormat = false;
            }
4

I was able to read the Swisstopo tiff file by using standard ImageIO instead of the org.apache.xmlgraphics codec (which is used in CgiarProvider) and upgrading the compiler from 1.8 to 1.17 (I have not tried intermediate versions).

What is the minimal SDK version for graphhopper nowaday? The pom file states it is java 1.8 but I remember there were some discussions to use more modern baseline.

The CgiarProvider is using some MMAPed heights DataAccess structure + a map in main memory. I have GBytes of elevation in my case. Is it something I need to do? Or can I simply keep in main memory the last used Rasters?

Thanks in advance for your advices.

5

JDK 8

Probably not, but you need to anticipate that height data from different tiles will be fetched randomly. If the file loading is somehow too slow it might be worth to optimize this, but otherwise there is no real reason you need to do anything special. GraphHopper more or less simply reads one OSM node after the other and fetches the associated elevation values. To make sure each tile is only loaded once, or a few times you would have to somehow order the nodes by their geographic location.

Or can I simply keep in main memory the last used Rasters?

Yes, this might be sufficient, but since the node positions can be very random (?!) maybe this won’t really speed up anything. You can probably start without any such optimizations and later decide if it needs to be improved.

Another approach you could take: Modify the OSM file and assign an elevation tag for each node (something like gh:ele=651). Then use this in a very simple ElevationProvider that implements double getEle(ReaderNode). I’m not sure what kind of resolution you are looking at. The current elevation providers work with something like a ~30x30m or ~90x90m resolution (1 pixel of the elevation tile corresponds to an area of ~30x30m (the exact size of the pixel also depends on the coordinates)). But if are using elevation data at a higher resolution enriching the OSM file like this could be beneficial.

6

Hi, thank you @easbar for your inputs.

I created an early draft in my fork: Add swissalti3d elevation provider (WIP) by gberaudo · Pull Request #1 · gberaudo/graphhopper · GitHub

The raster precision of the data is really high (2x2m, ~0.5m in the z direction in the good cases). I want to keep that precision with graphhoppper.

Indeed, it would be great if graphopper was iterating the nodes in a geographically aware way. I doubt I will have time to look into that, but it sounds useful for the existing and future elevation providers.

I will now look at integrating proj4j and test on a real scenario.

7

I was able to add the small proj4j dependency and test with a real life case (an OSM extract of a montagneous area in Switzerland). Regarding caching, I added a basic mechanism to keep only the last 100 tiles in memory. In addition, I configured sampling so that 3d points are added along the ways.

This did work fine, consumming directly the elevation tiles from the remote opendata source. Cool!

I don’t know the impact of sampling on the generated network though. Does it have an impact on routing performance? On the other hand it allows quite accurate routing taking elevation into account. Nice for cycling!

  # To increase elevation profile resolution, use the following two parameters to tune the extra resolution you need
  # against the additional storage space used for edge geometries. You should enable bilinear interpolation when using
  # these features (see #1953 for details).
  # - first, set the distance (in meters) at which elevation samples should be taken on long edges
  # graph.elevation.long_edge_sampling_distance: 2
  # - second, set the elevation tolerance (in meters) to use when simplifying polylines since the default ignores
  #   elevation and will remove the extra points that long edge sampling added
  # graph.elevation.way_point_max_distance: 1
8

This topic was automatically closed 90 days after the last reply. New replies are no longer allowed.