Elevation raster generation (*.bef files)

afischerdev commented 8 months ago

According to the issues #558 and #539 I added some new classes with new naming conventionElevationRaster. This generates now also bef raster files in 1sec format from hgt formatted tiles in 1sec raster. The 1sec 5x5 raster generation has a fallback, when no 1sec source is found, a 3sec hgt file is called.

As a fallback the ESRI asc format is also supported. But please note this format has a 6000x6000 raster, the 3sec from hgt source have a size of 6001x6001 points.

Arguments for single file generation:

ElevationRasterTileConverter <srtm-filename | all> <hgt-data-dir> <srtm-output-dir> [arc seconds (1 or 3,default=3)] [hgt-fallback-data-dir]

Samples

$ ... ElevationRasterTileConverter srtm_34_-1 ./srtm/hgt3sec ./srtm/srtm3_bef
$ ... ElevationRasterTileConverter srtm_34_-1 ./srtm/hgt1sec ./srtm/srtm1_bef 1
$ ... ElevationRasterTileConverter srtm_34_-1 ./srtm/hgt1sec ./srtm/srtm1_bef 1 ./srtm/hgt3sec

Arguments for multi file generation (world wide):

$ ... ElevationRasterTileConverter all ./srtm/hgt3sec ./srtm/srtm3_bef
$ ... ElevationRasterTileConverter all ./srtm/hgt1sec ./srtm/srtm1_bef 1
$ ... ElevationRasterTileConverter all ./srtm/hgt1sec ./srtm/srtm1_bef 1 ./srtm/hgt3sec

To work with the new raster size the rd5 generatator needs an extra paramater for the 3sec fallback

e.g. $ ... PosUnifier nodes55 unodes55 bordernids.dat bordernodes.dat ../srtm/srtm1_bef ../srtm/srtm3_bef

As we have no hgt data on server I generated the Europa 1sec bef files locally from Sonny source and copied it to server into /srtm1_bef. There is already a folder /srtm3_bef that contains all the current elevation raster files.

The source already contains the routines for a 5x5 raster naming as it is used for the rd5 tiles. But as we have no world wide hgt data on server, we can't use it at the moment.

zod commented 8 months ago

Thanks for your work, as I only proposed the changes but currently lack the time to do any work. I hope I can find some time to review your changes this week.

I remember @abrensch stating that he didn't find improvements in elevation accuracy by using the 1asec files compared to 3asec. Did you notice any improvements?

I think the most sensible way is to use LIDAR data where it's available (hopefully most of europe) and SRTM data for all other areas.

afischerdev commented 8 months ago

@zod

I hope I can find some time to review your changes this week.

We have all the time we need. There already is an elevation pool.

Did you notice any improvements?

No, I didn't. But I used Iceland for my testing - is nice and isolated. But it has no old values. This reminds me, I wanted to do tests with tunnel entrances. But this is the only idea for a quality check at the moment.

use LIDAR data where it's available (hopefully most of europe)

Yes, that is the plan.

quaelnix commented 7 months ago

I think the most sensible way is to use LIDAR data where it's available and SRTM data for all other areas.

Yes, that is the plan.

I would love to see this, but how do we solve the problem that the current elevation filter logic won't work with both: https://github.com/abrensch/brouter/blob/ef73d468c057fc72208b5a32756ce73cda92c780/brouter-core/src/main/java/btools/router/RoutingEngine.java#L840-L845 The SRTM data requires a threshold of 10 m while the LIDAR data apparently requires a threshold of 5 m to get resonable results in the total ascent calculation. I guess we would have to come up with a new filter method that works reasonably well with both data sets?

afischerdev commented 7 months ago

@quaelnix Thanks for that pointer. I'll see where is the way to get that info forward to routing .

a threshold of 5 m

How did you get this value?

quaelnix commented 7 months ago

By comparing the calculated total ascent of several routes with the expected total ascent (current BRouter, Komoot, Garmin).

afischerdev commented 7 months ago

