Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Elevation raster generation (*.bef files) #644

Merged
merged 15 commits into from Apr 3, 2024

Conversation

afischerdev
Copy link
Collaborator

According to the issues #558 and #539 I added some new classes with new naming convention ElevationRaster.
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
Copy link
Collaborator

zod commented Nov 7, 2023

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
Copy link
Collaborator Author

@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
Copy link
Collaborator

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:

if (ehb > 0) {
ascend += ehb;
ehb = 0;
} else if (ehb < -10) {
ehb = -10;
}

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
Copy link
Collaborator Author

@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
Copy link
Collaborator

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

@afischerdev
Copy link
Collaborator Author

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
Copy link
Collaborator

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
Copy link
Collaborator Author

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

@quaelnix
Copy link
Collaborator

quaelnix commented Nov 29, 2023

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

@abrensch
Copy link
Owner

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
Copy link
Collaborator Author

@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
Copy link
Contributor

nrenner commented Feb 21, 2024

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
Copy link
Contributor

nrenner commented Feb 22, 2024

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:

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
    

@nrenner
Copy link
Contributor

nrenner commented Feb 23, 2024

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
  • released 2016-2017, unmaintained
  • Data sources, probably a few useful LiDAR DTMs outside Europe, e.g. USA

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

@abrensch
Copy link
Owner

@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
Copy link
Contributor

nrenner commented Feb 24, 2024

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
Copy link
Collaborator

@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
Copy link
Collaborator Author

@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
Copy link
Collaborator

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
Copy link
Collaborator Author

@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
Copy link
Collaborator

the question was about interests in an area to create

Germany.

@afischerdev
Copy link
Collaborator Author

@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
Copy link
Collaborator Author

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
Copy link
Collaborator Author

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
Copy link
Collaborator

EssBee59 commented Mar 7, 2024

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:
image
Production:
image

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

@EssBee59
Copy link
Collaborator

EssBee59 commented Mar 11, 2024

@quaelnix
Copy link
Collaborator

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
Copy link
Collaborator

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

@EssBee59
Copy link
Collaborator

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

@afischerdev afischerdev merged commit 6330325 into abrensch:master Apr 3, 2024
1 check passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

6 participants