der-stefan / OpenTopoMap

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

contour lines are missing #338

Closed a14stoner closed 1 year ago

a14stoner commented 1 year ago

Hi, I am rendering a map with DEMs from here: https://e4ftl01.cr.usgs.gov/ASTT/ASTGTM.003/2000.03.01/ --> ASTER DEMs from 2019.

I followed the tutorial from here: https://github.com/der-stefan/OpenTopoMap/blob/master/mapnik/HOWTO_DEM.md

The result is the following:

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 renderingit but the contour lines are missing Here is an example of one database entry:

image

I changed nothing at the opentopomap.xml or in the styles-otm/contours.xml. Here a little part from opentopomap.xml:

<Layer name="contours">
    <StyleName>contours</StyleName>
    <Datasource>
        <Parameter name="table">(SELECT way,ele FROM planet_osm_line) AS contours </Parameter>
        <!--<Parameter name="table">(SELECT way,ele FROM contours) AS contours </Parameter>-->
        <Parameter name="dbname">contours</Parameter>
        <Parameter name="type"><![CDATA[postgis]]></Parameter>
        <Parameter name="host"><![CDATA[postgis]]></Parameter>
        <Parameter name="port"><![CDATA[5432]]></Parameter>
        <Parameter name="password"><![CDATA[]]></Parameter>
        <Parameter name="user"><![CDATA[gis]]></Parameter>
        <Parameter name="estimate_extend"><![CDATA[false]]></Parameter>
        &postgis-settings;
        </Datasource>
</Layer>

As you can see i use the table planet_osm_line.

Thanks in advance.

sletuffe commented 1 year ago

To build your contour lines, you are using : https://github.com/der-stefan/OpenTopoMap/blob/master/mapnik/README.md with the phyghtmap + pbf + osm2pgsql tool ?

Well, I've not used this step, I prefered the more common gdal_contour tool, it can open a tif DEM file and insert directly into postgres the resulting contour lines in one step.

a14stoner commented 1 year ago

Okay then ill try it this way and report back the results ;) thank you

gdal_contour -f PostgreSQL -i 1 -nln contours warp-90.tif "PG:host=postgis user=gis dbname=contours"

trying it with these settings. Lets wait :D it takes a long time.. :D

after it i will update opentopomap.xml and set the configuration to that:

<Layer name="contours">
    <StyleName>contours</StyleName>
    <Datasource>
        <!--<Parameter name="table">(SELECT way,ele FROM planet_osm_line) AS contours </Parameter>-->
        <Parameter name="table">(SELECT way,ele FROM contours) AS contours </Parameter>
        <Parameter name="dbname">contours</Parameter>
        <Parameter name="type"><![CDATA[postgis]]></Parameter>
        <Parameter name="host"><![CDATA[postgis]]></Parameter>
        <Parameter name="port"><![CDATA[5432]]></Parameter>
        <Parameter name="password"><![CDATA[]]></Parameter>
        <Parameter name="user"><![CDATA[gis]]></Parameter>
        <Parameter name="estimate_extend"><![CDATA[false]]></Parameter>
        &postgis-settings;
        </Datasource>
</Layer>

That the contours table is used.

sletuffe commented 1 year ago

And you might or might not be interested, but I have an PostgreSQL dump of worldwide contours (It saves me a lot of pain every time I want to re-import my db) It is based on a mix or SRTM and ASTER, in the EPSG:3857 projection. https://maps.refuges.info/download/

It was generated during dozen of hours with a command like the one you used (but with -i 10 : you are too greedy !)

a14stoner commented 1 year ago

ohh ! nice to know! ill try it by myself before - it it doesnt work ill use your dump! THANKS!!!

a14stoner commented 1 year ago

Okay i imported it with gdal_contours with -i 10 because it took so long. I had to rename the colum wbk_geometry to ele because the column is named like this in the opentopomap.xml.

The result is the following:

image

no contours lines... can it be that the styles or even the opentopomap.xml is wrong? i made no other changes in the file...

