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

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.}}]


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.

• Here is another comparison shot showing the slight difference between the regions: i.stack.imgur.com/5kL5b.jpg – Jason B. Aug 19 '16 at 12:44
• Some information in general on projections from TU Delft here. For transportation modeling I once had a similar challenge which could be solved by Gauss-Krüger - there is a general link there to the epsg which may be helpful. Here in Germany all cities should rather use a similar projection - should that not be rather standardized for NL as well? – gwr Aug 22 '16 at 20:24
• For completeness this is the official (bilingual) description of the geodetic systems used in the Netherlands. – gwr Aug 23 '16 at 11:27
• I thank gwr, JasonB, and jose for their time and efforts. This was more difficult than I had anticipated. – Romke Bontekoe Aug 25 '16 at 10:29
• You may have noticed that my map of Amsterdam (above) has coordinates (via Frame->True) which do not match the latitude, and perhaps longitude, of the city (also above). Any suggestions? – Romke Bontekoe Aug 25 '16 at 10:34

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
]


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.

• The parameter offsetWGS84 was of course derived using the conversion tool cited in my answer. One might try more data points but I found that it seems reasonably close. – gwr Aug 23 '16 at 14:53
• this looks impressive! – Romke Bontekoe Aug 23 '16 at 20:13
• I have now implemented options as to allow Grid- and LatLon-Offsets as suggested by @jose. The differences remain invisible to my eye though. – gwr Aug 24 '16 at 11:06
• Please see my Edit 2. – Romke Bontekoe Aug 26 '16 at 8:20

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.

• @RomkeBontekoe - I was able to get the proper transformation. – Jason B. Aug 22 '16 at 21:26
• This looks very good! – Romke Bontekoe Aug 23 '16 at 8:01

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?

• Thanks for this solution. It is certainly close, but, as you also noted, it appears somewhat displaced. I fail to see where your WGS84 comes from. – Romke Bontekoe Aug 18 '16 at 5:52
• The WL maps use lat, lon coordinates on the ITRF00 datum, which is nearly identical to WGS84. Any other datum may use an ellipsoid with different orientation or center, and that can affect the meaning of latitude and longitude. Then, separately, we need the cartographic projection to convert those lat, lon values into the x, y coordinates of the flat map. The projection specification I gave handles the second part, but assumes that there is no datum change. – jose Aug 18 '16 at 9:43
• Comparing with the code given by JasonB, the results coincide to a few parts per million for the test polygon, which corresponds to sub-meter precision. How large is the displacement you see, Romke? – jose Aug 18 '16 at 9:46
• @jose - This is great, the result is visibly indistinguishable to what I put below. I wish there was more information on the parameters you can and can't feed to the projection argument of GeoGridPosition. BTW, would there be any easier way of making that plot than using replacement rules on the Polygon and FilledCurve[{{Line[... objects? This is what I came up with – Jason B. Aug 18 '16 at 14:32
• Using the fact that we are dealing with a small region of the world, I think the simplest solution here is to implement a constant offset as suggested by @gwr. However, I would implement a constant offset in projected space, not in latlon space. Hence, change the "GridOrigin" parameter in the projection I gave above to {155025.14, 463116.91}. This approximates the results given by JasonB. – jose Aug 23 '16 at 18:10