Here comes a workaround for a filter logic needed in RoutingEngine - please see above. The information is passed on via the file name over 1 or 3 sec generation. Not very elegant, I think, but should work. The rd5 tiles then have a feature about the generation and this can then be used to fill the necessary filter with suitable values.

quaelnix commented 7 months ago

When I was playing around with Sonny's LiDAR data back in June I found that the 3asec outperformed the 1asec data.

Either I did something wrong when generating the *.bef files, or the 1asec data is simply more prone to artifacts.

Would it be possible to provide a download link for the .rd5 files from Germany that contain the 1asec LiDAR data?

afischerdev commented 7 months ago

@quaelnix Sorry, there are no rd5 files with 1sec around. But I could provide the four *.bef files for the generation.

quaelnix commented 7 months ago

I guess I must have done something wrong because the 1asec data looks fabulous: brouter_lidar_1asec_sub brouter_lidar_3asec_sub

abrensch commented 4 months ago

I'm coming back to the subject "BEF generation" as part of my cleanup-project (branch "cleanupbranch")

Meanwhile I adapted all binary-encoding related stuff to use the "statcoding" library.

Unfortunatly, that used the reverse bit order compared to the current encoding code, and that affects also BEF generation.

So I adapted the BEF format (purely technically, just reverse bit order and a magic-header word to detect wrong format) and at the same time I adapted the naming scheme (srtm_38_03.bef -> E5_N45.bef )

I also converted the tiny befs in the unit test that way.

Then I converted the original SRTM3 set to that new format ( /home/mapcrator/srtm3_bef )

But I feel that is the right time to do the full job (=multi source merge from lidar, srtm1, srtm3 to create a global 1-arc-second data set).

So I would like to understand the current status of research and data availability

thanx for any pointers,

Arndt

afischerdev commented 4 months ago

@abrensch Please see my repo map-creator. The three files starting with Elevation do the job. And at the end PosUnifier uses the new logic with 1sec and fallback to 3sec

If you want me I could do the transfer to the reverse bit order - but it will take some time. In the past I generated the 1sec bef localy and transfered it to the server.

nrenner commented 4 months ago

There is "NASADEM Merged DEM Global 1 arc second", released in 2020:

NASADEM extends the legacy of the Shuttle Radar Topography Mission (SRTM) by improving the digital elevation model (DEM) height accuracy and data coverage as well as providing additional SRTM radar-related data products. The improvements were achieved by reprocessing the original SRTM radar signal data and telemetry data with updated algorithms and auxiliary data not available at the time of the original SRTM processing.

https://lpdaac.usgs.gov/news/release-nasadem-data-products/

nrenner commented 4 months ago

Another interesting global dataset instead of SRTM/NASADEM might be:

Copernicus DEM (GLO-30 Public)

https://spacedata.copernicus.eu/collections/copernicus-digital-elevation-model

"GLO-30 offers global coverage at a resolution of 30 metres"
"Data were acquired through the TanDEM-X mission between 2011 and 2015"
"The GLO-30 dataset is available worldwide with a free license available here."

The license link is broken, but I found the document with a "_latest" added:

broken: CSCDA_ESA_Mission-specific+Annex_31_Oct_22.pdf
found: CSCDA_ESA_Mission-specific+Annex_31_Oct_22_latest.pdf
- see "Licence for Copernicus DEM instance COP-DEM-GLO-30-F Global 30m Full, Free & Open" on last pages 21-23 (looks good to me)
and also INFO/eula_F.pdf in download below

Example:

wget https://prism-dem-open.copernicus.eu/pd-desk-open-access/prismDownload/COP-DEM_GLO-30-DTED__2023_1/Copernicus_DSM_10_N47_00_E009_00.tar
tar -xf Copernicus_DSM_10_N47_00_E009_00.tar

using DTED format here, because DGED would need adoption for 32 bit float values
for different locations replace N47_00_E009_00, see any 1° grid e.g. http://dwtkns.com/srtm30m/

