Message Boards Message Boards


Mismatch with GeoGraphics plotting of SHP data

Posted 5 months ago
4 Replies
4 Total Likes

I'm trying to work with geographic city-data from Belgium Statistics Service. I can extract all the different elements from the SHP set. But I had great trouble to plot some some Geographic data. Some of the definitions for the geo transformations are given below:


{"LambertConformalConic", "Centering" -> GeoPosition[{90., 4.36749}], 
"GridOrigin" -> {150000., 5.40009*10^6}, 
"StandardParallels" -> {51.1667, 49.8333}, 
"ReferenceModel" -> {Quantity[6.35691*10^6, "Meters"], 
Quantity[6.37839*10^6, "Meters"]}}

I used the geo data of the township of Aartselaar for these tets. The polygon coordinates are in the form of:

Polygon[{{151660.03130000085, 203013.60550000146},{.....

The polygon data must be transformed to adjust tot the coordinate system of the background maps.


But its difficult to find the all the meanings and syntaxes of the parameters to build a projection definition in GeoGridPosition. I tried to find a match between the elements of "Projection" (above) and

In[128]:= GeoProjectionData["LambertConicConformal"]

Out[128]= {"LambertConicConformal", {"Centering" -> {0, 0}, 
  "CentralScaleFactor" -> 1, "GridOrigin" -> {0, 0}, 
  "ReferenceModel" -> 1, "StandardParallels" -> {33, 45}}}

And merged those two in "projn2". This result in an error during an evaluation:

GeoGridPosition::invset: "ReferenceModel"->{6.35691*10^6m,6.37839*10^6m} is not a valid setting.

So I changed the definition of the ReferenceModel to "Bessl1841". I found this in an example with coordinates around Amsterdam. Resulting in below:

projb2 = {"LambertConformalConic", 
"Centering" -> GeoPosition[{90., 4.36748666666667}], 
"GridOrigin" -> {150000.013, 5.400088438*^6}, 
"StandardParallels" -> {51.1666672333333, 49.8333339}, 
"ReferenceModel" -> "Bessel1841"}

This works, but there is an offset between the polygons and the map (see att). I did a lot of trials on Belgium and Amsterdam data. I put everything together in SHP_4.nb. The output is deleted to limit the file length.


4 Replies

I found a correction after a lot of trial and error that solved the mismatch. I did some corrections on the GridOrigin in projb2. The problem is that I have no idea what "Bessel1841" means.

projb2 = {"LambertConformalConic", 
  "Centering" -> GeoPosition[{90., 4.36748666666667}], 
  "GridOrigin" -> {149910.013, 5.399265000*^6}, 
  "StandardParallels" -> {51.1666672333333, 49.8333339}, 
  "ReferenceModel" -> "Bessel1841"}
Posted 4 months ago

Does this help?

"Bessel1841" is a reference ellipsoid for the Earth. You can get its size from GeodesyData, for example:

In[39]:= GeodesyData["Bessel1841", "EllipsoidParameters"]
Out[39]= {Quantity[6.3774*10^6, "Meters"], {299.153}}

Note however that the information you provided above mentions another ellipsoid: "International1924", so I'd suggest to use that instead of "Bessel1841".

I used ellipsoid: "International1924". The result without Origin correction was much better. I measured the difference on the map:

GeoDisplacement[{115.944, -58.0622}]

Result was almost perfect with some GridOrigin correction, see att.:

belgProj = {"LambertConformalConic", 
  "Centering" -> GeoPosition[{90., 4.36748666666667}], 
  "GridOrigin" -> {149908.06924394477`, 5.400146500244812`*^6}, 
  "StandardParallels" -> {51.1666672333333, 49.8333339}, 
  "ReferenceModel" -> "International1924"}

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract