GeoPosition coordinates from Shape file (.SHP) with unknown Datum

Question

The City of Amsterdam supports an Open Data policy. Much information can be retrieved at the level of boroughs (stadsdelen). Unfortunately, GeoGraphics does not provide outlines of these smaller areas. A bundle of Shape files (.SHP) of Administrative Areas can be downloaded.

I want to overlay these areas with the city map, but the .SHP data are in a coordinate system which is not recognized. (Related to 42533.)

Overall, the images of the administrative areas look good (no reverse coordinates)

but the coordinates are not in geographic {Latitude,Longitude}, since Amsterdam is at

CityData[Entity["City",{"Amsterdam","NoordHolland","Netherlands"}],"Coordinates"]

gives {52.37, 4.89}.

The coordinate system

Import[amsterdamDir, "CoordinateSystemInformation"]

gives

"PROJCS" -> {"RD_New", 
"GEOGCS" -> {"GCS_Amersfoort", 
"DATUM" -> {"Amersfoort", 
  "SPHEROID" -> {"Bessel_1841", 6.3774*10^6, 299.153}}, 
"PRIMEM" -> {"Greenwich", 0}, 
"UNIT" -> {"Degree", 0.0174532925199432955}}, 
"PROJECTION" -> {"Double_Stereographic"}, 
"PARAMETER" -> {"False_Easting", 155000}, 
"PARAMETER" -> {"False_Northing", 463000}, 
"PARAMETER" -> {"Central_Meridian", 5.38764}, 
"PARAMETER" -> {"Scale_Factor", 0.999908}, 
"PARAMETER" -> {"Latitude_Of_Origin", 52.1562}, 
"UNIT" -> {"Meter", 1}}

However, its Datum (GCSAmersfoort) is unknown. (Amersfoort is another city in The Netherlands.)

How can I transform these coordinates to GeoPosition coordinates?

For testing, a tiny polygon from the Shape file is

Polygon[{{115050., 485415.}, {115614., 485596.}, {115601., 
485634.}, {115887., 485725.}, {116049., 485777.}, {116052., 
485768.}, {116124., 485793.}, {116136., 485798.}, {116154., 
485809.}, {116162., 485816.}, {116171., 485825.}, {116262., 
485750.}, {116263., 485687.}, {116263., 485666.}, {116236., 
485665.}, {116206., 485662.}, {116148., 485649.}, {116122., 
485642.}, {115834., 485549.}, {115846., 485509.}, {115684., 
485456.}, {115676., 485481.}, {115578., 485449.}, {115633., 
485279.}, {115564., 485257.}, {115494., 485229.}, {115211., 
485138.}, {115146., 485116.}, {115061., 485378.}, {115050., 
485415.}}]

Thanks in advance.

EDIT 1

Following @jose's solution I obtain this map, with the above polygon overlaid

This shows a small area (Calandlaan/Lelylaan buurt) in Amsterdam-West. A screenshot from maps of the City of Amsterdam show this outline

The correspondence is close, but there is a small displacement to be resolved.

EDIT 2

The shift between the two solutions of @gwr is about 21 cm. Just by taking the GeoDistance[ ] between the two projections of the same coordinate.

pos1=rdPosition[{115061, 485378 },OffsetGrid->{+25.14, +116.91}, OffsetWGS84->{0, 0}]
pos2=rdPosition[{115061, 485378 },OffsetGrid->{+0, +0}, OffsetWGS84 ->{-0.001053524113977744`, -0.000352533797287613`}]

GeoDistance[{pos1,pos2}]

OUT:
GeoPosition[{52.3548, 4.80096}, "WGS84"]
GeoPosition[{52.3548, 4.80097}, "WGS84"]
Quantity[0.206247, "Meters"]

Indeed, too small to discern at the scale of the maps.

Answer 1

Since there is a rather constant offset between ITRF00 and WGS84 we might expand upon @jose's solution:

Clear[rdPosition];

rdPosition::usage = "\
rdPosition[{x,y}] takes a geocentric position according to the new \
Dutch RD system and returns a GeoPosition using WGS84";

Options[rdPosition] = {
    "OffsetWGS84" -> {0, 0}, (* {Δϕ,Δλ} to adjust ITRF00 -> WGS84 *)
    "OffsetGrid"  -> {+25.14, +116.91} (* {ΔEast,ΔNorth} to adjust GridOrigin *)
};

(* OffsetWGS84 of original answer was {-0.001053524113977744`, -0.000352533797287613`} *)

rdPosition[ pos : {x_?NumericQ, y_?NumericQ}, opts:OptionsPattern[rdPosition] ] := With[
  {
    falseEasting = 155000,
    falseNorthing = 463000,
    centering = {52.1562, 5.38764},
    centralScaleFactor = 0.999908,
    referenceEllipsoid = "Bessel1841",
    offsetGrid = OptionVAlue["OffsetGrid"],
    offsetWGS84 = OptionValue["OffsetWGS84"]
  },
  GeoGridPosition[
    pos,
    {
      "Stereographic", (* projection *)
      "ReferenceModel" -> referenceEllipsoid,
      "GridOrigin" -> {falseEasting, falseNorthing} + offsetGrid,
      "Centering" -> centering,
      "CentralScaleFactor" -> centralScaleFactor
    }
  ] // RightComposition[
         GeoPosition, (* using WL standard ITRF00 *)
         ReplaceAll[ GeoPosition[{ϕ_, λ_}] :> GeoPosition[{ϕ, λ} + offsetWGS84, "WGS84"] ]
       ]
]

