More detailed DEMS #3
Comments
|
I am open to suggestions on how to make it more customizable and maybe have some inputs. for different sources usually I just make a copy of the script and tweak it, for example I had made one for GEBCO bathymetry which is basically the same with a few tweaks https://github.com/acalcutt/GEBCO_to_MBTiles/blob/main/create_terrainrgb.sh The real trouble for me is the space each layer takes. if you want to do higher resolution, each zoom layer is basically 4x bigger than the one before it. I do like the idea of using higher resolution with fallback to jaxa, but I have no idea how you would accomplish that. it seems like you would need layers or masking in the gdal vrt portion of the script, and i'm not sure that is an option. We had discussed something similar filtering by height at https://osmus.slack.com/archives/C01G3D28DAB/p1692976515173599 , and masking was a suggestion, which is similat to what this would need I think. |
|
Agree with what you're saying but see below the difference on our sample over France, this looks really great and I'm really looking being able to offer such nice rendering for the planet even with storage needs are over a TB (as long as the script can do it, we can use cloud compute power to calculate it once and make it available (torrent ?) for everyone like planetiler does it for the OSM data in mbtile)
A couple of findings :
|
|
Maybe some ideas here : https://github.com/tilezen/joerd/blob/master/docs/data-sources.md |
Now with the Webp compression I guess it opens new opportunity size wise ? Also to continue the discussion :
Another approach would be to keep the terrain data at Z12 and work on a hillshading mbtile till Z15 ? |
|
I don't think I had ever mentioned it, but I did test making a swisstopo Z0-Z17 in terrarium format just to see if I could make it It looks pretty good up untill zoom 16/17 where it starts to get a bit blocky I do wonder though if I went back to terrainrgb format and found some better defaults for base and interval if it would look better. maplibre support for other base/intervals has been added in maplibre/maplibre-style-spec#326 . I also have not yet tested how it looks in just the default mapbox terrainrgb defaults, but it sounds like they are not good for high zoom levels |
|
After that I also tried to combine TIN Italy with swisstopo with something like the following. this seemed like it was working but eventually the system froze on me so I never finished. I originally tried to use the tifs from both sets directly, but it was complaining of different projections, so I added a intermediate step of creating vrts in the projection I wanted, then combing them into a file list to combine into a combined vrt later. |
|
Yes you have to normalize projection before doing that. That was probably not necessary using only one dataset, but using multiple sources will make that step mandatory. You also have to do the dataset fusions at borders BEFORE the rendering takes place else it's generating unwanted artefacts. |
|
Jaxa has realeased a new updated version of their data (v 4.1) : https://www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/aw3d30_e.htm |
|
One way for the overlapp of DEM tiles haveing different resolution would be to have a routine to interpolate the DEMs datasets to find overlapping tiles. Since they may be with different resolution they woul need to be merged with a smoothing algorithm. |
|
Just an FYI, I did make a new JAXA 2024 webp file I was also made a new swissalti file With generating these I was thinking above the original topic of creating a merged file. I am testing something like this right now |
|
This looks very promising. When you have validated the algorithm, feel free to check if this works with https://geoservices.ign.fr/rgealti (1m resolution for France). A good area to check if this works is in the Alps at the border between France and Italy (near Val d'Isère / Grande Sassière). |
|
Used your latest JAXA 2024 webp file here : https://maps.utagawavtt.com/styles/terrain-rgb/#12.3/45.42516/6.00561 I guess you have not yet used the cubic mode as wer can see the linear artefacts between sampleng bands.
Have you also checked on implemenintg better comprssion using this ? #7 ? |
|
I just finished the test that combined the jaxa, swiss, and itally dems and did this one with 'cubic' instead. The areas with swiss/itally dems seem improved to me, however the areas with jaxa look much worse.. I am not sure if it looks worse because of the 'cubic' instead of 'cubicspline ' or if this has to do with me having to convert the jaxa dataset to float32 so it could be used with the other datasets in a VRT You can see the new jaxa/swiss/italy example I made with cubic at Compared to the new jaxa webp you just downloaded it made in cubicspline The issue i saw looks like this in the jaxa area |
|
Thanks for the test.
The new jaxa/swiss/italy is definitely showing more details so this is great from that perpsective.
But the tiles junctions are definitely showing artefacts.
@xfischer, @nbozon, @xycarto, @ThomasHalwax, @rouault any thoughts/recommendations ?
|
|
I've tried countless times to fix this. I thought it was an threshold/corner-case issue with my own resampling algorithm, but it seems more related to a general extent of the way resampling is done with DEM of different resolutions. Maybe on GDAL forums there are questions about that ? |
|
Hey @xfischer ! In some projects I just use proper maximal zoom level for a dem resolution. For instance for webMercator for 30 m jaxa, maximal proper zoom for 256x256 dem image is 12(also it depends on tile Y coordinates). Bigger zooms of the data will cause stairs artefacts. I do not handle this in the @openglobus, yet…but the solution could be some proper interpolation (gdal “trilinearsmoothaverage” made some noise on the raw data) elevation data, or using high resolution dem:) Cheers! |
|
I made a few smaller area versions of the jaxa_swiss_italy combination in 'bilinear' and 'cubicspline' with but it didn't seem to make much difference for the combined file. I made these demo pages so maybe we can make a post on the gdal mailing list. It shows the 3 countries that I had datasets for in 1 spot The first demo page is just with the jaxa dataset (webp, cubicspline) The next one is the combined file, the jaxa/swiss/italy datasets (webp, cubicspline) The combined file was made with the was made with this process I agree it seems to be something related to the differing resolutions, becaue I have made a few of the datasets by themselfs and this effect is not as bad. |
|
Thanks @acalcutt , poking also @bertt and @frodrigo as they may also have been through that challenge and found a solution. One method is also described in this article : https://wwwtyro.net/2019/03/21/advanced-map-shading.html |
@utagawal this method uses map tiles as data source, those same kind of tiles we are trying to generate right now. |
|
I already see this issue. Not sure to remember what was the origin. at procedure-construction.sh.zip (using a rio-rgbify new optimisation https://medium.com/@frederic.rodrigo/optimization-of-rgb-dem-tiles-for-dynamic-hill-shading-with-mapbox-gl-or-maplibre-gl-55bef8eb3d86 ) |
|
@acalcutt I saw your posting on the GDAL forum, and unfortunatly no answer yet, so I gave a shot at asking ChatGPT, Claude and Gemini about this, here are their answers, maybe it would provide with some leads to an accepable solution ? Google Gemini's answer :
1. Reprojecting to a Common Reference System:
2. Resampling with Higher Order Interpolation:
3. Advanced Techniques:
Additional Considerations:
By trying these techniques and considering the additional points, you should be able to reduce or eliminate the lines in your final terrain RGB tiles. |
|
@acalcutt I gave a shot asking LLMs about the issue. I do not have the expertise to judge the validity of their answer but it may put you on the path for solving the challenge. I'll paste their answers so you can have a look. |
MistralAIThe issue you're encountering with the added lines when combining different source files is likely due to slight misalignments or differences in the elevation data between the datasets. Here are a few suggestions to help mitigate this issue: Use a different resampling method: Although you mentioned that you've tried different resampling methods, it might still be worth experimenting with other algorithms available in gdalwarp. For example, you can try using -r lanczos, -r average, or -r mode instead of -r cubic. Feathering: To create a smoother transition between datasets, you can apply feathering to the edges of the higher-resolution datasets (TINITALY and SWISSALTI). Feathering gradually reduces the weight of the dataset towards its edges, allowing for a more seamless blend with the underlying dataset (JAXA). You can use the -dstalpha option in gdalwarp to create an alpha band, and then manipulate it to create a feathering effect. Here's an example of how to create a feathering effect for the SWISSALTI dataset: Mosaicking: Instead of using VRTs, you can try creating a mosaic of the datasets using gdal_merge.py. This approach might help reduce the appearance of lines between datasets. However, keep in mind that this method might increase the processing time and memory requirements. Post-processing: After generating the terrain RGB tiles, you can apply post-processing techniques to reduce the visibility of the lines. For example, you can use a blur or smoothing filter to minimize the appearance of seams between datasets. By implementing these suggestions, you should be able to reduce or eliminate the added lines when combining different source files. Keep in mind that finding the optimal solution might require some trial and error, as well as adjusting the parameters to suit your specific use case. |
GeminiThe lines you're seeing in your final terrain RGB tiles likely arise from misalignment between the source DEMs. Here are some potential solutions to fix these lines: 1. Reprojecting to a Common Reference System:
2. Resampling with Higher Order Interpolation:
3. Advanced Techniques:
Additional Considerations:
By trying these techniques and considering the additional points, you should be able to reduce or eliminate the lines in your final terrain RGB tiles. |
ClaudeI'll help you analyze and solve the issue with the artifacts appearing when merging multiple DEM sources. Since you can't access the individual file links, could you paste the relevant code sections or error outputs directly in our conversation? That would help me diagnose the specific issues. From what you've described, it sounds like you're experiencing edge artifacts where different DEM sources meet. This is a common issue when merging elevation data from different sources. Here are some potential solutions:
# Add -blend parameter to gdalbuildvrt
gdalbuildvrt -blend 100 output.vrt file1.tif file2.tif file3.tif
gdalbuildvrt -resolution highest -r cubic -srcnodata -9999 -vrtnodata -9999 output.vrt filelist.txt
# Resample all inputs to the same resolution
gdalwarp -tr 0.0002777777777778 0.0002777777777778 -r cubic input.tif aligned_output.tif
# Specify NoData values explicitly
gdalbuildvrt -srcnodata -9999 -vrtnodata -9999 output.vrt filelist.txtLooking at your script, I can suggest a few modifications that might help eliminate those edge artifacts when merging the DEMs. Here's an enhanced version with explanations: #!/bin/bash
#custom version of rio rgbify which adds speed improvements is recommended https://github.com/acalcutt/rio-rgbify
#set max file limit
ulimit -s 65536
JAXA_INPUT_DIR=./jaxa
ITALY_INPUT_DIR=./tinitaly
SWISS_INPUT_DIR=./swissalti
OUTPUT_DIR=./output
jaxavrt1=${OUTPUT_DIR}/jaxa.vrt
jaxavrt2=${OUTPUT_DIR}/jaxa_warp.vrt
italy1=${OUTPUT_DIR}/italy.vrt
italy2=${OUTPUT_DIR}/italy_warp.vrt
swiss1=${OUTPUT_DIR}/swiss.vrt
swiss2=${OUTPUT_DIR}/swiss_warp.vrt
outvrt=${OUTPUT_DIR}/out.vrt
mbtiles=${OUTPUT_DIR}/jaxa_swiss_italy_terrainrgb_z0-Z12_webp.mbtiles
# Create output directory if it doesn't exist
[ -d "$OUTPUT_DIR" ] || mkdir -p $OUTPUT_DIR || { echo "error: $OUTPUT_DIR " 1>&2; exit 1; }
# Process JAXA data
gdalbuildvrt -overwrite \
-srcnodata -9999 -vrtnodata -9999 \
-r cubic \
-resolution highest \
${jaxavrt1} ${JAXA_INPUT_DIR}/*_DSM.tif
gdalwarp -r cubic \
-t_srs EPSG:3857 \
-dstnodata 0 \
-wo OPTIMIZE_SIZE=YES \
-wo NUM_THREADS=ALL_CPUS \
-multi \
-co COMPRESS=DEFLATE \
-co PREDICTOR=2 \
${jaxavrt1} ${jaxavrt2}
# Process Italy data
gdalbuildvrt -overwrite \
-srcnodata -9999 -vrtnodata -9999 \
-r cubic \
-resolution highest \
${italy1} ${ITALY_INPUT_DIR}/*.tif
gdalwarp -r cubic \
-t_srs EPSG:3857 \
-dstnodata 0 \
-wo OPTIMIZE_SIZE=YES \
-wo NUM_THREADS=ALL_CPUS \
-multi \
-co COMPRESS=DEFLATE \
-co PREDICTOR=2 \
${italy1} ${italy2}
# Process Swiss data
gdalbuildvrt -overwrite \
-srcnodata -9999 -vrtnodata -9999 \
-r cubic \
-resolution highest \
${swiss1} ${SWISS_INPUT_DIR}/*.tif
gdalwarp -r cubic \
-t_srs EPSG:3857 \
-dstnodata 0 \
-wo OPTIMIZE_SIZE=YES \
-wo NUM_THREADS=ALL_CPUS \
-multi \
-co COMPRESS=DEFLATE \
-co PREDICTOR=2 \
${swiss1} ${swiss2}
# Create ordered list of VRTs
rm -f filenames.txt
printf '%s\n' ${jaxavrt2} >filenames.txt
printf '%s\n' ${italy2} >>filenames.txt
printf '%s\n' ${swiss2} >>filenames.txt
# Create final VRT with blending
gdalbuildvrt -overwrite \
-resolution highest \
-r cubic \
-srcnodata 0 -vrtnodata 0 \
-blend 100 \
${outvrt} -input_file_list filenames.txt
# Create terrain mbtiles
rio rgbify -b -10000 -i 0.1 \
--min-z 0 --max-z 12 \
-j 24 --format webp \
${outvrt} ${mbtiles}
# Update mbtiles metadata
sqlite3 ${mbtiles} '
CREATE UNIQUE INDEX tile_index on tiles (zoom_level, tile_column, tile_row);
UPDATE metadata SET value = "jaxa_swiss_italy_terrainrgb_z0-Z12_webp" WHERE name = "name" AND value = "";
UPDATE metadata SET value = "JAXA AW3D30 (2024), swissALTI3D (2024), and Tinitaly (2023) converted with rio rgbify" WHERE name = "description";
UPDATE metadata SET value = "webp" WHERE name = "format";
UPDATE metadata SET value = "1" WHERE name = "version";
UPDATE metadata SET value = "baselayer" WHERE name = "type";
INSERT OR REPLACE INTO metadata (name,value) VALUES
("attribution","<a href=https://earth.jaxa.jp/en/data/policy/>JAXA AW3D30 (2024)</a> | <a href=https://www.swisstopo.admin.ch/en/height-model-swissalti3d>swissALTI3D (2024)</a> | <a href=https://tinitaly.pi.ingv.it/>Tinitaly (2023)</a>"),
("minzoom","0"),
("maxzoom","12"),
("bounds","-180,-90,180,90"),
("center","0,0,5");
'
Key changes I've made to help address the edge artifacts:
Additional suggestions if you're still seeing artifacts:
gdalwarp -r cubic -t_srs EPSG:3857 -dstnodata 0 -wo CUTLINE_ALL_TOUCHED=TRUE -wo CUTLINE_BLEND_DIST=100 ...
|
ChatGPTIt sounds like you're running into issues with "seams" or visible lines appearing where the datasets overlap, which can happen when DEM layers are combined in a VRT (Virtual Raster) with different resolutions or varying data values along the edges. Here are a few strategies to address this: 1. Adjust Blending or Merging StrategyWhen combining multiple DEM sources, blending is key to avoiding these lines. Try using the
Updated command example: gdalbuildvrt -resolution average -addalpha -srcnodata -9999 -vrtnodata -9999 ${jaxavrt1} ${JAXA_INPUT_DIR}/*_DSM.tif2. Use Smoothing Options with
|
|
I made a new jaxa/swiss/italay test and tried to take into account the AI suggestions above 1.) In any command I saw there was a place to add resampling algorithm (-r), I set it to lanczos. So basically all gdalbuildvrt, gdalwarp, and gdal_translate commands all had that option and '-r lanczos' has been added. 2.) In gdalbuildvrt commands I added '-resolution highest' Just with those changes, it seems to have fixed the lines I was seeing on the merged terrain. I will note I still see the band artifacts we were looking into in #5 but that is a different issue. I am now trying to increase the bounds to cover at least tiles 5/16/11, 5/17/11, 5/16/12, 5/17/12, so it works better at blending into the full planet terrain file I made. I was testing layering them with addprotocol like I am also testing france using RGE ALTI 1m tiles, but that has been the hardest dataset to work with so far. I was thinking about switching to the datacet the 'procedure-construction.sh.zip' above uses, but that looks like it is 25m Just on a side note, the AI either made up a few answers or they are outdated. |
|
You're getting close ! This is very good news. I challenged Claude about de
|
1 similar comment
|
You're getting close ! This is very good news. I challenged Claude about de
|
|
And about
I believe it is not worht using the 1m LIDAR version. The 5m version should be good enough and will produce less artefacts : https://geoservices.ign.fr/rgealti |
At 45° of latitude, 25m/pixel is equivalent tiles 11 resolution. Unless you want tiles 12 or more, at this latitude, 1m/pixel or 5m/pixel is not required. |
|
Claude released a new model, so I reasked the question and here is its answer : The issue seems to be that when you merge different DEM sources in the final VRT, you're getting visible artifacts/lines at the boundaries where they meet. The fact that these lines don't appear when using individual sources (like your JAXA-only example) suggests this is specifically a data integration issue. Key observations:
Instead of my previous overly complex suggestions, I'd focus on:
gdalwarp -r cubic -t_srs EPSG:3857 -dstnodata 0 -tap -tr 30 30 ${sourcevrt} ${outputvrt}The
Would you be able to share what the artifacts look like in more detail? Are they straight lines that follow tile boundaries, or do they follow the edges of where different datasets meet? |
|
Have you looked into the several vertical references of the source data? And what's the target vertical reference system of the resulting model? |
|
BTW, just saw https://app.reliefmaps.io/ and they nailed it. The creator is from France :) |
|
I have been working on a tool that lets you merge and layer terrrain raster-dems, which I think will make some similar results. Initially it started as a standalone rasterio script, but now I am in the process of trying to build it into rio rgbify in this branch The main code is in basically you supply it with a config file like this. if you supply bounds it gets all the tiles in the range and requests them from each source. if you don't supply bounds it uses the bounds of the last file in the list. if one source doesn't have the tile at that zoom level for the source it tries to upscale a lower tile. in the masks out certain elevations (0 by default) for all sources but the first one, then merges them together In my test it seemed to work pretty well. I used it to merge my jaxa and swisstopo and the results wern't bad for example, swiss topo looked like this and i used the tool to merge it with JAXA and it looks like this My idea was this could be used to make sparse datasets. So I could make high resultion datasets of all the sources we have, then use this tool to merge the datasets together later. this would avoid having to get all the sources in the same format and trying to process them all together. Just to note, the tool is a work in process and I don't have the main rgbify function working right now (I am working on trying to improve the memory use of that right now) |
|
This is it !! |
|
Just to note that when you do go beyond the zoom of the original raster-dem, it still needs to upscale it and that leads to a bit of blockyness. I am sure this can be somewhat improved by playing with the resampling used (like we talked about above using highest and using biliner or lanczos in all steps. including when creating the original raster-dem) For example
however when you are over zooming something I think that is expected. at least this allows you to minimize the effect since it only happens beyond the original zoom level and you can confine it to the are of the higer resolution file by making a sparse dataset. This is a more zoomed in example. i will upload it soon so it can be looked at |
|
@acalcutt thanks a lot for the amazing work you're doing. In the more zoomed example, it may be worth to try to avoid "banding" using a surface blur filter. |
|
On the topic of More Detailed DEMS, I have uploaded a torrent of both SwissAlti and TIN Italy. I am also trying to generate RGE Alti 5m right now, but the way the one is formatted it seems to take a bit longer I also added demo of the test merged terrain at SwissAlti, which merged with jaxa. you can see a demo of that here. Also, in the https://github.com/acalcutt/rio-rgbify/tree/merge I tried adding a gaussian_blur_sigma setting which gets multiplied by the zoom level difference when overzooming. It does seem to slightly work. my example above used a 1.5 gaussian_blur_sigma setting. I also started to re-organize some of these scripts to separate each dataset at https://github.com/acalcutt/terrain_merged/tree/reorg , but some of them are still untesed to I haven't merged it into the main branch yet. |
|
Getting better and better... @acalcutt can you plase add a map layer on the merged example with location names so we can identify where the border is (and spot potential artefacts) ? Here ius the datasources used by Sonny Dems to generate the European DEM : https://drive.google.com/drive/folders/0BxphPoRgwhnoWkRoTFhMbTM3RDA?resourcekey=0-wRe5bWl96pwvQ9tAfI9cQg It could be interesting in your script to have a config file to specify the list of DEM sources ? |
|
How about something like this for a demo, I added the new merged raster-dem as terrain and hillshade, but added the just swisstopo as a raster layer, so you can see the edges I also added this one which uses the demptiles, but the border isn't very accurate in that Edit: here is one with osm raster tiles added in also |
My merge script already uses a config file, so if we make the terrain raster-dem file from the sonny dems, we could merge and layer it in with the other raster-dems. I don't think there is a limit, though i have only tested up to 3 sources so far |
|
Spotted :
Here : https://web-platform-emyprgwe.stackblitz.io/#17.22/45.939009/7.036383/10/33 Maybe you need to add yet another filter to detect -at merging border- artefact when there is a "big difference" (=stairs step) in elevation and when it occurs reapply a smoothing algorithm ? What about looking at antialiasing methods (only here the pixel is bigger...) ? https://www.reddit.com/r/pcmasterrace/comments/4r8xze/quick_cheatsheet_for_antialiasing_types/ |
Hi All, I am very late to the game on this and not even sure how I ended up on this thread. I will have a look an see if there is anything I could add. |
|
I guess 'lanczos' is not the resampling setting to go with. I tried to make the SONNY DEM of europe with 'lanczos' and was really unhappy with how it came out, it had lots of that banding we talked about before I then re-ran it in 'bilinier' and it actually looks good with that setting We had talked about using bilinier before in #5 (comment) , but I had forgotten |
|
I did a 3rd test with 'cubic' resampling and i think that one looks the best, especially when zooming in a bit. this was the other one we discussed trying in #5 (comment) . I think i will try and use this going forward and see how it goes. |
|
Definitley the best rendering is with cubic. Did you find a way forward to fix the DEM boundaries differences in height/sampling definition ? |
|
I haven't looked at it, I am trying to get some source terrain files created so i can test the merging with them right now Really it is going to do something like that anywhere there is a big height difference between sources. I do think once we merge in some of the other higher resolution sources in those areas it will not be quite so bad I'm guessing if you had to apply anti-aliasing it would have to be over the final image, which might lower the quality of your higher res sources. |
|
According to the tests we made, the merge should happen on the sampled data before the rendering happens, else this will create some artefacts. |
Any possibiliyt to allow other inputs and make the config files more customizable ?
For example being able to specify DEMs with higher resolution (we have 5m now for most western countries and even LIDAR) and as a fallback use AW3D30 ?
This will certainly require to specify URL, file format and re-projection.
Challenge will be to be carefull on how to manage potential artefacts at tiles and country borders when multiple DEMS overlap and do not have the same resolution.
The text was updated successfully, but these errors were encountered: