der-stefan / OpenTopoMap

A topographic map from OpenStreetMap and SRTM data
https://opentopomap.org
Other
455 stars 118 forks source link

holes in the map at the hill shadows #337

Open a14stoner opened 1 year ago

a14stoner commented 1 year ago

Hello, I made the import the exact way as described here: https://github.com/der-stefan/OpenTopoMap/blob/master/mapnik/README.md

When rendering the map some holes appear at the shadow side of the hills - you can see it in the attached image.

image

image

i see the holes come from the hillshade tifs

image

I tried to recreate the hillshade tifs but no sucess. the holes are always there.

sletuffe commented 1 year ago

Have you forgotten the void filling step ?

Fill all voids

cd unpacked
for hgtfile in *.hgt;do gdal_fillnodata.py $hgtfile $hgtfile.tif; done

https://gdal.org/programs/gdal_fillnodata.html

Anyway, even if you did, the result is still varying in quality from places to places. As said here : https://github.com/der-stefan/OpenTopoMap/issues/334 there still is a lack of a world wide, open and high quality DEM.

AFAIK, both OTM and my project https://github.com/RefugesInfo/mri are using conflated tiff DEM made and assembled from https://www.opensnowmap.org

a14stoner commented 1 year ago

Hey @sletuffe yea i filled the voids. but didnt help. Is it possible that i get the conflated tiff DEM which you and OTM are using?

Unfortunately the NASA SRTM files are not possible to download: http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/N44E006.hgt.zip that link does not work anymore.

sletuffe commented 1 year ago

Can't find it anymore, looks like I deleted the raw one, I only have left the contour lines and hillshaded (multi zoom) versions

@yvecai is it ok for you that I give @a14stoner the link you sent me 3 years ago from the raw dem on your servers ? (it still works, but, well, that's a ~400GB dowload, so better ask)

a14stoner commented 1 year ago

Alright. This would be nice. Maybe its also possible to snip out a BBOX of [6,44,18,56] (EPSG:4326). Thats the area i need the data.

yvecai commented 1 year ago

I'll double check my current bandwidth later this evening and will send the download URL to @a14stoner

yvecai commented 1 year ago

That's OK, @a14stoner, you'll find my email address at the bottom of the 'About' page at opensnowmap.org, I can't put these urls on a public discussion.

a14stoner commented 1 year ago

Ill try my best with the data from @yvecai and report back if I managed it ;)

a14stoner commented 1 year ago

I identified which step produced the voids:

gdaldem hillshade -z 2 -co compress=lzw -co predictor=2 -co bigtiff=yes -compute_edges warp-90.tif hillshade-90.tif && gdal_translate -co compress=JPEG -co bigtiff=yes -co tiled=yes hillshade-90.tif hillshade-90-jpeg.tif

Before the second command it looks like this:

image

after it:

image

so i will continue with the not compressed file and check it its better now.

P.S.: when compressing with lzw at the last step no holes are there PP.S: rendering it with mapnik results in a nice map without holes!!:

image

Just an information for @yvecai and @sletuffe : I got ASTER DEMs from 2019 here: https://e4ftl01.cr.usgs.gov/ASTT/ASTGTM.003/2000.03.01/

sletuffe commented 1 year ago

It's been a long time since I've played with gdal tools, but I'm quitte surprised the jpeg compression step could add those holes. It's like if those holes where allready there but unseen (like nodata values beeing rendered as black in the shadows of moutains) but the jpeg compression made them "appear" as white.

Is is possible that you are missing the -compute_edges switch ? If [-compute_edges](https://gdal.org/programs/gdaldem.html#cmdoption-compute_edges) is specified, gdaldem will compute values at image edges or if a nodata value is found in the 3x3 window, by interpolating missing values.

The thing is that without jpeg and tiled=yes your hillshaded tiff will be much bigger and also more time consuming to extract tiles from.

Anyway, does the ASTER DEMs from 2019 are better in areas like https://opentopomap.org/#marker=13/60.77559/25.65411 ? I'll be interested to know if those areas have been improved ?

https://github.com/der-stefan/OpenTopoMap/issues/327

yvecai commented 1 year ago

Maybe the artifacts comes from the fact that the hillshade is fully black, but the hillshade information is in the alpha channel. I don't know how the JPEG compression works with an alpha channel.

Apparently the holes from #327 are not present in the Opensnowmap's DEM, these area being covered by EU-DEM, though.

What's new(ish) in the open DEM family is this one : https://www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/aw3d30_e.htm Haven't used it, but some hillshade extract can be obtained from here to have a look: https://portal.opentopography.org/raster?opentopoID=OTALOS.112016.4326.2

sletuffe commented 1 year ago

Apparently the holes from #327 are not present in the Opensnowmap's DEM, these area being covered by EU-DEM, though.

It seams you are better than us at sweeping it under the carpet by painting it white ;-) Here is a more visible place in Norway on opensnowmap : https://www.opensnowmap.org/#map=13/11.836/63.062&b=snowmap&m=false&h=false