Using the conversion tool here we can check this function:

rdPosition @ {115061, 485378 }

GeoPosition[{52.3548, 4.80096}, "WGS84"]

This seems close enough and we can now simply add a convenient function to convert shape data and use it:

Clear[ shapeConvert ];

shapeConvert::usage = "\
shapeConvert[data] converts the shape data to geographic primitives\
which can be used within GeoGraphics.";

shapeConvert[ shapeData_, opts:OptionsPattern[rdPosition] ] := ReplaceAll[
    shapeData,
    pos:{x_?NumericQ, y_?NumericQ} :> pos ~ rdPosition ~ opts
]

Administrative Areas

amstAreas = ReplaceAll[
    "Geometry", 
    Import[
      "http://maps.amsterdam.nl/gebiedsindeling/GEBIEDSINDELINGEN.ZIP",
      {"SHP", "Data"}
    ]
];

Row[{
  amstAreas[[-2, 1]] // shapeConvert[ #, 
      OffsetGrid -> {+25.14, +116.91}, 
      OffsetWGS84 -> {0, 0}
  ] & // GeoGraphics[{EdgeForm[Red], #}, ImageSize -> Large] & ,
  amstAreas[[-2, 1]] // shapeConvert[#, 
      OffsetGrid -> {+0, +0}, 
      OffsetWGS84 -> {-0.001053524113977744`, -0.000352533797287613`}
  ] & // GeoGraphics[{EdgeForm[Orange], #}, ImageSize -> Large] &
}]

Row[{
  amstAreas[[2, 1]] // shapeConvert[ #, 
      OffsetGrid -> {+25.14, +116.91}, 
      OffsetWGS84 -> {0, 0}
  ] & // GeoGraphics[{EdgeForm[Red], #}, ImageSize -> Large] & ,
  amstAreas[[2, 1]] // shapeConvert[#, 
      OffsetGrid -> {+0, +0}, 
      OffsetWGS84 -> {-0.001053524113977744`, -0.000352533797287613`}
  ] & // GeoGraphics[{EdgeForm[Orange], #}, ImageSize -> Large] &
}]

I fail to really spot a difference between using a Grid-Offset and a WGS84-Offset.

Update

I have now implemented two different possible offsets (Grid, WGS84) by using options to rdPosition. The options may also be passed to shapeConvert. The grid offset has been proposed by @jose in a comment to his answer.

Answer 2

I'm still hoping that jose might be able to make this work using GeoGridPosition or GeoGridPositionXYZ. Or maybe some ambitious person will use the equations on page 65 of this pdf to do the job manually.

But I'm just going to farm the job out to the PROJ.4 Cartographic Projections Library, which I installed by following the instructions here.

In this function, you will need to change the path to the cs2cs script for your local installation,

ClearAll[convertAmersfoort]
convertAmersfoort[data_?(ArrayQ[#, _, NumericQ] &)] := 
  Module[{cs2cs, file, result}, 
   cs2cs = "~/projects/land-registry/proj-4.9.1/bin/cs2cs";
   file = CreateTemporary[];
   Export[file, data, "Table"];
   result = 
    Import["!" <> cs2cs <> 
        " -f '%.7f' +proj=sterea +lat_0=52.15616055555555 \
+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 \
+ellps=bessel \
+towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 \
+units=m +no_defs +to +proj=latlong +datum=WGS84 -s " <> file, 
       "Table"] // Map@Most // GeoPosition;
   DeleteFile[file];
   result];
convertAmersfoort[pgons_] := 
  ReplaceAll[
   pgons, {Polygon[a__] :> Polygon@convertAmersfoort@a, 
    Line[a__] :> Line@convertAmersfoort@a}];

If pgon is the polygon defined at the end of the OP, then

GeoGraphics[{
  EdgeForm[Directive[Blue, Thick, Dashed]], convertAmersfoort[pgon]
  }
 ]

returns a faithful reproduction of the linked original,

You can apply this to all the polygons from the original SHP file. Here is another example of the comparison,

amstAreas = 
  "Geometry" /. 
   Import["http://maps.amsterdam.nl/gebiedsindeling/GEBIEDSINDELINGEN.ZIP", {"SHP", "Data"}];
GeoGraphics[{EdgeForm[Directive[Blue, Thick, Dashed]], 
  convertAmersfoort[amstAreas[[-2, 1]]]}]

versus that from the original site:

If there were method for using the 7-parameter "towgs84" datum-conversion, then using another program would not be necessary.

Answer 3

My approximation to that projection is

proj = {"Stereographic", "ReferenceModel" -> "Bessel1841", "GridOrigin" -> {155000, 463000}, "Centering" -> {52.1562, 5.38764}, "CentralScaleFactor" -> 0.999908}

However, I don't know how to relate the Amersfoort datum to WGS84, and this may result in a global displacement.

Calling pol the testing polygon given above, try this:

GeoGridPosition[First[pol], proj]

geopol = Polygon[GeoPosition[%]]

GeoGraphics[geopol]

Does the result look reasonable?