convert to HGT

gdal_translate -of SRTMHGT Copernicus_DSM_10_N47_00_E009_00/DEM/Copernicus_DSM_10_N47_00_E009_00_DEM.dt2 N47E009.HGT

or convert to .bil.zip

mkdir tmp
gdal_translate -of EHdr Copernicus_DSM_10_N47_00_E009_00/DEM/Copernicus_DSM_10_N47_00_E009_00_DEM.dt2 tmp/Copernicus_DSM_10_N47_00_E009_00_DEM.bil
zip -j Copernicus_DSM_10_N47_00_E009_00_DEM.bil.zip tmp/*
rm -r tmp

(drag&drop onto leaflet-raw-dem viewer to inspect)

nrenner commented 4 months ago

Regarding DTMs from LiDAR: apart from the Sonny compilation, I'm only aware of the older Terrain Tiles below. Might be worth doing some research and asking around.

Terrain Tiles (Skadi)

https://registry.opendata.aws/terrain-tiles/

Skadi format: 1° SRTMGL1 format tiles but with global coverage and compressed using gzip
- (BRouter-Web uses the Terrarium tiles for the Hillshading layer)
released 2016-2017, unmaintained
Data sources, probably a few useful LiDAR DTMs outside Europe, e.g. USA
- (but maybe not Canada, see joerd#189)
- misses newer sources, e.g. compared to OpenDemSearcher

Example: https://s3.amazonaws.com/elevation-tiles-prod/skadi/N40/N40W075.hgt.gz

abrensch commented 4 months ago

@abrensch Please see my repo map-creator. The three files starting with Elevation do the job. And at the end PosUnifier uses the new logic with 1sec and fallback to 3sec

thanx. I will merge it manually to the cleanup branch.

However, I like the idea to have only one single full 5x5 file in BEF format at map creation time. BEF format is made for fast decoding, and for 1 arc second data that matters. Different handling per raster type at mapcreation time should not ne neccesary, because the BEF files contain all meta infos needed.

I remember that last time I tried that (some years ago) I had a hard struggle to create such files with an "extended border" of some pixellines, in order to allow bigger interpolation areas also at the tile-border.

No I think, because it is so difficult, creating such a border should be a seperate processing step working only on 5x5 BEF files.

The links on the availability of Copernicus data looks promising. I checked for it at one time or the other, but my status was you must apply for access as a research institute.

nrenner commented 4 months ago

The links on the availability of Copernicus data looks promising. I checked for it at one time or the other, but my status was you must apply for access as a research institute.

That was my status as well, I just discovered that it's open now when browsing the elevation tag on the AWS open data registry.

quaelnix commented 4 months ago

@afischerdev, would it be possible to rebase this branch on master and then perhaps create a https://brouter.de/afischerdev test instance, with .rd5 files containing the LIDAR (1asec) elevation data?

afischerdev commented 4 months ago

@quaelnix The test server is https://brouter.de/essbee I could generate some test area (e.g. Germany) and we switch the segment path for a while. Special interests?

quaelnix commented 4 months ago

Special interests?

I would be very interested in a test instance with Sonny's LiDAR 1asec data. Especially after @rkflx mentioned a few days ago that it was quite easy to find routes where this dataset fails, which does not match my previous observations.

afischerdev commented 4 months ago

@quaelnix Yes, I understand that, the question was about interests in an area to create. I will not generate a fully set of Europe.

quaelnix commented 4 months ago

the question was about interests in an area to create

Germany.

afischerdev commented 4 months ago

@quaelnix Sorry, I'm not able to generate a full Germany data set. The lua script is running more then 24 hours. Several config infos about the postgres database didn't help. So the test data will come without estimated_*_classes.

Added an image creator for bef and hgt files. Calls:

java -cp %CLASSPATH% btools.mapcreator.CreateElevationRasterImage -15 65 ./srtm/srtm3_bef_iceland testimg3.png 1200 1200 1
java -cp %CLASSPATH% btools.mapcreator.CreateElevationRasterImage -15 65 ./srtm/hgt3sec_iceland testimg3_hgt.png 1200 1200 1 hgt
java -cp %CLASSPATH% btools.mapcreator.CreateElevationRasterImage -15 65 ./srtm/srtm1_bef_iceland testimg1.png 3600 3600 1
java -cp %CLASSPATH% btools.mapcreator.CreateElevationRasterImage -15 65 ./srtm/hgt1sec_iceland testimg1_hgt.png 3600 3600 1 hgt

To compare the images I used WinMerge.

afischerdev commented 4 months ago

For more test choices: added a color table to define output. With the inbuild color table and a tile in the Netherlands you will get only green areas. A new color table could look like this:

# to elev, r, g, b
-20, 0, 0, 255
-1, 128, 0, 255
0, 102, 153, 153
1, 0, 102, 0
30, 251, 255, 128
100, 224, 108, 31

Result compare with WinMerge 3sec bef alt and viewfinder (I didn't find a N52E004 file in Sonny data). hgt_compare_nl

Call for this

java -cp %CLASSPATH% btools.mapcreator.CreateElevationRasterImage 4 52 ./srtm/srtm3_bef testimg3.png 1200 1200 1 bef colors_nl.csv
java -cp %CLASSPATH% btools.mapcreator.CreateElevationRasterImage 4 52 ./srtm/hgt3sec testimg3_hgt.png 1200 1200 1 hgt colors_nl.csv

afischerdev commented 4 months ago

The test server is running - with German data only.

Please try this from #670

Please let us know when to change the data back to standard.

EssBee59 commented 4 months ago

The test server is running - with German data only.

Please try this from #670

Please let us know when to change the data back to standard.

Hello afischerdev!

I started some tests and the first results are very promising: as example in that location the elevation data looks much better on the test as on the production server!

https://brouter.de/brouter-test/#map=14/49.9931/8.9207/osm-mapnik-german_style&lonlats=8.892191,50.004407;8.948779,49.97002 Test: Production:

Note: be carrefull using the test-system: not all the tags are available (estimated_traffic_class a. e.)

EssBee59 commented 3 months ago

This case is strange for me as the elevation costs differ by more than 100%...and the routing on the test system is longer and higher... "normal" routing: https://brouter.de/brouter-web/#map=14/50.9043/10.6361/osm-mapnik-german_style&lonlats=10.623796,50.907014;10.684802,50.931672

"forced routing" to the same route as "test": https://brouter.de/brouter-web/#map=15/50.9130/10.6639/osm-mapnik-german_style&lonlats=10.623796,50.907014;10.667346,50.921202;10.684802,50.931672 Elev cost = 864

Test-system: https://brouter.de/brouter-test/#map=13/50.9128/10.6639/osm-mapnik-german_style&lonlats=10.623796,50.907014;10.684802,50.931672 Elev cost = 289

quaelnix commented 3 months ago

This case is strange for me as the elevation costs differ by more than 100%

If you decrease the downhillcutoff to 0.75 you get a similar elevation cost.

the routing on the test system is longer and higher...

The Ascend | Plain ascend which is ultimately displayed in BRouter-Web is calculated independently of the elevation cost, which is fed into the route cost. The filter logic is completely different.

quaelnix commented 3 months ago

@EssBee59, also keep in mind that placing waypoints will affect the content of the elevation filter around the waypoint.

EssBee59 commented 3 months ago

I just liked to understand, why the elevation costs here differ so much...

same profile
same route

https://brouter.de/brouter-test/#map=13/50.9129/10.6639/osm-mapnik-german_style&lonlats=10.657338,50.902141;10.684802,50.931672 (elev costs 289)

https://brouter.de/brouter-web/#map=13/50.9129/10.6640/osm-mapnik-german_style&lonlats=10.657338,50.902141;10.684802,50.931672 (elev cost 884)

I could not find why the costs so differ

abrensch / brouter