profile tool

[gubuntu]

allow pro users (members) to digitise a line and return a profile that can also be saved as a pdf.

What to use as base DEM?

supersedes #4

[gubuntu]

Here's the solution https://github.com/arno974/lizmap-altiProfil

Do we already have SRTM for SA in a PostGIS DB? We could use that directly else load it into the SAGTA PostGIS backend.

[vermeulendivan]

@gubuntu Do we have a DEM yet? From what it seems we don't. I've created a SRTM DEM mosaic which covers the entire SA a few months ago "for fun" on my own time. Its still at 30m spatial resolution, about 2gb ins size, but pretty sure we can use better storage/format.

[NyakudyaA]

@gubuntu Do we have a DEM yet? From what it seems we don't. I've created a SRTM DEM mosaic which covers the entire SA a few months ago "for fun" on my own time. Its still at 30m spatial resolution, about 2gb ins size, but pretty sure we can use better storage/format.

We already have this also in maps.kartoza but can you clip to SA extent and convert to COG and I will upload it to the cloud storage

[gubuntu]

We have various DEMs in various places, @NyakudyaA probably knows best. Read the plugin specifications carefully. If working in an area covered by IGN you would use that service, otherwise it needs a DEM loaded into PostGIS. If we already have one in PG then try to use that. I know we used to to have one to support a DEM service we ran for MyFarmWeb. But if not, you're welcome to use yours.

Are you able to set up the extension or do you need @NyakudyaA or @dominicschaff to help? Etienne, one of the developers, used to work for us so you can also ask him if we get stuck.

[gubuntu]

No CoG, from my reading of the docs it has to be in PG.

I think we should be able to use the one in the maps.kartoza.com backend if it is already in PG (or is it just a tiff?)

[NyakudyaA]

I haven't read the docs for the plugin but we do not have a raster in the DB. I would rather we spin up a PostgreSQL Db with this stack and ingest the raster because we want the request to be very fast

[vermeulendivan]

@NyakudyaA Its already clipped/masked to the coastal boundary (love using dark mode on qgis), but will perform clipping to the northern SA border as well.

[NyakudyaA]

@vermeulendivan https://maps.kartoza.com/geoserver/kartoza/wms?service=WMS&version=1.1.0&request=GetMap&layers=kartoza%3Asouth_africa&bbox=16.451890000000276%2C-34.834170000000086%2C32.944980000000214%2C-22.125029999999242&width=768&height=591&srs=EPSG%3A4326&format=application/openlayers

[NyakudyaA]

In the docs they mention the raster type should be   MNT. @gubuntu do you know how to configure this on ingestion to the DB

[gubuntu]

I think it's just French for DEM

[NyakudyaA]

Still debugging why the request are failing

[NyakudyaA]

Just update here. I am getting somewhere now.

Still the error is coming from the query that is generated in the backend to actually draw the profile

