13

I have points in geographic coordinate system and I wanted to convert them to swiss grid (CH1903+).

Sample data:

 id       lon      lat
  2 7.173500 45.86880
  3 7.172540 45.86887
  4 7.171636 45.86924
  5 7.180180 45.87158
  6 7.178070 45.87014
  7 7.177229 45.86923  
  8 7.175240 45.86808
  9 7.181409 45.87177
  10 7.179299 45.87020

Respected Results:

id       E              N      
2     2579408.2431  1079721.1499
3     2579333.7158  1079729.1852
4     2579263.6502  1079770.1125
5     2579928.0358  1080028.4605
6     2579763.6471  1079868.9218
7     2579698.0756  1079767.9762
8     2579543.1019  1079640.6512
9     2580023.6226  1080049.2672
10    2579859.1889  1079875.2740 
  • 2
    @Aaron, I am the same guy. I have not get proper answer there. Can u help me? I am really amazed how much you are picky. – Topdombili Jan 4 '13 at 21:17
  • 1
    @Top This is not pickiness, it's StackExchange policy. Cross-posting creates all kinds of inconsistencies and problems. (You may also have noticed that posting in the wrong forum can get some less than useful answers, too.) Please delete the SO version you posted. – whuber Jan 4 '13 at 21:19
  • @Topdombili, I just wanted to point out that according to whuber's answer, the input values are on WGS84 and are undergoing a datum transformation plus a projection to the CH1903+ / LV95 grid. – mkennedy Jan 4 '13 at 21:41
  • @mkennedy Thank you for that observation. I was remiss in not explaining that I have assumed the original (lat, lon) coordinates are WGS 84 (that assumption is buried in a comment in the code). If they are not, change the value of proj4string(d) accordingly. My attention was primarily drawn to the false easting and northing parameters, x0 and y0, because some popular references on the Web (such as in the first comment in the code) have dropped their most significant digits, thereby displacing all the points by a few thousand kilometers :-). – whuber Jan 4 '13 at 22:20
  • 1
    @whuber, ouch re: the references! I did see your comment about the input set to WGS 84, but wanted to make sure Topdombili saw it. – mkennedy Jan 4 '13 at 22:45
17

Use the RGDAL package. There's an issue of which CRS to use; RGDAL does not recognize the EPSG code. You have to provide the parameters explicitly, as shown here. (Apparently these are approximations, but they're supposed to be pretty good. They seem to be within about 0.1 m of the intended values.)

# References:
# http://lists.maptools.org/pipermail/proj/2001-September/000248.html (has typos)
# http://www.remotesensing.org/geotiff/proj_list/swiss_oblique_cylindrical.html
#
# Input coordinates.
#
x <- c(7.173500, 7.172540, 7.171636, 7.180180, 7.178070, 7.177229, 7.175240, 7.181409, 7.179299)
y <- c(45.86880, 45.86887, 45.86924, 45.87158, 45.87014, 45.86923, 45.86808, 45.87177, 45.87020)
#
# Define the coordinate systems.
#
library(rgdal)
d <- data.frame(lon=x, lat=y)
coordinates(d) <- c("lon", "lat")
proj4string(d) <- CRS("+init=epsg:4326") # WGS 84
CRS.new <- CRS("+proj=somerc +lat_0=46.9524056 +lon_0=7.43958333 +ellps=bessel +x_0=2600000 +y_0=1200000 +towgs84=674.374,15.056,405.346 +units=m +k_0=1 +no_defs")
# (@mdsumner points out that
#    CRS.new <- CRS("+init=epsg:2056")
# will work, and indeed it does. See http://spatialreference.org/ref/epsg/2056/proj4/.)
d.ch1903 <- spTransform(d, CRS.new)
#
# Plot the results.
#
par(mfrow=c(1,3))
plot.default(x,y, main="Raw data", cex.axis=.95)
plot(d, axes=TRUE, main="Original lat-lon", cex.axis=.95)
plot(d.ch1903, axes=TRUE, main="Projected", cex.axis=.95)
unclass(d.ch1903)

Plots

        lon        lat  

[1,] 2579408.24 1079721.15
[2,] 2579333.69 1079729.18
[3,] 2579263.65 1079770.55
[4,] 2579928.04 1080028.46
[5,] 2579763.65 1079868.92
[6,] 2579698.00 1079767.97
[7,] 2579543.10 1079640.65
[8,] 2580023.55 1080049.26
[9,] 2579859.11 1079875.27

  • I wanted to ask the result of conversion accuracy might be less while there are without decimal values while the available coordinates in respected results are in nearest 10 degree, I meant 2 digits decimal. – Topdombili Jan 4 '13 at 21:39
  • 1
    I don't understand your comment. Are you asking about how precise the projected coordinates are when the original (lat, lon) values are specified with limited precision? If so, you may find this answer to be useful. – whuber Jan 4 '13 at 21:41
  • 1
    rgdal does recognize EPSG:2056, FWIW – mdsumner Jan 4 '13 at 22:01
  • @md Thank you! I had found a reference that states this is EPSG:9814, but RGDAL did not recognize it. – whuber Jan 4 '13 at 22:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.