Search code examples
pythongiscoordinate-systemspyproj

PROJ pyproj conversion of point EPSG 4326 (WSG 84) to (EPSG 28992)

Given a (lon, lat) point (5.068913, 52.067567) I would like to convert from EPSG 4326 to EPSG 28992 using pyproj.

The Proj, and transform functions in pyproj both seem to be suitable for such a task:

When I use the Proj function I get a different result to using transform, why?

For example

from shapely.geometry import Point
from pyproj import Proj, transform

from matplotlib import pyplot as plt

x1, y1 = 5.068913, 52.067567

in_proj = Proj(init='epsg:4326') 
out_proj = Proj(init='epsg:28992')

point1 = Point(out_proj(x1, y1))
point2 = Point(transform(in_proj, out_proj, x1 ,y1))

print(point1 == point2)

fig, ax = plt.subplots()
x, y = point1.xy
ax.plot(x, y, 'ro')
x, y = point2.xy
ax.plot(x, y, 'ro')

points

Solution

  • The conversion between

    EPSG:4326

    and

    EPSG:28992

    requires not only projecting the input coordinates, using the Oblique_Stereographic projection, in this case. That can be accomplished with the Proj function...

    ...but also performing a datum transformation, WGS84 to Amersfoort datum. That can only be accomplished with the transform function. Note that transform does everything that Proj does plus datum transformations. There is no need for a two-step transformation.

    Illustrating what have been said, we have for EPSG:28992 the following WKT,

    PROJCS["Amersfoort / RD New",
    GEOGCS["Amersfoort",
        DATUM["Amersfoort",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725],
            AUTHORITY["EPSG","6289"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4289"]],
    PROJECTION["Oblique_Stereographic"],
    PARAMETER["latitude_of_origin",52.15616055555555],
    PARAMETER["central_meridian",5.38763888888889],
    PARAMETER["scale_factor",0.9999079],
    PARAMETER["false_easting",155000],
    PARAMETER["false_northing",463000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","28992"]]

    The string

    TOWGS84[565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725]

    ...implies the need to perform a datum conversion... that cannot be acomplished with the Projfunction, but only with transform.

    Conclusion:

    Assuming the pyproj installation you are using is configured with the correct TOWGS84 strings (that is not always the case), then the transform result should be the one to be considered correct.