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 <[email protected]> On Behalf Of
[email protected]
Sent: Friday, November 6, 2020 3:00 PM
To: [email protected]
Subject: gdal-dev Digest, Vol 198, Issue 9
EXTERNAL
Send gdal-dev mailing list submissions to
[email protected]
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
[email protected]
You can reach the person managing the list at
[email protected]
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 <[email protected]>
To: [email protected]
Subject: Re: [gdal-dev] Spherical Mercator Projection mismatching
major/minor axis
Message-ID: <[email protected]>
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
> [email protected]
> 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 <[email protected]>
To: [email protected]
Subject: Re: [gdal-dev] Spherical Mercator Projection mismatching
major/minor axis
Message-ID: <[email protected]>
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 <[email protected]>
To: "[email protected]" <[email protected]>
Cc: Vladislav Martin <[email protected]>
Subject: [gdal-dev] GDAL(Technical support for Extracting Elevations
From Large DTED Terrain Datasets)
Message-ID:
<bn8pr18mb2482dea0964f4abd46a17aeb97...@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 <[email protected]>
To: [email protected]
Subject: Re: [gdal-dev] Spherical Mercator Projection mismatching
major/minor axis
Message-ID: <[email protected]>
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 <[email protected]>
To: [email protected]
Cc: Craig Bruce <[email protected]>
Subject: Re: [gdal-dev] Spherical Mercator Projection mismatching
major/minor axis
Message-ID: <5839595.iMEXY37nVm@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
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
------------------------------
End of gdal-dev Digest, Vol 198, Issue 9
****************************************
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
