[PROJ] OGRCreateCoordinateTransformation()

Even Rouault even.rouault at spatialys.com
Thu Mar 24 08:50:22 PDT 2022


Steve,

(adding proj mailing list since it's actually purely a PROJ topic

> Does the fact that there isn’t a consensus mean 3126 & 3127 aren’t 
> going to end up in a PROJ release?
>

Did you try them to confirm they help ? Maybe if you want to weight in 
the PRs (if they actually help) that could change the balance.

> I read through the PR’s. Are you thinking my issue is that a vertical 
> shift is being introduced because of a transformation to WGS84 (from 
> NAD83(CSRS)) to look up the geoid separation value?
>
That's probably an issue with the wrong interpolation CRS being used, 
but I can't generally predict how this complex code works without 
launching a debugger and step-by-step tracing what happens.
>
> PROJ4_GRIDS provided a way to specify exactly which geoid file(s) to 
> use.  Some vertical datums, such as NAVD88, only have a single EPSG 
> code for multiple geoids – what’s the replacement for PROJ4_GRIDS in 
> PROJ9?  If I could switch to WKT2 strings (not sure that’s really 
> feasible) does it fix these problems?
>
WKT2 hardly helps for that. The closest equivalent is the GEOIDMODEL[] 
sub node of a VERTCRS[] 
(http://docs.opengeospatial.org/is/18-010r7/18-010r7.html#83). But 
that's something that refer to an existing transformation . And in a 
VRTCRS[] you can actually have several GEOIDMODEL[] listed (e.g. for 
NAVD88, all the geoid models), since that's more a definition of the 
VertCRS than which particular geoid model to use. So that's not 
nominally appropriate for a "random" grid not known of PROJ, although we 
have a bit bent the concept to allow a GEOIDMODEL["PROJ 
name_of_the_grid_dot_tif_or_gtx"]. But I would suspect that in that 
situation the ambiguity with the geographic / interpolation CRS to use 
would be similar to WKT1 PROJ4_GRIDS.
>
> Thanks,
>
> steve
>
> *From:* Even Rouault <even.rouault at spatialys.com>
> *Sent:* Wednesday, March 23, 2022 6:06 PM
> *To:* Steve Riddell <sriddell at geocue.com>; gdal-dev at lists.osgeo.org
> *Subject:* Re: OGRCreateCoordinateTransformation()
>
> Steve,
>
> you may try PR https://github.com/OSGeo/PROJ/pull/3126 or 
> https://github.com/OSGeo/PROJ/pull/3127, that might potentially help, 
> but they don't make consensus (see discussions of those PRs). 
> Basically it is a nightmare to deal with the old TOWGS84 / PROJ4_GRIDS 
> approach, that poorly specifies some elements of transformations 
> (especially the vertical part), and doesn't fit well in the ISO-19111 
> model used by newer PROJ (we use and abuse the BoundCRS approach in 
> complex ways). One of the issue too is that we've realized late in the 
> process that we didn't emulate some aspects of the old transformation 
> pipeline, and there isn't general enthusiasm to change that now.
>
> Even
>
> Le 23/03/2022 à 23:42, Steve Riddell a écrit :
>
>     Hi Even,
>
>     I tried the fix with the CRS’s I sent, and it worked. But as I
>     mentioned, my real target CRS is a custom projected system
>     (below), and it still fails with what seems to be the same error
>     as in my original email.
>
>     COMPD_CS["TempTM + CGVD28 height - HT2_0",
>
>         PROJCS["Custom",
>
>             GEOGCS["NAD83(CSRS)",
>
>     DATUM["NAD83_Canadian_Spatial_Reference_System",
>
>                     SPHEROID["GRS 1980",6378137,298.257222101,
>
>     AUTHORITY["EPSG","7019"]],
>
>                     TOWGS84[0,0,0,0,0,0,0],
>
>     AUTHORITY["EPSG","6140"]],
>
>                 PRIMEM["Greenwich",0,
>
>     AUTHORITY["EPSG","8901"]],
>
>     UNIT["degree",0.0174532925199433,
>
>     AUTHORITY["EPSG","9122"]],
>
>                 AUTHORITY["EPSG","4617"]],
>
>     PROJECTION["Transverse_Mercator"],
>
>     PARAMETER["latitude_of_origin",49.351346659616],
>
>     PARAMETER["central_meridian",-123.20266499149],
>
>             PARAMETER["scale_factor",1],
>
>     PARAMETER["false_easting",15307.188],
>
>     PARAMETER["false_northing",6540.975],
>
>             UNIT["Meters",1],
>
>             AXIS["Easting",EAST],
>
>             AXIS["Northing",NORTH]],
>
>         VERT_CS["CGVD28 height - HT2_0",
>
>             VERT_DATUM["Canadian Geodetic Vertical Datum of 1928",2005,
>
>     EXTENSION["PROJ4_GRIDS","HT2_0.gtx"],
>
>                 AUTHORITY["EPSG","5114"]],
>
>             UNIT["metre",1,
>
>                 AUTHORITY["EPSG","9001"]],
>
>             AXIS["Gravity-related height",UP],
>
>             AUTHORITY["EPSG","5713"]]]
>
>     Experimenting, I removed the “TOWGS84” directive from BOTH the
>     source and target systems, and I do get an answer instead of an error:
>
>     For test point: (49.351346658, -123.202664992, 289.809)
>
>     I get: (15307.1880,6540.9748,307.744)
>
>     But since geoid separation at the test point is -18.238, Z should
>     be 308.047.  In fact, if I code up the following pipeline and
>     initialize a transform object:
>
>     strPipe ="+proj=pipeline\
>
>             +step +proj=axisswap +order=2, 1\
>
>             +step +proj=unitconvert +xy_in=deg +xy_out=rad\
>
>             +step +inv +proj=vgridshift +grids=HT2_0.gtx +multiplier=1\
>
>             +step +proj=tmerc +lat_0=49.351346659616
>     +lon_0=-123.20266499149 +k=1\
>
>             +x_0=15307.188 +y_0=6540.975 +ellps=GRS80";
>
>     The call to Transform() gives: (15307.1880, 6540.9748, 308.042)  
>     -- XY matches above, but now I get the expected Z (-0.005).
>      However, I’d like to avoid assembling pipelines for every
>     conversion/transform.
>
>     Any arrows left in your quiver?
>
>     Thanks,
>
>     Steve
>
>     *From:* Steve Riddell
>     *Sent:* Friday, March 18, 2022 6:17 PM
>     *To:* Even Rouault <even.rouault at spatialys.com>
>     <mailto:even.rouault at spatialys.com>; gdal-dev at lists.osgeo.org
>     *Subject:* RE: OGRCreateCoordinateTransformation()
>
>     Thanks very much, Even.  For reporting the problem, I “packaged
>     up” a case with CRS’s with EPGS codes. For the actual problem I’m
>     working, the target CRS is a custom projection –no EPSG.
>
>     *From:* Even Rouault <even.rouault at spatialys.com>
>     *Sent:* Friday, March 18, 2022 5:05 PM
>     *To:* Steve Riddell <sriddell at geocue.com>; gdal-dev at lists.osgeo.org
>     *Subject:* Re: OGRCreateCoordinateTransformation()
>
>     Steve,
>
>     Fix in PROJ queued in https://github.com/OSGeo/PROJ/pull/3123
>
>     The complexity of dealing with legacy features TOWGS84[] and
>     PROJ4_GRIDS is highly stressing PROJ pipeline computation engine.
>
>     Using the plain EPSG codes  EPSG:4955 -> EPSG:4617+5713 would be
>     much preferable here
>
>     Even
>
>     Le 18/03/2022 à 21:14, Steve Riddell a écrit :
>
>         Hi,
>
>         I’m getting an error (below) using GDAL trying to set up a
>         transformation between NAD83(CSRS) geodetic, w/ ellipsoidal
>         heights, and NAD83(CSRS) geodetic, CGVD28 heights. Similar
>         transformations work fine on some datums,  but fail on others.
>          Hope someone can offer some insight.  Note the if I don't use
>         the PROJ4_GRIDS extension, I don't get the error, but I don't
>         get the geoid adjustment either.
>
>         COMPD_CS["NAD83(CSRS) + Ellipsoid (Meters)",
>
>             GEOGCS["NAD83(CSRS)",
>
>         DATUM["NAD83_Canadian_Spatial_Reference_System",
>
>                     SPHEROID["GRS 1980",6378137,298.257222101,
>
>         AUTHORITY["EPSG","7019"]],
>
>                     TOWGS84[0,0,0,0,0,0,0],
>
>                     AUTHORITY["EPSG","6140"]],
>
>                 PRIMEM["Greenwich",0,
>
>                     AUTHORITY["EPSG","8901"]],
>
>         UNIT["degree",0.0174532925199433,
>
>                     AUTHORITY["EPSG","9122"]],
>
>                 AUTHORITY["EPSG","4617"]],
>
>             VERT_CS["Ellipsoid (Meters)",
>
>                 VERT_DATUM["Ellipsoid",2002],
>
>                 UNIT["metre",1,
>
>                     AUTHORITY["EPSG","9001"]],
>
>                 AXIS["ellipsoidal height",UP]]]
>
>         COMPD_CS["NAD83(CSRS) + CGVD28 height - HT2_0",
>
>             GEOGCS["NAD83(CSRS)",
>
>         DATUM["NAD83_Canadian_Spatial_Reference_System",
>
>                     SPHEROID["GRS 1980",6378137,298.257222101,
>
>         AUTHORITY["EPSG","7019"]],
>
>                     TOWGS84[0,0,0,0,0,0,0],
>
>                     AUTHORITY["EPSG","6140"]],
>
>                 PRIMEM["Greenwich",0,
>
>                     AUTHORITY["EPSG","8901"]],
>
>         UNIT["degree",0.0174532925199433,
>
>                     AUTHORITY["EPSG","9122"]],
>
>                 AUTHORITY["EPSG","4617"]],
>
>             VERT_CS["CGVD28 height - HT2_0",
>
>                 VERT_DATUM["Canadian Geodetic Vertical Datum of
>         1928",2005,
>
>         EXTENSION["PROJ4_GRIDS","HT2_0.gtx"],
>
>                     AUTHORITY["EPSG","5114"]],
>
>                 UNIT["metre",1,
>
>                     AUTHORITY["EPSG","9001"]],
>
>                 AXIS["Gravity-related height",UP],
>
>                 AUTHORITY["EPSG","5713"]]]
>
>         ERROR 6: Cannot find coordinate operations from `EPSG:4955' to
>         `COMPOUNDCRS["NAD83(CSRS) + CGVD28 height -
>         HT2_0",BOUNDCRS[SOURCECRS[GEOGCRS["NAD83(CSRS)",DATUM["NAD83
>         Canadian Spatial Reference System",ELLIPSOID["GRS
>         1980",6378137,298.257222101,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",4617]]],TARGETCRS[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["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]],ABRIDGEDTRANSFORMATION["Transformation
>         from NAD83(CSRS) to WGS84",METHOD["Position Vector
>         transformation (geog2D
>         domain)",ID["EPSG",9606]],PARAMETER["X-axis
>         translation",0,ID["EPSG",8605]],PARAMETER["Y-axis
>         translation",0,ID["EPSG",8606]],PARAMETER["Z-axis
>         translation",0,ID["EPSG",8607]],PARAMETER["X-axis
>         rotation",0,ID["EPSG",8608]],PARAMETER["Y-axis
>         rotation",0,ID["EPSG",8609]],PARAMETER["Z-axis
>         rotation",0,ID["EPSG",8610]],PARAMETER["Scale
>         difference",1,ID["EPSG",8611]]]],BOUNDCRS[SOURCECRS[VERTCRS["CGVD28
>         height - HT2_0",VDATUM["Canadian Geodetic Vertical Datum of
>         1928"],CS[vertical,1],AXIS["gravity-related
>         height",up,LENGTHUNIT["metre",1]],ID["EPSG",5713]]],TARGETCRS[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,3],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal
>         height",up,ORDER[3],LENGTHUNIT["metre",1]],ID["EPSG",4979]]],ABRIDGEDTRANSFORMATION["CGVD28
>         height - HT2_0 to WGS84 ellipsoidal
>         height",METHOD["GravityRelatedHeight to
>         Geographic3D"],PARAMETERFILE["Geoid (height correction) model
>         file","HT2_0.gtx",ID["EPSG",8666]]]]]'
>
>           
>
>         Thanks in advance for any assistance!
>
>         Best regards,
>
>         Steve
>
>         Steve Riddell
>
>         Software Engineering
>
>         *GeoCue Group, Inc.*
>
>         520 6th Street | Madison, AL 35756 USA
>
>         *Phone:*256.461.8289 | *Fax:* 256.461.8249
>
>         /LIDAR/Drone Mapping Software & Services //– True View® 3D
>         Imaging Sensors //–Image/LIDAR Data Management Solutions/
>
>         www.geocue.com <http://www.geocue.com/>support.geocue.com
>         <https://support.geocue.com/>
>
>     -- 
>
>     http://www.spatialys.com
>
>     My software is free, but my time generally not.
>
> -- 
> http://www.spatialys.com
> My software is free, but my time generally not.

-- 
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20220324/87ccbf45/attachment-0001.html>


More information about the PROJ mailing list