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.

Here is another comparison shot showing the slight difference between the regions: http://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
Regarding the Frame-Coordinates given for the map you should take a look at the tutorial for GeoGraphics. There is a section called Projected Coordinates which explains what is going on. If you press the right button on you mouse, the menu gives get coordinates. Using this feature will show you, that everything in fact is in order: lat,lon - values will match. — gwr, Aug 26 '16 at 10:40

Jason B. · Answer 1 · 2016-08-22T19:12:55.653

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,

Mathematica graphics

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

Mathematica graphics

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

gwr · Accepted Answer · 2016-08-24T11:09:42.770

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.

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
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

score 7 · Answer 3 · answered Aug 17 '16 at 22:51

7

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?

answered Aug 17 '16 at 22:51

jose

6,328
1
14
24

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 05: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 09: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 09: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
For example, could it be possible for something like GeoGraphics[First@First@amstAreas, GeoProjection -> proj] to work? – Jason B. Aug 18 '16 at 14:42
@jose - Please see my edit. – Romke Bontekoe Aug 19 '16 at 10:38
@jose Is it true, that WL cannot transfer ITRF (or some other datum) to WGS84 - which is among the most relevant geodetic systems? The information how to do it seems to be public here it seems. What other conversions would people be interested in? :) – gwr Aug 23 '16 at 13:58
The Wolfram Language knows how to transform between many pairs of datums, using 7-parameter Helmert transformations (3 parameters for a translation, 3 parameters for a rotation and 1 parameter for a global rescaling). GeodesyData knows the relevant information, and GeoPosition/GeoPositionXYZ can access it. The difference between WGS84 and ITRF00 is small, and can be ignored. However, the parameters needed here are not small, judging by the towgs84 values given by @JasonB. – jose Aug 23 '16 at 17:48
@jose - is there any mechanism to provide the 7 parameters to GeoPositionXYZ manually, in addition to the centering, scaling, and origin values? – Jason B. Aug 23 '16 at 17:51
3

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

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

3 Answers3

Administrative Areas

Update

Linked