Coordinate offset after transformation from 4267 to 4326

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Coordinate offset after transformation from 4267 to 4326

Betsy Emmons
Hello, 

I'm using Proj.4 v4.9.1 and GeoTools v13.x. When a WMS layer is brought into QGIS as EPSG:4326 there is a noticeable offset compared to the same local layer. The WMS layer source is EPSG:4267. 

A WMS 4326 sample point has a coordinate of -97.2341220, 28.0611696 and the same local 4326 sample point has a coordinate of -97.2339531, 28.0608752. This is roughly a 40m difference. Some local calculations have shown the error margin to be <=111m.

Here is the source 4267 projection definition: GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]. 

Here's what geotools is using:
GEOGCS["NAD27", 
  DATUM["North American Datum 1927", 
    SPHEROID["Clarke 1866", 6378206.4, 294.9786982138982, AUTHORITY["EPSG","7008"]], 
    TOWGS84[2.478, 149.752, 197.726, 0.526, -0.498, 0.501, 0.685], 
    AUTHORITY["EPSG","6267"]], 
  PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
  UNIT["degree", 0.017453292519943295], 
  AXIS["Geodetic longitude", EAST], 
  AXIS["Geodetic latitude", NORTH], 
  AUTHORITY["EPSG","4267"]]

Proj.4 is really a great library! Thank you for maintaining and contributing to it! 

Thank you, 
Betsy Emmons


Betsy Emmons
Community Support Manager | Boundless
[hidden email]
<a href="tel:908-208-9692" value="+17322996769" style="color:rgb(17,85,204)" target="_blank">908.208.9692


_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate offset after transformation from 4267 to 4326

Martin Desruisseaux-3
Hello Betsy

The difference that you get may depend on whether a datum shift has been
applied or not. If both Proj.4 and GeoTools have applied a datum shift,
then it may depend on the values in the TOWGS84 expression. There is
many different possible TOWGS84 values for a transformations between
NAD27 and WGS84, depending on the geographic area. If you are
transforming coordinates in USA, you should have TOWGS84[-8, 160, 176]
(for an error below 10 metres). If you are transforming coordinates in
Canada, you should have TOWGS84[-10, 158, 187] (for an error below 20
metres). The TOWGS84 parameters in the GeoTools CRS provided with your
email are for Cuba (operation code EPSG:15978) except the sign of 0.526
which should be -0.526.

In order to compare with Proj.4, we would need to know which parameters
Proj.4 has selected. It may also use the NADCON grids in USA, or NTv2
grids in Canada. I do not know how to get this information from Proj.4.

Alternatively you could also try Apache SIS 0.7 for coordinate
transformations. It takes in account the geographic area of the points
to transform, support GML and WKT version 1 and 2, tell you which
parameters it selected, in which geographic area they are valid and what
accuracy to expect, etc. Version 0.7 should be released this week.

    Martin


Le 24/05/16 à 18:11, Betsy Emmons a écrit :

> Hello,
>
> I'm using Proj.4 v4.9.1 and GeoTools v13.x. When a WMS layer is
> brought into QGIS as EPSG:4326 there is a noticeable offset compared
> to the same local layer. The WMS layer source is EPSG:4267.
>
> A WMS 4326 sample point has a coordinate of -97.2341220, 28.0611696
> and the same local 4326 sample point has a coordinate of -97.2339531,
> 28.0608752. This is roughly a 40m difference. Some local calculations
> have shown the error margin to be <=111m.
>
> Here is the source 4267 projection
> definition: GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]].
>
> Here's what geotools is using:
> GEOGCS["NAD27",
>   DATUM["North American Datum 1927",
>     SPHEROID["Clarke 1866", 6378206.4, 294.9786982138982,
> AUTHORITY["EPSG","7008"]],
>     TOWGS84[2.478, 149.752, 197.726, 0.526, -0.498, 0.501, 0.685],
>     AUTHORITY["EPSG","6267"]],
>   PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
>   UNIT["degree", 0.017453292519943295],
>   AXIS["Geodetic longitude", EAST],
>   AXIS["Geodetic latitude", NORTH],
>   AUTHORITY["EPSG","4267"]]
>
> Proj.4 is really a great library! Thank you for maintaining and
> contributing to it!
>
> Thank you,
> Betsy Emmons

