14

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 
Topdombili
  • 251
  • 1
  • 2
  • 4
  • 3
    @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 04 '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 04 '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 04 '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 04 '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 04 '13 at 22:45
  • how can we do that reverse direction? from X,Y to Latitude,Longitude converting – user142576 May 06 '19 at 11:49
  • This does not really answer the question. If you have a different question, you can ask it by clicking Ask Question. You can also add a bounty to draw more attention to this question once you have enough reputation. - From Review – TomazicM May 06 '19 at 12:12
  • This does not provide an answer to the question. Once you have sufficient reputation you will be able to comment on any post; instead, provide answers that don't require clarification from the asker. - From Review – whyzar May 06 '19 at 12:43

1 Answers1

18

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

whuber
  • 69,783
  • 15
  • 186
  • 281
  • 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 04 '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 04 '13 at 21:41
  • 1
    rgdal does recognize EPSG:2056, FWIW – mdsumner Jan 04 '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 04 '13 at 22:07