[gdal-dev] GDAL(Technical support for Extracting Elevations

Sun Nov 8 16:50:02 PST 2020

Greetings,

Our team decided to use GDAL in our code to support certain business requirements. Each time we do an analysis we need the elevation profile lines, our analyses run through 10,000s of such cases across a large area. Which means we need to extract these elevation profile lines 10,000s per session.
Currently, we have a partially-working solution in C# (we're using C# bindings for GDAL) but when we load several DTED files and increase their resolution we run into "OutOfMemoryExceptions". We are  "manually" extract all of the elevations out of these files into large lists in memory (10's of millions of elevations, spread over large area). How would you go about making sure that you could do this quickly and efficiently? We need technical guidance on how to support the objectives we outlined below within the GDAL/OGR framework.

Objectives:
1) We have a large set of terrain files in DTED format, we want to load these into the application
2) We want to potentially warp the terrain data (increase resolution or crop the extents of the terrain data) which they have loaded into the application
3) Once the terrain data is loaded into the application, we want to:
                                                                               * extract elevation profile lines from the terrain data (i.e., the list of elevations, in meters, between two input lat/lon coordinates).
                                                                               * extract elevation, in meters, at a single input lat/lon coordinate.

Regards,

-----Original Message-----
From: gdal-dev <gdal-dev-bounces at lists.osgeo.org> On Behalf Of gdal-dev-request at lists.osgeo.org
Sent: Friday, November 6, 2020 3:00 PM
To: gdal-dev at lists.osgeo.org
Subject: gdal-dev Digest, Vol 198, Issue 9

EXTERNAL

Send gdal-dev mailing list submissions to
        gdal-dev at lists.osgeo.org

To subscribe or unsubscribe via the World Wide Web, visit
        https://lists.osgeo.org/mailman/listinfo/gdal-dev
or, via email, send a message with subject or body 'help' to
        gdal-dev-request at lists.osgeo.org

You can reach the person managing the list at
        gdal-dev-owner at lists.osgeo.org

When replying, please edit your Subject line so it is more specific
than "Re: Contents of gdal-dev digest..."

Today's Topics:

   1. Re: Spherical Mercator Projection mismatching major/minor
      axis (jratike80)
   2. Re: Spherical Mercator Projection mismatching major/minor
      axis (Tobias Wendorff)
   3. GDAL(Technical support for Extracting Elevations From Large
      DTED Terrain Datasets) (Nabil Adechoubou)
   4. Re: Spherical Mercator Projection mismatching major/minor
      axis (Craig Bruce)
   5. Re: Spherical Mercator Projection mismatching     major/minor
      axis (Even Rouault)

----------------------------------------------------------------------

Message: 1
Date: Fri, 6 Nov 2020 04:07:43 -0700 (MST)
From: jratike80 <jukka.rahkonen at maanmittauslaitos.fi>
To: gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] Spherical Mercator Projection mismatching
        major/minor axis
Message-ID: <1604660863882-0.post at n6.nabble.com>
Content-Type: text/plain; charset=us-ascii

Hi,

This does not really belong to my knowledge area but I'll have a try anyway.

Check what you have after reading the proj string instead. Here with Python

>>> from osgeo import osr
>>> spatialRef = osr.SpatialReference()
>>> spatialRef.ImportFromProj4("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
>>> +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null
>>> +towgs84=0,0,0,0,0,0,0 +wktext +no_defs")
>>> print(spatialRef)
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0
+x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"]]

It means that GDAL makes a clever guess and considers that the Proj string
that you provided means EPSG:3857. You probably mean the same. And in the
EPSG database EPSG:3857 is based on WGS84 ellipsoid
https://epsg.org/crs/wkt/id/3857

 ELLIPSOID["WGS 84",6378137,298.257223563,

and that it behaves in spherical way is handled by a special conversion

CONVERSION["Popular Visualisation Pseudo-Mercator",

I am not sure what to do if you definitely want to keep the ball with
diameter of 6378137 meters but perhaps dropping +nadgrids=@null could do it.
If EPSG:3857 is what you really need then it takes less writing to import by
the EPSG code but your way gives the same result.

-Jukka Rahkonen-

Andreas Roth wrote
> Hi,
>
> I call OSRImportFromProj4 with
> "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
> +k=1.0 +units=m +nadgrids=@null +towgs84=0,0,0,0,0,0,0 +wktext +no_defs"
>
> and with the returned handle i get the major and minor axis (
> OSRGetSemiMajor/OSRGetSemiMinor). I would expect the both axis to be
> equal,
> but they aren't.
> Instead i'm getting the following result:
> semi_major_axis=6378137
> semi_minor_axis=6356752.314
>
> Build on Ubuntu 20.10 (groovy; amd64) with libgdal-dev,
> version 3.1.3+dfsg-1ubuntu2) and libproj-dev:amd64 with version 7.1.0-1
>
> The full example application code:
> https://pastebin.com/xf8QpVhM
>
> Regards,
> Andreas
>
> _______________________________________________
> gdal-dev mailing list

> gdal-dev at .osgeo

> https://lists.osgeo.org/mailman/listinfo/gdal-dev

--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html

------------------------------

Message: 2
Date: Fri, 6 Nov 2020 12:43:39 +0100
From: Tobias Wendorff <tobias.wendorff at tu-dortmund.de>
To: gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] Spherical Mercator Projection mismatching
        major/minor axis
Message-ID: <62743fff-48ce-3a5e-d361-bac9f29d49ee at tu-dortmund.de>
Content-Type: text/plain; charset=utf-8; format=flowed

Am 06.11.2020 um 12:07 schrieb jratike80:
> This does not really belong to my knowledge area but I'll have a try anyway.

To add a bit of my knowledge, EPSG:3857 is the king/queen of
unconformness (better unconform-mess):

The geographic reference is based on geographic coordinates in WGS84
(ellipsoid), the projection is based on a sphere. The background is that
the spherical math is way easier and it was therefor faster to render
(and easiert to code). You don't need to solve any integrals etc.

Wikipedia describes the problem:
https://en.wikipedia.org/wiki/Web_Mercator_projection

A discussion about that is appear at least once a year in OpenStreetMap
community, when a GIS specialist or scientist uses the data for the
first time ;)

