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

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.

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?

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.