Re: OGRCreateCoordinateTransformation()
Even Rouault <[email protected]> Thu, 24 Mar 2022 16:50:22 +0100
| Newsgroups | gmane.comp.gis.gdal.devel,gmane.comp.gis.proj-4.devel |
|---|---|
| Message-ID | <[email protected]> |
This is a multi-part message in MIME format. --===============7111350555620532915== Content-Type: multipart/alternative; boundary="------------yyBcjwL2KFtxKOylpMmxep9u" Content-Language: en-US This is a multi-part message in MIME format. --------------yyBcjwL2KFtxKOylpMmxep9u Content-Type: text/plain; charset=UTF-8; format=flowed Content-Transfer-Encoding: 8bit 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 <[email protected]> > *Sent:* Wednesday, March 23, 2022 6:06 PM > *To:* Steve Riddell <[email protected]>; [email protected] > *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 <[email protected]> > <mailto:[email protected]>; [email protected] > *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 <[email protected]> > *Sent:* Friday, March 18, 2022 5:05 PM > *To:* Steve Riddell <[email protected]>; [email protected] > *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. --------------yyBcjwL2KFtxKOylpMmxep9u Content-Type: text/html; charset=UTF-8 Content-Transfer-Encoding: 8bit <html> <head> <meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> </head> <body> <p>Steve,</p> <p>(adding proj mailing list since it's actually purely a PROJ topic<br> <o:p></o:p></p> <blockquote type="cite" cite="mid:BN6PR01MB2404D77116699E1C7F55B982A6199@BN6PR01MB2404.prod.exchangelabs.com"> <div class="WordSection1"> <p class="MsoNormal"><o:p> </o:p></p> <p class="MsoNormal">Does the fact that there isn’t a consensus mean 3126 & 3127 aren’t going to end up in a PROJ release?</p> </div> </blockquote> <p> <br> </p> <p>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.<br> </p> <blockquote type="cite" cite="mid:BN6PR01MB2404D77116699E1C7F55B982A6199@BN6PR01MB2404.prod.exchangelabs.com"> <div class="WordSection1"> <p class="MsoNormal"><o:p></o:p></p> <p class="MsoNormal"><o:p> </o:p></p> <p class="MsoNormal">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? <br> </p> </div> </blockquote> 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. <br> <blockquote type="cite" cite="mid:BN6PR01MB2404D77116699E1C7F55B982A6199@BN6PR01MB2404.prod.exchangelabs.com"> <div class="WordSection1"> <p class="MsoNormal"><o:p></o:p></p> <p class="MsoNormal"><o:p> </o:p></p> <p class="MsoNormal">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?</p> </div> </blockquote> WKT2 hardly helps for that. The closest equivalent is the GEOIDMODEL[] sub node of a VERTCRS[] (<a class="moz-txt-link-freetext" href="http://docs.opengeospatial.org/is/18-010r7/18-010r7.html#83">http://docs.opengeospatial.org/is/18-010r7/18-010r7.html#83</a>). 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. <br> <blockquote type="cite" cite="mid:BN6PR01MB2404D77116699E1C7F55B982A6199@BN6PR01MB2404.prod.exchangelabs.com"> <div class="WordSection1"> <p class="MsoNormal"><o:p></o:p></p> <p class="MsoNormal"><o:p> </o:p></p> <p class="MsoNormal">Thanks,<o:p></o:p></p> <p class="MsoNormal">steve<o:p></o:p></p> <p class="MsoNormal"><o:p> </o:p></p> <div> <div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in"> <p class="MsoNormal"><b>From:</b> Even Rouault <a class="moz-txt-link-rfc2396E" href="mailto:[email protected]"><[email protected]></a> <br> <b>Sent:</b> Wednesday, March 23, 2022 6:06 PM<br> <b>To:</b> Steve Riddell <a class="moz-txt-link-rfc2396E" href="mailto:[email protected]"><[email protected]></a>; <a class="moz-txt-link-abbreviated" href="mailto:[email protected]">[email protected]</a><br> <b>Subject:</b> Re: OGRCreateCoordinateTransformation()<o:p></o:p></p> </div> </div> <p class="MsoNormal"><o:p> </o:p></p> <p>Steve,<o:p></o:p></p> <p>you may try PR <a href="https://github.com/OSGeo/PROJ/pull/3126" moz-do-not-send="true" class="moz-txt-link-freetext">https://github.com/OSGeo/PROJ/pull/3126</a> or <a href="https://github.com/OSGeo/PROJ/pull/3127" moz-do-not-send="true" class="moz-txt-link-freetext">https://github.com/OSGeo/PROJ/pull/3127</a>, 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. <o:p></o:p></p> <p>Even<o:p></o:p></p> <div> <p class="MsoNormal">Le 23/03/2022 à 23:42, Steve Riddell a écrit :<o:p></o:p></p> </div> <blockquote style="margin-top:5.0pt;margin-bottom:5.0pt"> <p class="MsoNormal">Hi Even,<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">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. <o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">COMPD_CS["TempTM + CGVD28 height - HT2_0",<o:p></o:p></p> <p class="MsoNormal"> PROJCS["Custom",<o:p></o:p></p> <p class="MsoNormal"> GEOGCS["NAD83(CSRS)",<o:p></o:p></p> <p class="MsoNormal"> DATUM["NAD83_Canadian_Spatial_Reference_System",<o:p></o:p></p> <p class="MsoNormal"> SPHEROID["GRS 1980",6378137,298.257222101,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","7019"]],<o:p></o:p></p> <p class="MsoNormal"> TOWGS84[0,0,0,0,0,0,0],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","6140"]],<o:p></o:p></p> <p class="MsoNormal"> PRIMEM["Greenwich",0,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","8901"]],<o:p></o:p></p> <p class="MsoNormal"> UNIT["degree",0.0174532925199433,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","9122"]],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","4617"]],<o:p></o:p></p> <p class="MsoNormal"> PROJECTION["Transverse_Mercator"],<o:p></o:p></p> <p class="MsoNormal"> PARAMETER["latitude_of_origin",49.351346659616],<o:p></o:p></p> <p class="MsoNormal"> PARAMETER["central_meridian",-123.20266499149],<o:p></o:p></p> <p class="MsoNormal"> PARAMETER["scale_factor",1],<o:p></o:p></p> <p class="MsoNormal"> PARAMETER["false_easting",15307.188],<o:p></o:p></p> <p class="MsoNormal"> PARAMETER["false_northing",6540.975],<o:p></o:p></p> <p class="MsoNormal"> UNIT["Meters",1],<o:p></o:p></p> <p class="MsoNormal"> AXIS["Easting",EAST],<o:p></o:p></p> <p class="MsoNormal"> AXIS["Northing",NORTH]],<o:p></o:p></p> <p class="MsoNormal"> VERT_CS["CGVD28 height - HT2_0",<o:p></o:p></p> <p class="MsoNormal"> VERT_DATUM["Canadian Geodetic Vertical Datum of 1928",2005,<o:p></o:p></p> <p class="MsoNormal"> EXTENSION["PROJ4_GRIDS","HT2_0.gtx"],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","5114"]],<o:p></o:p></p> <p class="MsoNormal"> UNIT["metre",1,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","9001"]],<o:p></o:p></p> <p class="MsoNormal"> AXIS["Gravity-related height",UP],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","5713"]]]<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">Experimenting, I removed the “TOWGS84” directive from BOTH the source and target systems, and I do get an answer instead of an error:<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">For test point: (<span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:black",serif">49.351346658, -123.202664992, 289.809)</span><o:p></o:p></p> <p class="MsoNormal"><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:black",serif">I get: (15307.1880,6540.9748,307.744)</span><o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">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:<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal" style="text-autospace:none"><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:black",serif">strPipe </span><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:teal",serif">=</span><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:black",serif"> </span><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:#A31515",serif">"+proj=pipeline\</span><o:p></o:p></p> <p class="MsoNormal" style="text-autospace:none"><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:#A31515",serif"> +step +proj=axisswap +order=2, 1\</span><o:p></o:p></p> <p class="MsoNormal" style="text-autospace:none"><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:#A31515",serif"> +step +proj=unitconvert +xy_in=deg +xy_out=rad\</span><o:p></o:p></p> <p class="MsoNormal" style="text-autospace:none"><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:#A31515",serif"> +step +inv +proj=vgridshift +grids=HT2_0.gtx +multiplier=1\</span><o:p></o:p></p> <p class="MsoNormal" style="text-autospace:none"><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:#A31515",serif"> +step +proj=tmerc +lat_0=49.351346659616 +lon_0=-123.20266499149 +k=1\</span><o:p></o:p></p> <p class="MsoNormal"><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:#A31515",serif"> +x_0=15307.188 +y_0=6540.975 +ellps=GRS80"</span><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:black",serif">;</span><o:p></o:p></p> <p class="MsoNormal"><span style="font-size:9.5pt;font-family:"Cascadia Mono ;color:black",serif"> </span><o:p></o:p></p> <p class="MsoNormal">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.<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">Any arrows left in your quiver?<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">Thanks,<o:p></o:p></p> <p class="MsoNormal">Steve<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <div> <div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in"> <p class="MsoNormal"><b>From:</b> Steve Riddell <br> <b>Sent:</b> Friday, March 18, 2022 6:17 PM<br> <b>To:</b> Even Rouault <a href="mailto:[email protected]" moz-do-not-send="true"><[email protected]></a>; <a href="mailto:[email protected]" moz-do-not-send="true" class="moz-txt-link-freetext">[email protected]</a><br> <b>Subject:</b> RE: OGRCreateCoordinateTransformation()<o:p></o:p></p> </div> </div> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">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.<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <div> <div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in"> <p class="MsoNormal"><b>From:</b> Even Rouault <<a href="mailto:[email protected]" moz-do-not-send="true" class="moz-txt-link-freetext">[email protected]</a>> <br> <b>Sent:</b> Friday, March 18, 2022 5:05 PM<br> <b>To:</b> Steve Riddell <<a href="mailto:[email protected]" moz-do-not-send="true" class="moz-txt-link-freetext">[email protected]</a>>; <a href="mailto:[email protected]" moz-do-not-send="true" class="moz-txt-link-freetext">[email protected]</a><br> <b>Subject:</b> Re: OGRCreateCoordinateTransformation()<o:p></o:p></p> </div> </div> <p class="MsoNormal"> <o:p></o:p></p> <p>Steve,<o:p></o:p></p> <p>Fix in PROJ queued in <a href="https://github.com/OSGeo/PROJ/pull/3123" moz-do-not-send="true" class="moz-txt-link-freetext">https://github.com/OSGeo/PROJ/pull/3123</a><o:p></o:p></p> <p>The complexity of dealing with legacy features TOWGS84[] and PROJ4_GRIDS is highly stressing PROJ pipeline computation engine.<o:p></o:p></p> <p>Using the plain EPSG codes EPSG:4955 -> EPSG:4617+5713 would be much preferable here<o:p></o:p></p> <p>Even<o:p></o:p></p> <div> <p class="MsoNormal">Le 18/03/2022 à 21:14, Steve Riddell a écrit :<o:p></o:p></p> </div> <blockquote style="margin-top:5.0pt;margin-bottom:5.0pt"> <p class="MsoNormal">Hi,<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">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.<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">COMPD_CS["NAD83(CSRS) + Ellipsoid (Meters)",<o:p></o:p></p> <p class="MsoNormal"> GEOGCS["NAD83(CSRS)",<o:p></o:p></p> <p class="MsoNormal"> DATUM["NAD83_Canadian_Spatial_Reference_System",<o:p></o:p></p> <p class="MsoNormal"> SPHEROID["GRS 1980",6378137,298.257222101,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","7019"]],<o:p></o:p></p> <p class="MsoNormal"> TOWGS84[0,0,0,0,0,0,0],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","6140"]],<o:p></o:p></p> <p class="MsoNormal"> PRIMEM["Greenwich",0,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","8901"]],<o:p></o:p></p> <p class="MsoNormal"> UNIT["degree",0.0174532925199433,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","9122"]],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","4617"]],<o:p></o:p></p> <p class="MsoNormal"> VERT_CS["Ellipsoid (Meters)",<o:p></o:p></p> <p class="MsoNormal"> VERT_DATUM["Ellipsoid",2002],<o:p></o:p></p> <p class="MsoNormal"> UNIT["metre",1,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","9001"]],<o:p></o:p></p> <p class="MsoNormal"> AXIS["ellipsoidal height",UP]]]<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">COMPD_CS["NAD83(CSRS) + CGVD28 height - HT2_0",<o:p></o:p></p> <p class="MsoNormal"> GEOGCS["NAD83(CSRS)",<o:p></o:p></p> <p class="MsoNormal"> DATUM["NAD83_Canadian_Spatial_Reference_System",<o:p></o:p></p> <p class="MsoNormal"> SPHEROID["GRS 1980",6378137,298.257222101,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","7019"]],<o:p></o:p></p> <p class="MsoNormal"> TOWGS84[0,0,0,0,0,0,0],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","6140"]],<o:p></o:p></p> <p class="MsoNormal"> PRIMEM["Greenwich",0,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","8901"]],<o:p></o:p></p> <p class="MsoNormal"> UNIT["degree",0.0174532925199433,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","9122"]],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","4617"]],<o:p></o:p></p> <p class="MsoNormal"> VERT_CS["CGVD28 height - HT2_0",<o:p></o:p></p> <p class="MsoNormal"> VERT_DATUM["Canadian Geodetic Vertical Datum of 1928",2005,<o:p></o:p></p> <p class="MsoNormal"> EXTENSION["PROJ4_GRIDS","HT2_0.gtx"],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","5114"]],<o:p></o:p></p> <p class="MsoNormal"> UNIT["metre",1,<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","9001"]],<o:p></o:p></p> <p class="MsoNormal"> AXIS["Gravity-related height",UP],<o:p></o:p></p> <p class="MsoNormal"> AUTHORITY["EPSG","5713"]]]<o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal">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]]]]]'<o:p></o:p></p> <pre> <o:p></o:p></pre> <pre><span style="font-family:"Calibri",sans-serif">Thanks in advance for any assistance!</span><o:p></o:p></pre> <pre><span style="font-family:"Calibri",sans-serif">Best regards,</span><o:p></o:p></pre> <pre><span style="font-family:"Calibri",sans-serif">Steve</span><o:p></o:p></pre> <pre><span style="font-family:"Calibri",sans-serif"> </span><o:p></o:p></pre> <p class="MsoNormal"> <o:p></o:p></p> <p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Tahoma",sans-serif">Steve Riddell</span><o:p></o:p></p> <p class="MsoNormal"><span style="font-family:"Tahoma",sans-serif">Software Engineering</span><o:p></o:p></p> <p class="MsoNormal"><b><span style="color:#525252">GeoCue Group, Inc.</span></b><o:p></o:p></p> <p class="MsoNormal"><span style="font-size:10.0pt">520 6th Street | Madison, AL 35756 USA</span><o:p></o:p></p> <p class="MsoNormal"><b><span style="font-size:10.0pt">Phone:</span></b><span style="font-size:10.0pt"> 256.461.8289 | <b>Fax:</b> 256.461.8249 </span><o:p></o:p></p> <p class="MsoNormal"><i><span style="color:#0070C0">LIDAR/Drone Mapping Software & Services </span></i><i><span style="font-size:9.0pt">– </span><span style="color:#538135">True View® 3D Imaging Sensors </span></i><i><span style="font-size:9.0pt">–</span><span style="color:#C45911"> Image/LIDAR Data Management Solutions</span></i><o:p></o:p></p> <p class="MsoNormal"><span style="font-size:10.0pt"> </span><a href="http://www.geocue.com/" moz-do-not-send="true"><span style="font-size:10.0pt;color:#0070C0">www.geocue.com</span></a><span style="font-size:10.0pt;color:#0070C0"> <a href="https://support.geocue.com/" moz-do-not-send="true">support.geocue.com</a> </span><o:p></o:p></p> <p class="MsoNormal"><span style="font-size:10.0pt"> </span><o:p></o:p></p> <p class="MsoNormal"> <o:p></o:p></p> <pre><span style="font-size:11.0pt;font-family:"Calibri",sans-serif"> </span><o:p></o:p></pre> </blockquote> <pre>-- <o:p></o:p></pre> <pre><a href="http://www.spatialys.com" moz-do-not-send="true" class="moz-txt-link-freetext">http://www.spatialys.com</a><o:p></o:p></pre> <pre>My software is free, but my time generally not.<o:p></o:p></pre> </blockquote> <pre>-- <o:p></o:p></pre> <pre><a href="http://www.spatialys.com" moz-do-not-send="true" class="moz-txt-link-freetext">http://www.spatialys.com</a><o:p></o:p></pre> <pre>My software is free, but my time generally not.<o:p></o:p></pre> </div> </blockquote> <pre class="moz-signature" cols="72">-- <a class="moz-txt-link-freetext" href="http://www.spatialys.com">http://www.spatialys.com</a> My software is free, but my time generally not.</pre> </body> </html> --------------yyBcjwL2KFtxKOylpMmxep9u-- --===============7111350555620532915== Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev --===============7111350555620532915==--