------------------------------

Message: 3
Date: Fri, 6 Nov 2020 15:55:25 +0000
From: Nabil Adechoubou <Nabil.Adechoubou at validatek.com>
To: "gdal-dev at lists.osgeo.org" <gdal-dev at lists.osgeo.org>
Cc: Vladislav Martin <Vladislav.Martin at validatek.com>
Subject: [gdal-dev] GDAL(Technical support for Extracting Elevations
        From Large DTED Terrain Datasets)
Message-ID:
        <BN8PR18MB2482DEA0964F4ABD46A17AEB97ED0 at BN8PR18MB2482.namprd18.prod.outlook.com>

Content-Type: text/plain; charset="utf-8"

Greetings,

Our team decided to use GDAL in our code to support certain business requirements. Each time we do an analysis we need the elevation profile lines, our analyses run through 10,000s of such cases across a large area. Which means we need to extract these elevation profile lines 10,000s per session.
Currently, we have a partially-working solution in C# (we're using C# bindings for GDAL) but when we load several DTED files and increase their resolution we run into "OutOfMemoryExceptions". We are  "manually" extract all of the elevations out of these files into large lists in memory (10's of millions of elevations, spread over large area). How would you go about making sure that you could do this quickly and efficiently? We need technical guidance on how to support the objectives we outlined below within the GDAL/OGR framework.

Objectives:
1) We have a large set of terrain files in DTED format, we want to load these into the application
2) We want to potentially warp the terrain data (increase resolution or crop the extents of the terrain data) which they have loaded into the application
3) Once the terrain data is loaded into the application, we want to:
                                                                               * extract elevation profile lines from the terrain data (i.e., the list of elevations, in meters, between two input lat/lon coordinates).
                                                                               * extract elevation, in meters, at a single input lat/lon coordinate.

Regards,

Nabil Adechoubou
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201106/1878e694/attachment-0001.html>

------------------------------

Message: 4
Date: Fri, 6 Nov 2020 14:27:04 -0400
From: Craig Bruce <csbruce at cubewerx.com>
To: gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] Spherical Mercator Projection mismatching
        major/minor axis
Message-ID: <b7856f2a-f776-e0f0-eee4-d9da56073b2f at cubewerx.com>
Content-Type: text/plain; charset="utf-8"

An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201106/29efd27c/attachment-0001.html>

------------------------------

Message: 5
Date: Fri, 06 Nov 2020 20:24:22 +0100
From: Even Rouault <even.rouault at spatialys.com>
To: gdal-dev at lists.osgeo.org
Cc: Craig Bruce <csbruce at cubewerx.com>
Subject: Re: [gdal-dev] Spherical Mercator Projection mismatching
        major/minor axis
Message-ID: <5839595.iMEXY37nVm at even-i700>
Content-Type: text/plain; charset="us-ascii"

Craig,

> The way I handle it is (slightly older notation):
>
> PROJCS["WGS84 / Spherical Mercator",
>     GEOGCS["WGS84basedSpheric_GCS",
>         DATUM["WGS84basedSpheric_Datum",
>             SPHEROID["WGS84based_Sphere", 6378137, 0]],
>         PRIMEM["Greenwich", 0],
>         UNIT["degree", 0.0174532925199433]],
>     PROJECTION["Mercator"],
>     PARAMETER["latitude_of_origin", 0],
>     PARAMETER["central_meridian", 0],
>     PARAMETER["false_easting", 0],
>     PARAMETER["false_northing", 0],
>     UNIT["meter", 1],
>     AUTHORITY["EPSG", "3857"]]
>

The above WKT using plain Mercator and a sphere definition could potentially
result in incorrect datum transformation being done by some systems if
transforming coordinates from a non-WGS84 based CRS to this definition
(actually just testing it with PROJ >= 6 or GDAL <= 2.4, it rejects its
because it doesn't know the 'Mercator' projection method. Should be
'Mercator_1SP' for the WKT1 dialect that GDAL/PROJ speaks)

Jukka was quite correct: EPSG:3857 is defined on the WGS 84 datum, but the
projection method to use is spherical mercator, that is forcing flattening to
0. See https://epsg.org/crs/wkt/id/3857

> Will GDAL accept a "0" for the inverse of flattening to represent a sphere?

Yes. That's the canonical way of expressing a sphere

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

------------------------------

Subject: Digest Footer

_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

------------------------------

End of gdal-dev Digest, Vol 198, Issue 9
****************************************