sletuffe commented 1 year ago

Hard to guess what's wrong.

Check that your contour lines in PG are good : select ele,st_astext(way) from contours limit 10 ?

Is it empty ? Is the projection the expected spherical mercator EPSG:3857 ?

a14stoner commented 1 year ago

Since i imported the warp-90.tif with gdal_contours the table looks like this:

image

ogc_fid, id, wkb_geometry (which i renamed to ele).

image

made some checks and it seems i have the wrong coordinate system.

ill try to recreate it with that is has EPSG:3857

#################################################################################

Just another info: Thats the gdalsrsinfo output from the raw.tif:

PROJ.4 : +proj=longlat +datum=WGS84 +no_defs

OGC WKT2:2018 :
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

so EPSG:4326

Then when i execute the gdal_warp:

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

i tell gdalwarp that the t_srs is the proj4 string from EPSG:3857 butr the result is this:

PROJ.4 : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs

OGC WKT2:2018 :
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["unknown",
            ELLIPSOID["unknown",6378137,0,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["Mercator (variant A)",
        METHOD["Mercator (variant A)",
            ID["EPSG",9804]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

is that okay? shouldnt there stand 3857?

#################################################################################

i made a test - i did the gdalwarp with -t_srs EPSG:3857 and not with the proj4 string. Then when i imported it with gdal_contours --- the result of st_srid:

image

maybe that helps.

a14stoner commented 1 year ago

I think i found the problem,

when importing the data with gdal_contours the field wkb_geometry has to be renamed to way. But i dont have the values for the field ele. so the filters dont work.

ill try to download your dump from https://maps.refuges.info/download/ . Maybe it works with this..

a14stoner commented 1 year ago

Another thing i tried

  1. removing filters from styles-otm/contours.xml
  2. Renaming wkb_geometry to way
  3. Adding column ele to contours table

image

now they are here. to just importing it with gdal_contours is not enough. @sletuffe did you always import it with gdal_contours and it worked out of the box?

I dont know from where i should get the elevation tag which is now missing.

sletuffe commented 1 year ago

ogc_fid, id, wkb_geometry (which i renamed to ele).

That's a first problem : ele stands for "elevation" : it should hold the elevation of the line. The geometry of the (multi-) polyline being in the wkb_geometry column. For some reason I don't remember, mine is renamed to "geom" but I think you can rename the column to the one in the opentopomap.xml (way)

sletuffe commented 1 year ago

now they are here. to just importing it with gdal_contours is not enough. @sletuffe did you always import it with gdal_contours and it worked out of the box?

I dont know from where i should get the elevation tag which is now missing.

You are almost there ! Nice ;-)

To have the elevation, you forgot the "-a ele" switch to gdal_contour (I haven't paid attention at your first command)

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

-a <name>

    Provides a name for the attribute in which to put the elevation. If not provided no elevation attribute is attached.
sletuffe commented 1 year ago

As a side note : You seam to allways show screen captures of an openlayer slippy map, wich mean that you most likely need to :

For testing purpose I prefere https://github.com/springmeyer/nik2img (it's old, but still does the job) : it generates a png file from the mapnik style. Much faster than with all other components. And of course, a sample of contours and tiff that is small enough so all tasks run very fast.

a14stoner commented 1 year ago

Thanks for the suggestion. What i do: edit files, import db - test

I have a mapproxy docker container running which gets the data via mapnik and postgresql. - it doesnt cache so its good i dont have to delete anything.

That works like a charm also : D just the importing of the data takes long.

########################################################################

I reloaded the tif with gdal_contours (this time with -a ele ...:D) and the map looks like this now..

image

The only good thing about when something doesnt work is that i finally understand how this all works:

Big thanks to @sletuffe !!! really nice that you take time for me an help me finish making this map..

sletuffe commented 1 year ago

\o/

You are right, understanding what is behind the scene makes it easier to later understand the possible problems.

But some times downloading preprocess data that takes aaaaages to generate is also a good way to save time that others have spent :-) While still knowning what was needed to do it.