2022-04-14 16:59:11 172.20.0.1  error   Spatialite is not available
2022-04-14 16:59:36 172.20.0.1  error   2022-04-14 16:59:36 [403]   invalid query (ERROR:  step size cannot equal zero(
            WITH
                line AS(
                    -- From an arbitrary line
                    SELECT
                        ST_MakeLine(
                            ST_Transform(ST_SetSRID(ST_MakePoint(24.558229, -28.676768),4326), 4326),
                            ST_Transform(ST_SetSRID(ST_MakePoint(24.607667, -28.695742),4326), 4326)
                        )
                    AS geom
                ), 
                linemesure AS(
                    -- Add a mesure dimension to extract steps
                    SELECT
                        ST_AddMeasure(line.geom, 0, ST_Length(line.geom)) as linem,
                        generate_series(
                            0,
                            ST_Length(line.geom)::int,
                            --for very long line we reduce the steps
                            CASE
                                WHEN ST_Length(line.geom)::int < 1000 THEN 0
                                ELSE 0*5
                            END
                        ) as i,
                        CASE
                            WHEN ST_Length(line.geom)::int < 1000 THEN 0
                            ELSE 0*5
                        END as resolution
                    FROM line
                ), 
                points2d AS (
                    SELECT ST_GeometryN(ST_LocateAlong(linem, i), 1) AS geom, resolution FROM linemesure
                ),
                cells AS (
                    -- Get DEM elevation for each
                    SELECT 
                        p.geom AS geom, 
                        ST_Value(dem.rast, 1, p.geom) AS val,                      
                        resolution
                    FROM dem, points2d p
                    WHERE ST_Intersects(dem.rast, p.geom)
                ),           
                -- Instantiate 3D points
                points3d AS (
                    SELECT ST_SetSRID(
                                ST_MakePoint(ST_X(geom), ST_Y(geom), val),
                                4326
                            ) AS geom, resolution FROM cells
                ),
                line3D AS(
                    SELECT ST_MakeLine(geom)as geom, MAX(resolution) as resolution FROM points3d
                ),
                xz AS(
                    SELECT (ST_DumpPoints(geom)).geom AS geom,
                    ST_StartPoint(geom) AS origin, resolution
                    FROM line3D
                )
            -- Build 3D line from 3D points
            SELECT ST_distance(origin, geom) AS x, ST_Z(geom) as y, ST_X(geom) as lon, ST_Y(geom) as lat, resolution FROM xz))  /www/lib/jelix/plugins/db/pgsql/pgsql.dbconnection.php  187

[NyakudyaA]

Update: Partially working now but since to be an issue with using a DEM in 4326

[NyakudyaA]

Update this seems to be working now but the SQL queries for generating the profile seem to be a bit slow:

TODO

• It seems the profile tool is globally activated for the whole project and we only need to activate it for the members
• Add option for generating the profile into the print template

[NyakudyaA]

Update: I have optimised the SQL that generates points along the elevation

WITH
        profile_line as (
            select 
            ST_MakeLine(
                            ST_Transform(ST_SetSRID(ST_MakePoint(24.505014, -28.634891),4326), 3857),
                            ST_Transform(ST_SetSRID(ST_MakePoint(24.829797, -28.730670),4326), 3857)) as geom
                ),
                linemesure AS(
                    -- Add a mesure dimension to extract steps
                    SELECT

                        CASE
                            WHEN ST_Length(geom)::int < 1000 THEN 30
                            ELSE 30*5
                        END as resolution, geom 
                    from profile_line

                ),

                transect AS (
                     SELECT ST_Segmentize(
                        geom,
                       resolution
                        )::geometry
                      AS geom,resolution
                    from linemesure
                    ),
                    rast AS (
                     SELECT ST_Union(south_africa.rast) AS rast, resolution
                      FROM south_africa
                      JOIN transect
                      ON ST_Intersects(transect.geom, ST_ConvexHull(south_africa.rast))
                        group by transect.resolution
                    ),
                    z AS (
                      SELECT rast.resolution,(ST_DumpPoints(ST_SetZ(rast.rast, transect.geom, resample => 'bilinear'))).*
                     FROM rast
                      CROSS JOIN transect
                    ),
                    cells AS (
                        SELECT round(ST_Z(geom)) AS val, geom, resolution
                    FROM z
                    ),
                    -- Instantiate 3D points
                    points3d AS (
                        SELECT ST_SetSRID(
                                    ST_MakePoint(ST_X(geom), ST_Y(geom), val),
                                    3857
                                ) AS geom, resolution FROM cells
                    ),
                    line3D AS(
                        SELECT ST_MakeLine(geom)as geom, MAX(resolution) as resolution FROM points3d
                    ),
                    xz AS(
                        SELECT (ST_DumpPoints(geom)).geom AS geom,
                        ST_StartPoint(geom) AS origin, resolution
                        FROM line3D
                    )
                -- Build 3D line from 3D points
                SELECT ST_distance(origin, geom) AS x, ST_Z(geom) as y, ST_X(ST_Transform(geom, 4326)) as lon, ST_Y(ST_Transform(geom, 4326)) as lat, resolution FROM xz

This should speed up generating profiles no matter how big the line is. Will add this to the profile tool

[NyakudyaA]

This has been deployed to both staging and production