_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate offset after transformation from 4267 to 4326

jody.garnett
G'Day Martin:

The GeoTools has been configured with GridShft files for NAD27 - this is not something I have personally tried before so I am not quite sure how it works.

I had a look with Betsey, the data being projected is in Texas (so slightly different from Cuba so the TOWGS84 parameters sound legit).


--
Jody Garnett

On 24 May 2016 at 17:51, Martin Desruisseaux <[hidden email]> wrote:
Hello Betsy

The difference that you get may depend on whether a datum shift has been
applied or not. If both Proj.4 and GeoTools have applied a datum shift,
then it may depend on the values in the TOWGS84 expression. There is
many different possible TOWGS84 values for a transformations between
NAD27 and WGS84, depending on the geographic area. If you are
transforming coordinates in USA, you should have TOWGS84[-8, 160, 176]
(for an error below 10 metres). If you are transforming coordinates in
Canada, you should have TOWGS84[-10, 158, 187] (for an error below 20
metres). The TOWGS84 parameters in the GeoTools CRS provided with your
email are for Cuba (operation code EPSG:15978) except the sign of 0.526
which should be -0.526.

In order to compare with Proj.4, we would need to know which parameters
Proj.4 has selected. It may also use the NADCON grids in USA, or NTv2
grids in Canada. I do not know how to get this information from Proj.4.

Alternatively you could also try Apache SIS 0.7 for coordinate
transformations. It takes in account the geographic area of the points
to transform, support GML and WKT version 1 and 2, tell you which
parameters it selected, in which geographic area they are valid and what
accuracy to expect, etc. Version 0.7 should be released this week.

    Martin


Le 24/05/16 à 18:11, Betsy Emmons a écrit :
> Hello,
>
> I'm using Proj.4 v4.9.1 and GeoTools v13.x. When a WMS layer is
> brought into QGIS as EPSG:4326 there is a noticeable offset compared
> to the same local layer. The WMS layer source is EPSG:4267.
>
> A WMS 4326 sample point has a coordinate of -97.2341220, 28.0611696
> and the same local 4326 sample point has a coordinate of -97.2339531,
> 28.0608752. This is roughly a 40m difference. Some local calculations
> have shown the error margin to be <=111m.
>
> Here is the source 4267 projection
> definition: GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]].
>
> Here's what geotools is using:
> GEOGCS["NAD27",
>   DATUM["North American Datum 1927",
>     SPHEROID["Clarke 1866", 6378206.4, 294.9786982138982,
> AUTHORITY["EPSG","7008"]],
>     TOWGS84[2.478, 149.752, 197.726, 0.526, -0.498, 0.501, 0.685],
>     AUTHORITY["EPSG","6267"]],
>   PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
>   UNIT["degree", 0.017453292519943295],
>   AXIS["Geodetic longitude", EAST],
>   AXIS["Geodetic latitude", NORTH],
>   AUTHORITY["EPSG","4267"]]
>
> Proj.4 is really a great library! Thank you for maintaining and
> contributing to it!
>
> Thank you,
> Betsy Emmons

_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs


_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate offset after transformation from 4267 to 4326

Martin Desruisseaux-3

Hello Jody

The NADCON grid files are from NAD27 to NAD83. The coordinate transformations that Betsey wanted to do are from NAD27 to WGS84, which is considered as a different operation in the EPSG database. This difference was often ignored because we were used to consider NAD83 as equivalent to WGS84. This was approximatively true 20 years ago, but there is now a difference of about 1.5 metres between NAD83 and latest realizations of WGS84.