This coloring by altitude + holes in white gives a better idea of missing parts in the DEM : https://maps.refuges.info/?zoom=5&lat=65.63407&lon=41.15005&layers=B0

That being said, and I apologize in advance to my beloved Finns, Norwegians or Russians (that's where the holes are), but the inconvenience is rather limited to mostly non-hiking areas. The might to solve this issue is therefore rather low...

yvecai commented 1 year ago

Ah yes, most certainly the Aster parts. Needs to find some disk for Opensnowmap's server this summer to try and fix this with ALOS or else, I hope it rains a lot.

Le 31 janvier 2023 18:19:03 GMT+01:00, sly - sylvain letuffe @.***> a écrit :

Apparently the holes from #327 are not present in the Opensnowmap's DEM, these area being covered by EU-DEM, though.

It seams you are better than us at sweeping it under the carpet by painting it white ;-) Here is a more visible place in Norway on opensnowmap : https://www.opensnowmap.org/#map=13/11.836/63.062&b=snowmap&m=false&h=false

This coloring by altitude + holes in white gives a better idea of missing parts in the DEM : https://maps.refuges.info/?zoom=5&lat=65.63407&lon=41.15005&layers=B0

That being said, and I apologize in advance to my beloved Finns, Norwegians or Russians (that's where the holes are), but the inconvenience is rather limited to mostly non-hiking areas. The might to solve this issue is therefore rather low...

-- Reply to this email directly or view it on GitHub: https://github.com/der-stefan/OpenTopoMap/issues/337#issuecomment-1410772084 You are receiving this because you were mentioned.

Message ID: @.***>

sletuffe commented 1 year ago

What's new(ish) in the open DEM family is this one : https://www.eorc.jaxa.jp/ALOS/en/dataset/aw3d30/aw3d30_e.htm Haven't used it, but some hillshade extract can be obtained from here to have a look: https://portal.opentopography.org/raster?opentopoID=OTALOS.112016.4326.2

I was curious and got a shot at it. I'm not interested in filling the nothern holes but I'm a bit disappointed with the current DEM when we look at contour lines where there are vertical cliffs. Near my home town, there is such a very high cliff and that's where I look at to evaluate the relative positionning of summits and cliff.

The result (N45 E6) is just as it was, no progress, the Dent de Crolles summit is still hanging in the middle of the cliff :

Sans titre

(This is the tiff tile + contour generation + OSM extract of ftp://ftp.eorc.jaxa.jp/pub/ALOS/ext1/AW3D30/release_v2012_single_format/N045E005/N045E005.zip )

yvecai commented 1 year ago

At Opensnowmap I added a layer with the French hi-res Dem where you can clearly see the cliff at the Dent de Crolles. Is it a simple issue of resolution (30m)? A flaw in gdal contour algorithm?

a14stoner commented 1 year ago

Thank you for this nice convo - i lean things from you :D

With the not compressed tif my map looks like this now:

image

All is good except the contour lines are missing. I imported everything according to the tutorial.

took warp-90.tif > converted it to a pbf > imported it with osm2pgsql

No errors when loading it but just no lines. Here is an example of one database entry:

image

Its the first time i make such a map so sorry if my questions are weird :D

a14stoner commented 1 year ago

I dont know if just skip the compress part to JPEG is good enough to close this issue? Can it be that stefan used another source for his DEMs where compressing with jpeg also works?

for my other issue i created #338

sletuffe commented 1 year ago

Yes, sorry, I guess you are right : If you are not making any mistake, that means the Readme needs improvement or updates

Maybe let it open

a14stoner commented 1 year ago

since #338 is done now ill focus again on this one because rendering is slow with the LZW compressed tiffs.

I googled a little bit and and came to no result. But is has to do something with how transparency is handled with JPEG compressed TIFF files.

Interesting thread here https://gis.stackexchange.com/questions/114370/compression-artifacts-and-gdal

"The truth is that transparency can't be handled with transparent pixels with JPEG compression not other lossy compression methods."

So I would need a basis with no "nodata" parts. or I just let i be as it is and the seeding of the whole map takes just more time.

a14stoner commented 1 year ago

I found something out. When editing the warp file with gdal_edit.py -a_nodata 0 warp-1000.tif we can really see the spots where the tif is transparent: Before:

image

After:

image

The holes are the same. there is no issue with transparency?

What i can also see it that there is no alpha channel in the file before and after the gdaldem:

Before:

Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
  PREDICTOR=2
Corner Coordinates:
Upper Left  (  556581.993, 7361893.069) (  4d59'59.50"E, 55d 0' 0.50"N)
Lower Left  (  556581.993, 5311893.069) (  4d59'59.50"E, 42d59'58.14"N)
Upper Right ( 1892581.993, 7361893.069) ( 17d 0' 4.87"E, 55d 0' 0.50"N)
Lower Right ( 1892581.993, 5311893.069) ( 17d 0' 4.87"E, 42d59'58.14"N)
Center      ( 1224581.993, 6336893.069) ( 11d 0' 2.19"E, 49d21'51.27"N)
Band 1 Block=256x256 Type=Int16, ColorInterp=Gray

After:

Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=BAND
  JPEGTABLESMODE=1
  JPEG_QUALITY=75
Corner Coordinates:
Upper Left  (  556581.993, 7361893.069) (  4d59'59.50"E, 55d 0' 0.50"N)
Lower Left  (  556581.993, 5311893.069) (  4d59'59.50"E, 42d59'58.14"N)
Upper Right ( 1892581.993, 7361893.069) ( 17d 0' 4.87"E, 55d 0' 0.50"N)
Lower Right ( 1892581.993, 5311893.069) ( 17d 0' 4.87"E, 42d59'58.14"N)
Center      ( 1224581.993, 6336893.069) ( 11d 0' 2.19"E, 49d21'51.27"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Gray
  NoData Value=0

But only the NoData Value is new which is 0

yvecai commented 1 year ago

"rendering is slow with the LZW compressed tiffs"

You seem to focus on rendering speed, but I doubt you will ever be satisfied with lossy Jpeg compression : even if you fix these holes artifacts, you'll probably see others coming from the jpeg compression. Maybe you need to check other things: ssd or nvme disk, more downsampled tiff, one for each zoom layer, geotiff block size, etc...

a14stoner commented 1 year ago

talking about downsampled tiffs for each zoom layer: Right now at zoom9-17 i have one hillshade tiff:

<Style name="hillshade-90">
    <Rule>
        &maxscale_zoom9;
        &minscale_zoom17;
        <RasterSymbolizer comp-op="grain-merge" opacity="0.83" scaling="bilinear" />
    </Rule>
</Style>
<Layer name="hillshade-90">
    <StyleName>hillshade-90</StyleName>
    <Datasource>
        <Parameter name="type">gdal</Parameter>
        <Parameter name="file">dem/hillshade-90.tif</Parameter>
        <!--Parameter name="file">dem/hillshade-30m-jpeg.tif</Parameter-->
        <Parameter name="format">tiff</Parameter>
    </Datasource>
</Layer>

that one is made with

gdalwarp -co BIGTIFF=YES -co TILED=YES -co COMPRESS=LZW -co PREDICTOR=2 -t_srs "+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m" -r bilinear -tr 90 90 raw.tif warp-90.tif

and

gdaldem hillshade -z 2 -co compress=lzw -co predictor=2 -co bigtiff=yes -compute_edges warp-90.tif hillshade-90.tif

When i want to make a tiff for every zoomlevel:

9,10,11,12,13,14,15,16,17 what -tr, and -z values do i have to take?

a14stoner commented 1 year ago

Just to show you my setup here:

image

I have these containers running. i cache everything in mapcache because i have other 3-4 other wms sources which i all want to store at one place. The application then provides wms/wmts from all avaliable layers.

The server is a VM. it has 24 cores, 64gb ram and a 2 tb disk - the storage is based on NVMe so the speed should be good.

I checked the postgres db X times to change configuration and make the performance better. for more optimisations (changing mapnik config, ...) i think i am too unexperienced :D

yvecai commented 1 year ago

Basically you want to have the bigger tiff for the zoom level closest to your Dem resolution, normally zoom 11. Then another tiff with 1/2 resolution, for z10, 1/4 for z9, 1/8 for z8... Zoom level 12 and above, keep the full resolution for Mapnik to upsample. This is described in the readme.

Le 3 février 2023 07:50:24 GMT+01:00, a14stoner @.***> a écrit :

talking about downsampled tiffs for each zoom layer: Right now at zoom9-17 i have one hillshade tiff:

<Style name="hillshade-90">
  <Rule>
      &maxscale_zoom9;
      &minscale_zoom17;
      <RasterSymbolizer comp-op="grain-merge" opacity="0.83" scaling="bilinear" />
  </Rule>
</Style>
<Layer name="hillshade-90">
  <StyleName>hillshade-90</StyleName>
  <Datasource>
      <Parameter name="type">gdal</Parameter>
      <Parameter name="file">dem/hillshade-90.tif</Parameter>
      <!--Parameter name="file">dem/hillshade-30m-jpeg.tif</Parameter-->
      <Parameter name="format">tiff</Parameter>
  </Datasource>
</Layer>

that one is made with

gdalwarp -co BIGTIFF=YES -co TILED=YES -co COMPRESS=LZW -co PREDICTOR=2 -t_srs "+proj=merc +ellps=sphere +R=6378137 +a=6378137 +units=m" -r bilinear -tr 90 90 raw.tif warp-90.tif

and

gdaldem hillshade -z 2 -co compress=lzw -co predictor=2 -co bigtiff=yes -compute_edges warp-90.tif hillshade-90.tif

When i want to make a tiff for every zoomlevel:

9,10,11,12,13,14,15,16,17 what -tr, and -z values do i have to take?

-- Reply to this email directly or view it on GitHub: https://github.com/der-stefan/OpenTopoMap/issues/337#issuecomment-1415149174 You are receiving this because you were mentioned.

Message ID: @.***>

a14stoner commented 1 year ago

i only render above z10 -> z10,11,12,13,14,15,16,17

so the current configuration and tifs should be okay. 👍

a14stoner commented 1 year ago

seeded 151144 tiles (4696 skipped), now at z5 x368 y82 (i started from z10 because i dont need z0-9). So i am at z15 right now. until its finished it will take ages..... when it seeded since 2 weeks.

i see at docker stats:

CONTAINER ID   NAME          CPU %      MEM USAGE / LIMIT     MEM %     NET I/O           BLOCK I/O         PIDS
ff0eae54793a   mapcache      0.03%      104.7MiB / 62.74GiB   0.16%     4.42GB / 8.51GB   73.8GB / 7.81GB   130
9a58e2365cf8   postgis       0.00%      16.41GiB / 62.74GiB   26.15%    1.22TB / 155TB    2.1TB / 86.1GB    31
b59586da1bd7   mapproxy      1602.39%   40.04GiB / 62.74GiB   63.82%    155TB / 1.04TB    2.63TB / 572kB    311

that mapproxy (which includes mapnik) has the most resources. When i see this i think rendering images with mapnik takes the most time - not selecting data from database.

The question is now how to increase the speed of rendering the tiles.

Kradenko commented 2 months ago

Hi,

I have just finished building my server but when I run my sample set of data I also get these holes? I went through this thread and as I understand my data is bad and that is the reason I get these holes in the shadows?

Am I using incorrect data or am I doing something incorrectly? I am pretty stuck at the moment.

Any help/advice would be greatly appreciated.

Regards, D.

Kradenko commented 1 month ago

Hi,

I have just finished building my server but when I run my sample set of data I also get these holes? I went through this thread and as I understand my data is bad and that is the reason I get these holes in the shadows?

Am I using incorrect data or am I doing something incorrectly? I am pretty stuck at the moment.

Any help/advice would be greatly appreciated.

Regards, D.

Hi again,

To add to the above, I have also now tried to see maybe if I could somehow manipulate gdal_fillnodata to possibly fill these gaps. Without any success. I have tried to use the "-md " flag and increase the distance it checks but the resulting renders stay the same.

The only thing I can think of now is to get the other data you mention in your posts. Is there any way I can somehow make this work correctly?

Thanks!

Kradenko commented 1 month ago

@a14stoner could you perhaps guide me as to how you resolved this problem? I am completely stuck :(

Thanks in advance!

a14stoner commented 1 month ago

Hi @Kradenko I used some different data for hillshading. Since its a year ago i need to check which data I used ..

Kradenko commented 1 month ago

Hi @Kradenko I used some different data for hillshading. Since its a year ago i need to check which data I used ..

@a14stoner Awesome, thank you!

Kradenko commented 1 month ago

Hi @Kradenko I used some different data for hillshading. Since its a year ago i need to check which data I used ..

Hi @a14stoner any luck with those data files?

Kradenko commented 1 week ago

Hi @Kradenko I used some different data for hillshading. Since its a year ago i need to check which data I used ..

Hi @a14stoner @sletuffe @yvecai is there anyone that can help me with this alternate dataset? I have had a go at it again and cannot find any way to make this work.

Please help :)

yvecai commented 1 week ago

@Kradenko, this issue is quite old and honestly it's hard to figure out what you're doing, with which data, and what is your issue. How about describing your complete process and where you are stuck at on GisStackexchange for instance ?

Kradenko commented 1 week ago

@Kradenko, this issue is quite old and honestly it's hard to figure out what you're doing, with which data, and what is your issue. How about describing your complete process and where you are stuck at on GisStackexchange for instance ?

Thank you for the advice! Will go post there with what I have done and hopefully I get a solution :)