Difference between revisions of "MGI / Balkans coordinate systems"

m (Validation of CRS parameters)
m
 
Line 12: Line 12:
 
* [[Image:logo_coord.gif|18px]] [http://spatial-analyst.net/CRS/croatia.csy croatia.csy] : HTRS96/TM coordinate system in ILWIS format.
 
* [[Image:logo_coord.gif|18px]] [http://spatial-analyst.net/CRS/croatia.csy croatia.csy] : HTRS96/TM coordinate system in ILWIS format.
 
* [[Image:Icon_zip.png|18px]] [http://spatial-analyst.net/DATA/grid_croatia_etrs_laea_1k.zip grid_croatia_etrs_laea_1k.shp] : European (ETRS) grid nodes for Croatia.
 
* [[Image:Icon_zip.png|18px]] [http://spatial-analyst.net/DATA/grid_croatia_etrs_laea_1k.zip grid_croatia_etrs_laea_1k.shp] : European (ETRS) grid nodes for Croatia.
 +
* [http://spatial-analyst.net/PDF/HR_geoid_transformation_parameters_per_county.pdf HR geoid transformation pars] : Transformation parameters per each county in Croatia.
 
}}  
 
}}  
  
Line 100: Line 101:
 
<br>  
 
<br>  
  
The difference from the original Bessel's ellipsoid (so-called [http://www.mapref.org/GeodeticReferenceSystemsAT.html ''Hermannskogel datum'']) to the current referent system ([http://en.wikipedia.org/wiki/WGS84 WGS84]) is non-systematic. The national land surveying agencies would typically set-up a network of referent points and then calculate the mean transformation parameters (local geoid) that are in average accurate enough for the most of the country (see map of the errors for the HRG2000 in figure below). There are in fact several sources for the global transformation parameters that can be used to link the local datum with the WGS84 system. The most accurate global transformation parameters are the one provided by the national land survey agency. For example, in Croatia, the global transformation parameters (HRG2000) have been estimated using 1780 referent points (the average 3D error is 0.7 meters) by the State Geodetic Authority. In Serbia, the global transformation parameters were determined based on 1217 referent points (the average 3D error is 0.4 meters).  
+
The difference from the original Bessel's ellipsoid (so-called [http://www.mapref.org/GeodeticReferenceSystemsAT.html ''Hermannskogel datum'']) to the current referent system ([http://en.wikipedia.org/wiki/WGS84 WGS84]) is non-systematic. The national land surveying agencies would typically set-up a network of referent points and then calculate the mean transformation parameters (local geoid) that are in average accurate enough for the most of the area of interest. There are in fact several sources for the global transformation parameters that can be used to link the local datum with the WGS84 system. The most accurate global transformation parameters are the one provided by the national land survey agency (see for example the [http://spatial-analyst.net/PDF/HR_geoid_transformation_parameters_per_county.pdf transformation parameters for each county in Croatia]). There are also global transformation parameters that can be used for a wider area, but then the mean positional accuracy can exceed 1 m. For example, In Croatia, transformation parameters for HRG2000 have been estimated using 1780 referent points (the average 3D error is 0.7 meters; see map of the errors for the HRG2000 in figure below) by the State Geodetic Authority. In Serbia, the global transformation parameters were determined based on 1217 referent points (the average 3D error is 0.4 meters).  
  
 
<br>  
 
<br>  
Line 109: Line 110:
 
|}
 
|}
  
The MGI / Balkans Zone 6 in the [http://trac.osgeo.org/proj/wiki/GenParms PROJ.4] should be set as:  
+
The MGI / Balkans Zone 6 in the [http://trac.osgeo.org/proj/wiki/GenParms PROJ.4] should be set as (thanks to Vedran Stojnović for the correct reference):  
 +
 
 
<pre>+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel  
 
<pre>+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel  
   +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824
+
   +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,5.541764
 
   +units=m
 
   +units=m
 
</pre>  
 
</pre>  
 +
 
where <tt>tmerc</tt> is the [http://www.remotesensing.org/geotiff/proj_list/ Transverse Mercator] projection system, <tt>ellps</tt> is the Bessel 1841 ellipsoid, and <tt>+towgs84</tt> carries the [http://trac.osgeo.org/proj/wiki/GenParms#towgs84-DatumtransformationtoWGS84 7 datum transformation parameters] (translation + rotation + scaling) to WGS84; also known as the [http://spatial-analyst.net/ILWIS/htm/ilwisapp/find_datum_trans_params_methodpage.htm ''Bursa Wolf''] method.  
 
where <tt>tmerc</tt> is the [http://www.remotesensing.org/geotiff/proj_list/ Transverse Mercator] projection system, <tt>ellps</tt> is the Bessel 1841 ellipsoid, and <tt>+towgs84</tt> carries the [http://trac.osgeo.org/proj/wiki/GenParms#towgs84-DatumtransformationtoWGS84 7 datum transformation parameters] (translation + rotation + scaling) to WGS84; also known as the [http://spatial-analyst.net/ILWIS/htm/ilwisapp/find_datum_trans_params_methodpage.htm ''Bursa Wolf''] method.  
  
Line 160: Line 163:
  
 
To validate these parameters, we can use a referent point (T249) for which we should have both WGS84 geographic coordinates and the local coordinates:  
 
To validate these parameters, we can use a referent point (T249) for which we should have both WGS84 geographic coordinates and the local coordinates:  
<pre>&gt; T249 &lt;- data.frame(Lat=as(char2dms("45d33'46.3998\"N"), "numeric"),
+
 
 +
<pre>
 +
> T249 <- data.frame(Lat=as(char2dms("45d33'46.3998\"N"), "numeric"),
 
+    Lon=as(char2dms("18d41'47.80491\"E"), "numeric"))
 
+    Lon=as(char2dms("18d41'47.80491\"E"), "numeric"))
&gt; coordinates(T249) &lt;- ~Lon+Lat
+
> coordinates(T249) <- ~Lon+Lat
&gt; proj4string(T249) &lt;- CRS("+proj=longlat +datum=WGS84")
+
> proj4string(T249) <- CRS("+proj=longlat +datum=WGS84")
</pre>  
+
</pre>
 +
 
 
To derive the local coordinates we use the <tt>spTransform</tt> method available in the [http://cran.r-project.org/web/packages/rgdal/ rgdal] package:  
 
To derive the local coordinates we use the <tt>spTransform</tt> method available in the [http://cran.r-project.org/web/packages/rgdal/ rgdal] package:  
<pre>&gt; T249.xy &lt;- spTransform(T249, CRS("+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999
+
 
 +
<pre>
 +
> T249.xy <- spTransform(T249, CRS("+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999
 
     +x_0=6500000 +y_0=0 +ellps=bessel  
 
     +x_0=6500000 +y_0=0 +ellps=bessel  
     +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824
+
     +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,5.541764
 
     +units=m"))
 
     +units=m"))
&gt; T249.xy
+
> T249.xy
SpatialPoints:
+
#!
          Lon      Lat
+
#!  SpatialPoints:
[1,] 6554781.6 5046739.2
+
#!          Lon      Lat
 +
#! [1,] 6554781.6 5046739.2
 +
#!
 
# should be: X=6554781.07,  Y=5046738.03
 
# should be: X=6554781.07,  Y=5046738.03
</pre>  
+
</pre>
  
 
The difference between the true and transformed local coordinates is:  
 
The difference between the true and transformed local coordinates is:  
  
<pre>&gt; spDistsN1(matrix(c(6554781.07, 5046738.03), ncol=2), coordinates(T249.xy))
+
<pre>  
 +
> spDistsN1(matrix(c(6554781.07, 5046738.03), ncol=2), coordinates(T249.xy))
 
[1] 1.2709
 
[1] 1.2709
 
</pre>  
 
</pre>  
Line 186: Line 197:
  
 
Next, we can check if the point ''falls'' where it should in Google maps:  
 
Next, we can check if the point ''falls'' where it should in Google maps:  
<pre>&gt; system(paste('"c:/Program Files/firefox/firefox.exe"', ' -url http://maps.google.com/maps?&amp;ll=',  
+
<pre>  
 +
> system(paste('"c:/Program Files/firefox/firefox.exe"', ' -url http://maps.google.com/maps?&amp;ll=',  
 
+    coordinates(T249)[2],',',coordinates(T249)[1],'&amp;t=h&amp;z=19&amp;q=',  
 
+    coordinates(T249)[2],',',coordinates(T249)[1],'&amp;t=h&amp;z=19&amp;q=',  
 
+    coordinates(T249)[2],',',coordinates(T249)[1], sep=""), wait=F)
 
+    coordinates(T249)[2],',',coordinates(T249)[1], sep=""), wait=F)
 +
#!
 
# http://maps.google.com/maps?&ll=45.56288883,18.696612475&t=h&z=19&q=45.56288883,18.696612475
 
# http://maps.google.com/maps?&ll=45.56288883,18.696612475&t=h&z=19&q=45.56288883,18.696612475
 
</pre>  
 
</pre>  
 +
 +
Which shows the correct location in the [http://maps.google.com/maps?&ll=45.56288883,18.696612475&t=h&z=19&q=45.56288883,18.696612475 Google Maps].
  
 
Or use the [http://cran.r-project.org/web/packages/RgoogleMaps/ RgoogleMaps] library to plot the maps in R:
 
Or use the [http://cran.r-project.org/web/packages/RgoogleMaps/ RgoogleMaps] library to plot the maps in R:
Line 210: Line 225:
 
We can also try reprojecting/transforming the coordinates in [[Software|ILWIS]]:
 
We can also try reprojecting/transforming the coordinates in [[Software|ILWIS]]:
 
   
 
   
<pre>?TRANSFORM(COORD(6554781.07,5046738.03,gk_6), LatLonWGS84)
+
<pre>
 +
?TRANSFORM(COORD(6554781.07,5046738.03,gk_6), LatLonWGS84)
 
= 45°33'46.4''N, 18°41'47.8''E
 
= 45°33'46.4''N, 18°41'47.8''E
 
</pre>  
 
</pre>  
Line 217: Line 233:
  
 
A simplified (but also less precise) solution is to use the 3 parameter datum shift (so called [http://spatial-analyst.net/ILWIS/htm/ilwisapp/find_datum_trans_params_methodpage.htm ''Molodensky'' method]):  
 
A simplified (but also less precise) solution is to use the 3 parameter datum shift (so called [http://spatial-analyst.net/ILWIS/htm/ilwisapp/find_datum_trans_params_methodpage.htm ''Molodensky'' method]):  
<pre>&gt; spTransform(T249, CRS("+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000
+
 
 +
<pre>
 +
> spTransform(T249, CRS("+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000
 
     +y_0=0 +ellps=bessel +towgs84=682,-199,480 +units=m"))
 
     +y_0=0 +ellps=bessel +towgs84=682,-199,480 +units=m"))
SpatialPoints:
+
#!
          Lon      Lat
+
#! SpatialPoints:
[1,] 6554780.2 5046738.3</pre>  
+
#!          Lon      Lat
 +
#! [1,] 6554780.2 5046738.3
 +
</pre>  
  
 
Note that the definition of the MGI / Balkans coordinate systems ([http://spatialreference.org/ref/epsg/31275/ EPSG:31275], [http://spatialreference.org/ref/epsg/31276/ EPSG:31276], [http://spatialreference.org/ref/epsg/31277/ EPSG:31277]) available via [http://spatialreference.org http://spatialreference.org] is incomplete because this database does not contain any definition of the geodetic datum. If the geodetic datum is missing, the difference from the true location will be even few hundreds of meters. In fact, if only one datum shift parameter is missing or is in error, transformation of coordinates will lead to frustrations. On the other hand, even if you correctly set-up the (global) coordinate system parameters, you will not be able to achieve precision better than 0.5 meters (see above).  
 
Note that the definition of the MGI / Balkans coordinate systems ([http://spatialreference.org/ref/epsg/31275/ EPSG:31275], [http://spatialreference.org/ref/epsg/31276/ EPSG:31276], [http://spatialreference.org/ref/epsg/31277/ EPSG:31277]) available via [http://spatialreference.org http://spatialreference.org] is incomplete because this database does not contain any definition of the geodetic datum. If the geodetic datum is missing, the difference from the true location will be even few hundreds of meters. In fact, if only one datum shift parameter is missing or is in error, transformation of coordinates will lead to frustrations. On the other hand, even if you correctly set-up the (global) coordinate system parameters, you will not be able to achieve precision better than 0.5 meters (see above).  
  
[[Image:T249_in_The_World_Coordinate_Converter.jpg|thumb|left|500px|Referent point T249 shown in [http://twcc.free.fr/?l=en&sc=WGS84&dc=UD1&wgs84=18.696612475%252C45.56288883 The World Coordinate Converter]. This uses the geodetic datum parameters from [http://spatialreference.org/ref/epsg/31276/ http://spatialreference.org].]]
+
{| border="0"
 
+
|- align="left"
 +
| [[Image:T249_in_The_World_Coordinate_Converter.jpg|thumb|left|600px|Referent point T249 shown in [http://twcc.free.fr/?l=en&sc=WGS84&dc=UD1&wgs84=18.696612475%252C45.56288883 The World Coordinate Converter]. This uses the geodetic datum parameters from [http://spatialreference.org/ref/epsg/31276/ http://spatialreference.org].]]
 +
|}
  
 
<br>
 
<br>
Line 235: Line 257:
  
 
*HRG2000 geoid: dX, dY, dZ, Rx, Ry, Rz, M_BF = 550.499, 164.116, 475.142, 5.80967, 2.07902, -11.62386, 0.99999445824;  
 
*HRG2000 geoid: dX, dY, dZ, Rx, Ry, Rz, M_BF = 550.499, 164.116, 475.142, 5.80967, 2.07902, -11.62386, 0.99999445824;  
 +
*HTRS96 geoid: dX, dY, dZ, Rx, Ry, Rz, M_BF = 550.5670, 164.6118, 474.1386, -5.976766, -2.099773, 11.495481, 0.999994552075;
 
*estimated by the [http://www.rgz.gov.rs Serbian Geodetic Authority]: dX, dY, dZ, Rx, Ry, Rz, M_BF = 574.027, 170.175, 401.545, 4.88786, -0.66524, -13.24673, 0.99999311067;  
 
*estimated by the [http://www.rgz.gov.rs Serbian Geodetic Authority]: dX, dY, dZ, Rx, Ry, Rz, M_BF = 574.027, 170.175, 401.545, 4.88786, -0.66524, -13.24673, 0.99999311067;  
 
*estimated by the [http://www.mapref.org/GeodeticReferenceSystemsAT.html Austrian Militärgeographisches Institut] (MGI) and valid for the whole Western Balkans region: dX, dY, dZ, Rx, Ry, Rz, M_BF = -678.059, -179.019, -585.545, -14.43, -0.42, -18.02, 0.999997490;
 
*estimated by the [http://www.mapref.org/GeodeticReferenceSystemsAT.html Austrian Militärgeographisches Institut] (MGI) and valid for the whole Western Balkans region: dX, dY, dZ, Rx, Ry, Rz, M_BF = -678.059, -179.019, -585.545, -14.43, -0.42, -18.02, 0.999997490;
Line 245: Line 268:
 
== Other coordinate systems of interest  ==
 
== Other coordinate systems of interest  ==
  
A coordinate system used in Croatia to display the whole country is known as the HTRS96/TM, centered at '''16°30''' (Basic, 2001). Its parameters are:
+
The most recent coordinate system used in Croatia to display the whole country is known as the [http://spatialreference.org/ref/epsg/3765/ '''Croatian Terrestrial Reference System'''] (HTRS96/TM) centered at '''16°30'''. Its parameters are ([http://hrcak.srce.hr/file/20088 Lapaine, 2007]):
 
   
 
   
<pre>+proj=tmerc +lat_0=0 +lon_0=16.5 +k=0.9997 +x_0=2500000 +y_0=0 +ellps=bessel
+
<pre>
   +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824
+
  +proj=tmerc +lat_0=0 +lon_0=16.5 +k=0.9999 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
  +units=m
+
</pre>
</pre>  
+
 
 +
which is aligned on the [http://en.wikipedia.org/wiki/European_Terrestrial_Reference_System_1989 European Terrestrial Reference System] (ETRS89). The most recent estimate of the global transformation parameters from the old geodetic datum (from year 1901) to HRTS96/TM for Croatian territory is:
 +
 
 +
<pre>
 +
   +towgs84=550.5670,164.6118,474.1386,-5.976766,-2.099773,11.495481,5.541764
 +
</pre>
  
 
Another coordinate system used by Military and similar organizations is the Universal Transverse Mercator system (Balkans falls in two zones: 33N and 34N). The parameters for Croatia are:  
 
Another coordinate system used by Military and similar organizations is the Universal Transverse Mercator system (Balkans falls in two zones: 33N and 34N). The parameters for Croatia are:  
  
<pre>+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
+
<pre>
 +
  +proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
 
</pre>  
 
</pre>  
  
Line 273: Line 302:
 
> HRbbox[1,2] <- ceiling(HRbbox[1,2]/pixsize)*pixsize-pixsize/2; HRbbox[2,2] <- ceiling(HRbbox[2,2]/pixsize)*pixsize-pixsize/2
 
> HRbbox[1,2] <- ceiling(HRbbox[1,2]/pixsize)*pixsize-pixsize/2; HRbbox[2,2] <- ceiling(HRbbox[2,2]/pixsize)*pixsize-pixsize/2
 
> HRbbox
 
> HRbbox
      min    max
+
#!
x 4572500 5062500
+
#!      min    max
y 2141500 2625500
+
#! x 4572500 5062500
 
+
#! y 2141500 2625500
 +
#!
 
# STEP 2: generate 1 km grid nodes;
 
# STEP 2: generate 1 km grid nodes;
 
> HR1kmgrid <- expand.grid(X=seq(HRbbox[1,1], HRbbox[1,2], by=pixsize),  
 
> HR1kmgrid <- expand.grid(X=seq(HRbbox[1,1], HRbbox[1,2], by=pixsize),  
Line 304: Line 334:
 
*Annoni, A. (Edt) 2003. [http://www.ec-gis.org/sdi/publist/pdfs/annoni2005eurgrids.pdf European Reference Grids - EUR Report 21494 EN]. Proceedings of the "European Reference Grids" workshop, Ispra, 27-29 October 2003.
 
*Annoni, A. (Edt) 2003. [http://www.ec-gis.org/sdi/publist/pdfs/annoni2005eurgrids.pdf European Reference Grids - EUR Report 21494 EN]. Proceedings of the "European Reference Grids" workshop, Ispra, 27-29 October 2003.
 
*Lapaine, M. 2001. Topografska kartografija u Hrvatskoj, Geografski horizont, 2, 1-53.
 
*Lapaine, M. 2001. Topografska kartografija u Hrvatskoj, Geografski horizont, 2, 1-53.
 +
*Lapaine, M. 2007. [http://hrcak.srce.hr/file/20088 New Official Map Projection of Croatia - HTRS96/TM]. Kartografija i geoinformacije, Vol. 6 No. izv. / spec.
  
 
----
 
----

Latest revision as of 00:23, 19 November 2018