The key point is that GeoTools will not necessarily provide the same transformation for NAD27 to WGS84 than for NAD27 to NAD83. It depends on what GeoTools finds in the EPSG database, and how the referencing module is designed. Last time I touched to the GeoTools referencing module, it had a bug: when it found more than one operation for the same sourceCRS and targetCRS, GeoTools unconditionally selected the most accurate one, regardless its area of validity. This explain why the TOWGS84 parameters are for Cuba: this is the operation with the smallest "positional accuracy" found in the EPSG database, but this accuracy is true only in a relatively small geographic area.

The above bug exists (or existed) because at the time I wrote the referencing module, I though that a sourceCRS and a targetCRS were sufficient for uniquely identifying a coordinate operation. I didn't realized that we can have as much as 80 different operations for the same source and target CRS depending on the geographic area. The fix for GeoTools would be that the CRS.findMathTransform(...) method takes another argument, which is the geographic area where the user intend to apply the coordinate operation. This problem has been fixed in Apache SIS (among other issues discovered as I continue to learn new aspects).

Even if GeoTools uses the grids, the grids are not the same for USA than for Canada, so a hint about geographic area is still desirable. For Texas, if the desired accuracy is 1 metre then there is two possible operations depending on whether the points are on the west side or the east side of 100°W: EPSG:8624 and EPSG:8625. If an accuracy of 10 metres is sufficient then the TOWGS84 parameters for USA given in my previous emails should be okay.

As a side note: the above also explain why the TOWGS84 element in WKT has been removed in the new WKT 2 specification (ISO 19162): the TOWGS84 element can not handle the real-world complexity and is actually dangerous since it can be misleading when given without information about its domain of validity, like what happened with the TOWGS84 parameters for Cuba.

    Martin


Le 26/05/16 à 12:39, Jody Garnett a écrit :
G'Day Martin:

The GeoTools has been configured with GridShft files for NAD27 - this is not something I have personally tried before so I am not quite sure how it works.

I had a look with Betsey, the data being projected is in Texas (so slightly different from Cuba so the TOWGS84 parameters sound legit).


_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate offset after transformation from 4267 to 4326

Martin Desruisseaux-3
I forgot to said: in order to verify what GeoTools actually used - i.e.
whether it used the grids or the TOWGS84 parameters, and in the later
case which parameter values -, a System.out.println(theMathTransform)
would help.

    Martin

_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate offset after transformation from 4267 to 4326

jody.garnett
Thanks for your help Martin, I will see what I can do.

I would still like to learn if there is a way to check what PROJ4 is up to, as this is the first time I have worked with grid shift files I have no confidence I am doing so correctly.

--
Jody Garnett

On 26 May 2016 at 11:19, Martin Desruisseaux <[hidden email]> wrote:
I forgot to said: in order to verify what GeoTools actually used - i.e.
whether it used the grids or the TOWGS84 parameters, and in the later
case which parameter values -, a System.out.println(theMathTransform)
would help.

    Martin



_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate offset after transformation from 4267 to 4326

Martin Desruisseaux-3
Le 26/05/16 à 14:56, Finn, Michael a écrit :

> Thank you for this detailed explanation.

My pleasure :-)


Le 26/05/16 à 15:12, Jody Garnett a écrit :

> I would still like to learn if there is a way to check what PROJ4 is
> up to, as this is the first time I have worked with grid shift files I
> have no confidence I am doing so correctly.

This would indeed be a useful information. In my understanding (I may be
wrong), PROJ4 is what EPSG calls an "early-binding implementation",
meaning that it relies on the information attached to the CRS like the
TOWGS84 element except maybe in a few special cases. Early binding
implementations are often designed on the assumption that some CRS
(typically WGS84) can be used as a "hub" or a "pivot CRS" for all
coordinate transformations. Whether such assumption is okay or not
depends on the desired accuracy.

By contrast, a "late-binding implementation" (in EPSG terminology)
searches information about coordinate operations in the EPSG database
only when the full (sourceCRS, targetCRS,
other-aspects-like-area-of-interest) tupple is known. So the
transformation from A to B is not necessarily the same than the
transformation from A to WGS84 to B. GeoTools and Apache SIS are
late-binding implementations.

