[gdal-dev] 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/gdal-dev/attachments/20220324/87ccbf45/attachment-0001.html>
More information about the gdal-dev
mailing list