More detailed DEMS

[utagawal]

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 ?

Dataset	Viewport	Disk space	Resolution	Format	Comments	Link
ETOPO1	Globale	0.45 GB	2000m	GeoTiff	Contient la bathymétrie	https://github.com/DahnJ/Awesome-DEM#etopo1
GEBCO_2019	Globale	10 GB	200m	GeoTiff	Contient la bathymétrie	https://github.com/DahnJ/Awesome-DEM#gebco
NASADEM	Globale	120 GB	30m	GeoTiff ou HGT		https://portal.opentopography.org/datasetMetadata?otCollectionID=OT.112016.4326.2
IGN RGE Alti	France Dom/Tom	600 GB	5m et 1m	ASC	Certaines zones basse résolution (voir fichiers liés dans livraison)	https://geoservices.ign.fr/telechargement?bdd%5B0%5D=1060&bdd%5B1%5D=1060
TIN Italy	Italie	15 GB	10m	GeoTiff		https://tinitaly.pi.ingv.it/
Sonny DEM	Europe	?	20m	?	Composite de très bonne qualité	https://sonny.4lima.de/
SwissAlti	Suisse	700 GB	2m et 50cm	GeoTiff		STAC (https://data.geo.admin.ch/api/stac/v0.9/) ou https://www.swisstopo.admin.ch/fr/geodata/height/alti3d.html#musterdaten

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.

[acalcutt]

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. With the jaxa dataset, zoom 11 looks pretty much as good as it will get, and that takes 138GB. just adding 1 level to zoom 12 brings that up to 427GB. I'd imagine going one more zoom level up to 13 is going to bring that file to at least 1TB, possibly 2TB

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.

[utagawal]

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 :

• No need to go more detailed than Z15
• projection 3857 (Web Mercator) must be used for all datasets (hence maybe some pre-processing whit GDAL potentially required for some)
• Tiles should have a size of 512px for better compression rates
• Hillshading needs to be precalculated with raster tiles, it is easier to compress see : https://medium.com/@frederic.rodrigo/optimization-of-rgb-dem-tiles-for-dynamic-hill-shading-with-mapbox-gl-or-maplibre-gl-55bef8eb3d86
• GDAL merge could be used to take care of the overlapping DATAsets

Zoom level	Resolution (meters/px)	Best dataset	Dataset native resolution (m/px)	Coverage	Dataset SRID	Alternate datasets (*)
1	55,346	ETOPO1	1800	Global + Bathymetry	4326	N/A
2	27,673	ETOPO1	1800	Global + Bathymetry	4326	N/A
0	13,837	ETOPO1	1800	Global + Bathymetry	4326	N/A
3	6,918	ETOPO1	1800	Global + Bathymetry	4326	N/A
4	3,459	ETOPO1	1800	Global + Bathymetry	4326	N/A
5	1,730	GEBCO_2019	464	Global + Bathymetry	4326	N/A
6	865	GEBCO_2019	464	Global + Bathymetry	4326	N/A
7	432	GEBCO_2019	464	Global + Bathymetry	4326	N/A
8	216	NASADEM	30	Global	4326	N/A
9	108	NASADEM	30	Global	4326	N/A
10	54	NASADEM	30	Global	4326	N/A
11	27	IGN_5m	5	France	2154	SwissAlti2m, TINItaly, IGN_spain
12	14	IGN_5m	5	France	2154	SwissAlti2m, TINItaly, IGN_spain
13	7	IGN_5m	5	France	2154	SwissAlti2m, TINItaly, IGN_spain
14	3.4	IGN_1m	1	France	2154	SwissAlti2m, TINItaly, IGN_spain
15	1.7	IGN_1m	1	France	2154	SwissAlti50cm, TINItaly, IGN_spain
16	0.8	IGN_1m	1	France	2154	SwissAlti50cm, TINItaly, IGN_spain

[utagawal]

Maybe some ideas here : https://github.com/tilezen/joerd/blob/master/docs/data-sources.md

[utagawal]

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. With the jaxa dataset, zoom 11 looks pretty much as good as it will get, and that takes 138GB. just adding 1 level to zoom 12 brings that up to 427GB. I'd imagine going one more zoom level up to 13 is going to bring that file to at least 1TB, possibly 2TB

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.

Now with the Webp compression I guess it opens new opportunity size wise ?

Also to continue the discussion :

• When we have a dataset available till Z15 for a country => get the data and calculate the tiles
• When you don't, keep the dataset from let say Z12 and generate the tiles based on this. When compression happens, it should deduplicate that and not generate additional storage requirement.

Another approach would be to keep the terrain data at Z12 and work on a hillshading mbtile till Z15 ?

[acalcutt]

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 https://tiles.wifidb.net/data/swissalti3d_terrarium_0-17/#8/46.695/8.064 https://stackblitz.com/edit/web-platform-9f92n2?file=index.html

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 https://github.com/maplibre/maplibre-style-spec/issues/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

[acalcutt]

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.

#!/bin/bash

#Requires custom version of rio rgbify which adds terrarium encoding support ( https://github.com/acalcutt/rio-rgbify/ )

INPUT_DIR=./download
OUTPUT_DIR=./output
vrtfile=${OUTPUT_DIR}/swiss_italy_terrainrgb0-15.vrt
vrtfile2=${OUTPUT_DIR}/swiss_italy_terrainrgb0-15_warp.vrt
mbtiles=${OUTPUT_DIR}/swiss_italy_terrainrgb0-15.mbtiles
#vrtfile2=${OUTPUT_DIR}/swiss_italy_terrainrgb0-15_warp.vrt

[ -d "$OUTPUT_DIR" ] || mkdir -p $OUTPUT_DIR || { echo "error: $OUTPUT_DIR " 1>&2; exit 1; }

#rm rio/*

for i in /mnt/usb/tinitaly/input/*.tif;
  do gdalwarp -of VRT $i ${i%%.*}_out.vrt -t_srs "EPSG:3857";
done;

for i in /mnt/usb/swisstopo_to_mbtiles/download/*.tif;
  do gdalwarp -of VRT $i ${i%%.*}_out.vrt -t_srs "EPSG:3857";
done;

printf '%s\n' /mnt/usb/tinitaly/input/*.vrt >filenames.txt
printf '%s\n' /mnt/usb/swisstopo_to_mbtiles/download/*.vrt >>filenames.txt

gdalbuildvrt -overwrite ${vrtfile} -input_file_list filenames.txt
gdalwarp -r cubicspline -t_srs EPSG:3857 -dstnodata 0 -co COMPRESS=DEFLATE ${vrtfile} ${vrtfile2}
rio rgbify -e mapbox -b -10000 -i 0.1 --min-z 0 --max-z 15 -j 16 --format png ${vrtfile2} ${mbtiles}

[utagawal]

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.

[utagawal]

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 Will you be able to make the mbtile available ?

[utagawal]

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. An option to spped up the discovery of these overlopping tiles would be to focus at the country borders (since usually detailed datasets are per country). Border lines can be found here : https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-boundary-lines/

[acalcutt]

Just an FYI, I did make a new JAXA 2024 webp file https://techidiots.net/downloads/open-base-map-torrents/mbtiles/jaxa_terrainrgb_z0-12_webp-torrent/view I am still working on other formats, but that one should e ready at least

I was also made a new swissalti file https://techidiots.net/downloads/open-base-map-torrents/mbtiles/swissalti3d_terrainrgb_z0-z16-torrent/view

With generating these I was thinking above the original topic of creating a merged file. I am testing something like this right now https://github.com/acalcutt/Terrain_Merged/ with that I am trying to merge jaxa, swissalti, and tinitally as vrts merged into one combined vrt. Do you think something like that should work? Right now I have it generating, and the combined vrt looks like it should work and had all the dataset it it

[utagawal]

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).

[utagawal]

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 ? https://github.com/acalcutt/jaxa_AW3D30_to_MBTiles/issues/7 ?

[acalcutt]

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 https://stackblitz.com/edit/web-platform-94jnfn?file=index.html

Compared to the new jaxa webp you just downloaded it made in cubicspline https://stackblitz.com/edit/web-platform-l2g2e5?file=index.html

The issue i saw looks like this in the jaxa area

[utagawal]

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 ?

[xfischer]

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 ? @nbozon may surely have the good persons to talk with (dreaming is not forbidden), or even @Zemledelec (sorry for poking) who has done an amazing job with OpenGlobus

[Zemledelec]

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!

[acalcutt]

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) https://stackblitz.com/edit/web-platform-al1zf5?file=index.html you can see the example looks pretty smooth (other than some small lines we have noticed in the past in)

The next one is the combined file, the jaxa/swiss/italy datasets (webp, cubicspline) https://stackblitz.com/edit/web-platform-wrcpum?file=index.html you can see the improved detail in the swiss(blue) and italy(yellow), but it introduced weird lines all over it seems. Even in France(orange), which should be using the jaxa dataset above has these lines.

The combined file was made with the was made with this process https://github.com/acalcutt/terrain_merged/blob/main/create_terrainrgb.sh

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.

[utagawal]

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

[xfischer]

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.

[frodrigo]

I already see this issue. Not sure to remember what was the origin. But there is a script to produce this data for France here

https://www.data.gouv.fr/fr/datasets/bd-alti-r-25-m-tuiles-pour-courbes-de-niveau-et-ombrage-dynamiques-1/

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 )

[utagawal]

@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 :

The 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:

• Currently, you're warping each VRT to EPSG:3857 (WGS84 Pseudo-Mercator). Ensure all the source DEMs have the same spatial reference system (SRS) before creating the VRTs. This eliminates any potential misalignment due to different projections. Check the metadata of the source files to determine their SRS and reproject them if necessary using gdalwarp before creating VRTs.

2. Resampling with Higher Order Interpolation:

• While you've tried different resampling methods, consider using a higher order interpolation like cubic spline or Lanczos resampling during warping with gdalwarp. These methods can better handle slight misalignments compared to nearest neighbor resampling.

3. Advanced Techniques:

• Feathering: Investigate techniques like feathering, where the edges of each DEM are blended with the underlying DEM to create a smoother transition. This can be achieved using specific tools depending on your programming language.
• Seamless Mosaicking: Explore tools designed for seamless mosaicking, which can account for small geometric offsets and radiometric differences between datasets. Libraries like Rasterio offer functionalities for this.

Additional Considerations:

• Data Quality: The lines might also be inherent to the data itself. Lower quality DEMs like JAXA might have inherent geometric or radiometric inconsistencies compared to the higher quality SwissALTI. This could lead to artifacts even after alignment.
• Data Preprocessing: Consider any preprocessing steps specific to your data formats. For instance, the script mentions converting JAXA data to Float32 before VRT creation. Explore any other preprocessing steps recommended for your specific DEM sources.

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.

[utagawal]

@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.

[utagawal]

MistralAI

The 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:

# Create an alpha band for the Swiss VRT
gdalwarp -r cubic -t_srs EPSG:3857 -dstnodata 0 -dstalpha ${swiss1} ${swiss2}_alpha.vrt

# Create a feathering effect (5km buffer in this example)
gdal_calc.py -A ${swiss2}_alpha.vrt --outfile=${swiss2}_feather.vrt --calc="A*(1-(distance(A)/5000))" --NoDataValue=0

# Warp the Swiss VRT with feathering applied
gdalwarp -r cubic -t_srs EPSG:3857 -dstnodata 0 ${swiss2}_feather.vrt ${swiss2}
Ordering of datasets: Ensure that the datasets are ordered correctly in the filenames.txt file, with the lowest priority (JAXA) at the top and the highest priority (SWISSALTI) at the bottom.

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.

[utagawal]

Gemini

The 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:

• Currently, you're warping each VRT to EPSG:3857 (WGS84 Pseudo-Mercator). Ensure all the source DEMs have the same spatial reference system (SRS) before creating the VRTs. This eliminates any potential misalignment due to different projections. Check the metadata of the source files to determine their SRS and reproject them if necessary using gdalwarp before creating VRTs.

2. Resampling with Higher Order Interpolation:

• While you've tried different resampling methods, consider using a higher order interpolation like cubic spline or Lanczos resampling during warping with gdalwarp. These methods can better handle slight misalignments compared to nearest neighbor resampling.

3. Advanced Techniques:

• Feathering: Investigate techniques like feathering, where the edges of each DEM are blended with the underlying DEM to create a smoother transition. This can be achieved using specific tools depending on your programming language.
• Seamless Mosaicking: Explore tools designed for seamless mosaicking, which can account for small geometric offsets and radiometric differences between datasets. Libraries like Rasterio offer functionalities for this.

Additional Considerations:

• Data Quality: The lines might also be inherent to the data itself. Lower quality DEMs like JAXA might have inherent geometric or radiometric inconsistencies compared to the higher quality SwissALTI. This could lead to artifacts even after alignment.
• Data Preprocessing: Consider any preprocessing steps specific to your data formats. For instance, the script mentions converting JAXA data to Float32 before VRT creation. Explore any other preprocessing steps recommended for your specific DEM sources.

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.

[utagawal]

Claude

I'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:

1. Try adding a feathering/blending zone:

   # Add -blend parameter to gdalbuildvrt
   gdalbuildvrt -blend 100 output.vrt file1.tif file2.tif file3.tif

2. Consider adding overlap handling:

   gdalbuildvrt -resolution highest -r cubic -srcnodata -9999 -vrtnodata -9999 output.vrt filelist.txt

3. You might want to pre-process the DEMs to ensure consistent resolution and alignment:

   # Resample all inputs to the same resolution
   gdalwarp -tr 0.0002777777777778 0.0002777777777778 -r cubic input.tif aligned_output.tif

4. Check for and handle NoData values consistently:

   # Specify NoData values explicitly
   gdalbuildvrt -srcnodata -9999 -vrtnodata -9999 output.vrt filelist.txt

Looking 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:

1. Added -resolution highest to all gdalbuildvrt commands to ensure consistent resolution handling

2. Added -blend 100 to the final VRT creation to smooth transitions between datasets

3. Standardized NoData handling across all steps

4. Added optimization parameters to gdalwarp commands:

   • -wo OPTIMIZE_SIZE=YES
   • -wo NUM_THREADS=ALL_CPUS
   • -multi
   • Added compression options

5. Made the NoData values consistent across all operations (using 0 and -9999 appropriately)

Additional suggestions if you're still seeing artifacts:

1. You might want to try different blending values (increase or decrease the 100 in -blend 100)

2. Consider adding a buffer zone around each dataset:

   gdalwarp -r cubic -t_srs EPSG:3857 -dstnodata 0 -wo CUTLINE_ALL_TOUCHED=TRUE -wo CUTLINE_BLEND_DIST=100 ...

3. If you're seeing specific artifacts at certain zoom levels, you might want to adjust the rio rgbify parameters or process different zoom levels separately with different parameters.

[utagawal]

ChatGPT

It 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 Strategy

When combining multiple DEM sources, blending is key to avoiding these lines. Try using the -resolution and -addalpha flags with gdalbuildvrt and gdalwarp. Here’s how:

• Use -resolution average in gdalbuildvrt commands to ensure GDAL picks the optimal resolution automatically by averaging across sources.

• Enable -addalpha in gdalwarp to include an alpha band, which can improve transparency handling across overlapping datasets.

  Updated command example:

  gdalbuildvrt -resolution average -addalpha -srcnodata -9999 -vrtnodata -9999 ${jaxavrt1} ${JAXA_INPUT_DIR}/*_DSM.tif

2. Use Smoothing Options with gdalwarp

With gdalwarp, consider -r bilinear or -r lanczos for smoother edges instead of -r cubic, which can sometimes create hard lines at boundaries.

   gdalwarp -r bilinear -t_srs EPSG:3857 -dstnodata 0 ${jaxavrt1} ${jaxavrt2}

3. Check Alignment and Nodata Value Consistency

Ensure all three datasets have the same nodata values set before you combine them. You can standardize nodata values using gdal_translate:

   gdal_translate -a_nodata 0 -of GTiff input.tif output_standardized.tif

4. Consider Using gdal_dem_mosaic (Optional Tool)

The gdal_dem_mosaic.py script (if available) is specifically designed for mosaicking elevation data and can help manage issues with seams by blending DEMs where overlap occurs.

5. Try a Custom Python Script for Smoother Results (Advanced)

If blending still poses issues, you could consider creating a custom Python script using GDAL or Rasterio, manually interpolating or smoothing the overlap areas.

Testing these adjustments should help to reduce or eliminate the lines, providing a more seamless VRT output for your tiles.