EPSG recommends late-binding implementations, but those implementations
are not easy to get right. In any cases, PROJ4 and GeoTools
implementations are different enough that in order to perform comparison
between them, we need to know the transformation paths (operation method
and parameters) that both of them selected.

    Martin


_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate offset after transformation from 4267 to 4326

jody.garnett
In reply to this post by jody.garnett
Just wanted to report back to this email thread, although I was never able to determine quite what PROJ4 was doing but we were able to use a series of points to make an educated guess:

The fix, at least for me, was to define a customer epsg_operations.properties file that combined the grid-shift transform (from NAD27 --> NAD83) and the shift from NAD83 --> WGS84):

      4267,4326=CONCAT_MT[ \
         PARAM_MT["NADCON", \
           PARAMETER["Latitude difference file", "conus.las"], \
           PARAMETER["Longitude difference file", "conus.los"]], \
         PARAM_MT["Molodenski", \
           PARAMETER["dim", 2], \
           PARAMETER["dx", 0.0], \
           PARAMETER["dy", 0.0], \
           PARAMETER["dz", 0.0], \
           PARAMETER["src_semi_major", 6378137.0], \
           PARAMETER["src_semi_minor", 6356752.314140356], \
           PARAMETER["tgt_semi_major", 6378137.0], \
           PARAMETER["tgt_semi_minor", 6356752.314245179]]]

This has output between GeoServer and QGIS lining up sufficient for the application I am working on.
--
Jody

--
Jody Garnett

On 26 May 2016 at 12:12, Jody Garnett <[hidden email]> wrote:
Thanks for your help Martin, I will see what I can do.

I would still like to learn if there is a way to check what PROJ4 is up to, as this is the first time I have worked with grid shift files I have no confidence I am doing so correctly.

--
Jody Garnett

On 26 May 2016 at 11:19, Martin Desruisseaux <[hidden email]> wrote:
I forgot to said: in order to verify what GeoTools actually used - i.e.
whether it used the grids or the TOWGS84 parameters, and in the later
case which parameter values -, a System.out.println(theMathTransform)
would help.

    Martin




_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate offset after transformation from 4267 to 4326

Martin Desruisseaux-3

Thanks for the update, this make sense.

For info, EPSG geodetic dataset defines a transformation from NAD27 to NAD83 using the "conus" files as you did, followed by a transformation from NAD83 to NAD83(HARN) (one of the new NAD83 definitions created after the original one) one the assumption that NAD83(HARN) is equivalent to WGS84 within the accuracy of the transformation. The later transformation is using other grids depending on the geographic area (I didn't saw a grid for the whole US for the NAD83 to NAD83(HARN) part, but maybe I didn't searched hard enough).

    Martin


Le 14/06/16 à 01:15, Jody Garnett a écrit :
Just wanted to report back to this email thread, although I was never able to determine quite what PROJ4 was doing but we were able to use a series of points to make an educated guess:

The fix, at least for me, was to define a customer epsg_operations.properties file that combined the grid-shift transform (from NAD27 --> NAD83) and the shift from NAD83 --> WGS84):

      4267,4326=CONCAT_MT[ \
         PARAM_MT["NADCON", \
           PARAMETER["Latitude difference file", "conus.las"], \
           PARAMETER["Longitude difference file", "conus.los"]], \
         PARAM_MT["Molodenski", \
           PARAMETER["dim", 2], \
           PARAMETER["dx", 0.0], \
           PARAMETER["dy", 0.0], \
           PARAMETER["dz", 0.0], \
           PARAMETER["src_semi_major", 6378137.0], \
           PARAMETER["src_semi_minor", 6356752.314140356], \
           PARAMETER["tgt_semi_major", 6378137.0], \
           PARAMETER["tgt_semi_minor", 6356752.314245179]]]

This has output between GeoServer and QGIS lining up sufficient for the application I am working on.
--
Jody


_______________________________________________
MetaCRS mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/metacrs