# Mismatch with GeoGraphics plotting of SHP data

Posted 5 months ago
920 Views
|
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: "CoordinateSystemInformation"-> PROJCS->{Belge_Lambert_1972,GEOGCS->{GCS_Belge_1972,DATUM-> {Belge_1972,SPHEROID->{International_1924,6.37839*10^6,297.}},PRIMEM-> {Greenwich,0.},UNIT->{Degree,0.0174533}},PROJECTION-> {Lambert_Conformal_Conic_2SP},PARAMETER->{False_Easting,150000.},PARAMETER-> {False_Northing,5.40009*10^6},PARAMETER->{Central_Meridian,4.36749},PARAMETER-> {Standard_Parallel_1,51.1667},PARAMETER->{Standard_Parallel_2,49.8333},PARAMETER-> {Latitude_Of_Origin,90.},UNIT->{Meter,1.}} Projection-> {"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.  Polygon[a__]:>Polygon@GeoPosition@GeoGridPosition[a,projb2] 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. Attachments:
4 Replies
Sort By:
Posted 4 months ago
 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"}  Attachments: