[R-sig-Geo] LL into national grid system, discrepancy between spTransform and Geocalc
Roger Bivand
Roger.Bivand at nhh.no
Wed Jan 16 19:07:41 CET 2008
On Wed, 16 Jan 2008, ahimsa campos-arceiz wrote:
> Dear sig-geo listeRs,
>
> I've been using an old version of the program geocalc (The geographic
> calculator) to transform my field GPS data from the original WGS84 lat-long
> to the local National_grid_system (Kandawala, in Sri Lanka)
>
> I recently decided to update myself and try to use spTransform from the
> library rgdal. It took me some time to figure out how to define Kandawala's
> CRS. Once I that I (apparently) got it, I found that the resulting data show
> a discrepancy with the results from geocalc. This difference is a constant
> of 216 and 43 m for the x and y axis, respectively. I guess there might
> something in my configuration.
>
> I will appreciate if anybody can help me to identify the origin of this
> difference or suggest a better way to proceed. Below is my script, with a
> few points as example, plus some documentation.
Thank you for a very complete and helpful question!
The answer is very close, in the +towgs84= argument. The input data are in
WGS84, and your output data are also, but should not be - although the
ellipse is changed, the datum isn't. Your table below gives one set of
values:
my.kandawala1 <- paste(my.kandawala, "+towgs84=-98.3,787,85,0,0,0,0")
k.points <- spTransform(SP_wgs84.LL, CRS(my.kandawala1))
> k.points
SpatialPoints:
X Y
[1,] 246809.7 111593.8
[2,] 246820.7 111656.9
[3,] 246821.8 111713.3
[4,] 246837.3 111728.8
Coordinate Reference System (CRS) arguments: +proj=tmerc +ellps=evrst30
+lat_0=7.000480277777778 +lon_0=80.771711111111130 +k=0.9999238418
+x_0=200000.000 +y_0=200000.00 +towgs84=-98.3,787,85,0,0,0,0
> cbind(X2, Y2)
X2 Y2
[1,] 246811.2 111592.8
[2,] 246822.3 111655.8
[3,] 246823.3 111712.2
[4,] 246838.8 111727.7
which is very close. The rgdal package includes the current EPSG list of
coordinate reference systems, and doing:
> EPSG <- make_EPSG()
to load the list into a data frame, and searching with grep gives slightly
different values:
> EPSG[grep("Kandawala", EPSG$note),]
code note
163 4244 # Kandawala
prj4
163 +proj=longlat +a=6377276.345 +b=6356075.41314024
+towgs84=-97,787,86,0,0,0,0 +no_defs
so:
> my.kandawala1 <- paste(my.kandawala, "+towgs84=-97,787,86,0,0,0,0")
> k.points <- spTransform(SP_wgs84.LL, CRS(my.kandawala1))
> k.points
SpatialPoints:
X Y
[1,] 246811.0 111592.9
[2,] 246822.0 111655.9
[3,] 246823.1 111712.3
[4,] 246838.6 111727.8
Coordinate Reference System (CRS) arguments: +proj=tmerc +ellps=evrst30
+lat_0=7.000480277777778 +lon_0=80.771711111111130 +k=0.9999238418
+x_0=200000.000 +y_0=200000.00 +towgs84=-97,787,86,0,0,0,0
> cbind(X2, Y2)
X2 Y2
[1,] 246811.2 111592.8
[2,] 246822.3 111655.8
[3,] 246823.3 111712.2
[4,] 246838.8 111727.7
is only out by some tens of centimetres.
Finding the +towgs84= is typically hard, with the best source the APRS
Grids & Datums columns (http://www.asprs.org/resources/GRIDS/), but
unfortunately Sri Lanka isn't covered (yet). The column is excellently
written (also humorous), and all users of spatial data needing to use
legacy maps would benefit from browsing there.
Hope this helps!
Roger
>
> Thank you very much,
>
> Ahimsa
>
>
> #=== start script
>
> #=== data ====
> # GPS points recorded in the field (Sri Lanka) using lat-long WGS84
> ID <- c(1:4)
> X <- c(81.19672, 81.19682, 81.19683, 81.19697)
> Y <- c(6.20115, 6.20172, 6.20223, 6.20237)
> wgs84.LL <- data.frame(ID,X,Y)
>
> # the same data after transforming it using geocalc
> X2 <- c(246811.2398, 246822.2554, 246823.3166, 246838.7971)
> Y2 <- c(111592.7887, 111655.8285, 111712.2253, 111727.7191)
> kand <- data.frame(X2,Y2)
> # (below are the details of geocalc version and proj data)
>
>
> #=== transformation process ===
> # I create spatial data,
> library(rgdal)
> SP_wgs84.LL <- SpatialPoints(cbind(X,Y), proj4string=CRS("+proj=longlat
> +datum=WGS84"))
>
> # I define Sri Lankan National grid CRS,
> my.kandawala <- "+proj=tmerc +ellps=evrst30 +lat_0=7.000480277777778 +lon_0=
> 80.771711111111130 +k=0.9999238418 +x_0=200000.000 +y_0=200000.00"
>
> # and then I change the projection:
> k.points <- spTransform(SP_wgs84.LL, CRS(test))
> k.points
>
> #=== comparing the points from geocalc and spTransform:
> kand - as.data.frame(k.points)
> # there is a difference of ~216m and ~43m on the x and y axis, respectively
>
> #============== end script
>
> Documentation
>
> About Geocalc:
> I'm using Geocalc, the geographic calculator version 3.05 (1992-1994). The
> original points in degrees (dd.ddddd) and configured as "system=Geodetic",
> "Datum=186 WGS 1984" are converted into a user defined Kandawala Datum. The
> details of this datum are stored in a DAT file. These are the contents in
> such file:
>
> 1000 Zone_1 0 4 0 4 85 0 200000.000000000000000 200000.000000000000000
> 80.771711111111130 7.000480277777778 0.999923840000000
>
>
> # Kandawala in ArcGIS
> ArcGIS contains a prj file named kandawala. These are its contents:
>
> GEOGCS["GCS_Kandawala",DATUM["D_Kandawala",SPHEROID["Everest_1830",
> 6377276.345,300.8017]],PRIMEM["Greenwich",0],UNIT["Degree",
> 0.017453292519943295]]
>
> # I found more details about Kandawala on the web:
>
> http://www.mail-archive.com/mapinfo-l@csn.net/msg01850.html
>
> # kandawala data
>
> # FUGRO SURVEY - GEO Version 2.38.01
>
> # GEODESY 1
> # Datum : Kandawala
> # Spheroid : Everest 1830C
> # Semi Major Axis : 6377276.345 m
> # Inverse Flattening : 300.801700000
> # Projection : Transverse Mercator (UTM)
> # Latitude Origin : 7ø 00' 01.7290" N
> # Longitude Origin : 80ø 46' 18.1600" E
> # False Easting : 200000.000 m
> # False Northing : 200000.000 m
> # Central Scale Factor : 0.9999238418
>
>
> # PARAMETERS FOR CONVERSION FROM WGS 84
> # Kandawala WGS 84 Geod 1 to Geod 2
> # dX : 98.300 m 0.000 m -98.300 m
> # dY : -787.500 m 0.000 m 787.500 m
> # dZ : -85.000 m 0.000 m 85.000 m
> # rX : 0.0000 " 0.0000 " 0.0000 "
> # rY : 0.0000 " 0.0000 " 0.0000 "
> # rZ : 0.0000 " 0.0000 " 0.0000 "
> # dS : 1.00000 ppm 0.00000 ppm -1.00000 ppm
>
> # TRANSFORMATIONS
> # Station Name : PW-5
> # Datum : Kandawala : WGS 84
> # Projection : Transverse Mercator : Transverse Mercator
> # Latitude : 6ø 57' 09.3865" N : 6ø 57' 10.5800" N
> # Longitude : 79ø 50' 37.0926" E : 79ø 50' 44.7687" E
> # Height : 8.097 m : -93.704 m
> # Easting : 97459.221 m : 372492.852 m
> # Northing : 194807.060 m : 768701.916 m
> # PSF : 1.0000539617 : 0.9998012206
> # Convergence : - 0ø 06' 44.47" : - 0ø 08' 23.08"
>
>
>
>
>
>
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no
More information about the R-sig-Geo
mailing list