Conversion du système de coordonnées géographiques dans R

J'ai des points dans le système de coordonnées géographiques et je voulais les convertir en grille suisse (CH1903 +).

Exemples de données:

 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

Résultats respectés:

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

convert coordinate-system r Topdombili
la source

@Aaron, je suis le même gars. Je n'ai pas obtenu de réponse appropriée là-bas. Pouvez-vous m'aider? Je suis vraiment étonné de voir combien vous êtes pointilleux.

Topdombili

@Top Ce n'est pas du pickiness, c'est la politique StackExchange. La publication croisée crée toutes sortes d'incohérences et de problèmes. (Vous avez peut-être également remarqué que la publication sur le mauvais forum peut également obtenir des réponses moins qu'utiles.) Veuillez supprimer la version SO que vous avez publiée.

whuber

@Topdombili, je voulais juste souligner que selon la réponse de whuber, les valeurs d'entrée sont sur WGS84 et subissent une transformation de référence plus une projection vers la grille CH1903 + / LV95.

mkennedy

@mkennedy Merci pour cette observation. J'étais négligent en n'expliquant pas que j'avais supposé que les coordonnées d'origine (lat, lon) étaient WGS 84 (cette hypothèse est enterrée dans un commentaire dans le code). Si ce n'est pas le cas, modifiez la valeur de en proj4string(d)conséquence. Mon attention a été principalement attirée sur les faux paramètres d'abscisse et d'ordonnée, x0et y0, parce que certaines références populaires sur le Web (comme dans le premier commentaire du code) ont perdu leurs chiffres les plus significatifs, déplaçant ainsi tous les points de quelques milliers de kilomètres :-).

whuber

@whuber, aie re: les références! J'ai vu votre commentaire sur l'entrée définie sur WGS 84, mais je voulais m'assurer que Topdombili l'a vu.

mkennedy

Réponses:

Utilisez le package RGDAL . Il y a un problème de savoir quel CRS utiliser; RGDAL ne reconnaît pas le code EPSG. Vous devez fournir les paramètres explicitement, comme indiqué ici. (Apparemment, ce sont des approximations, mais elles sont censées être assez bonnes. Elles semblent être à environ 0,1 m des valeurs prévues.)

# 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)

Parcelles

        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
la source

Je voulais demander que le résultat de la précision de conversion soit inférieur alors qu'il n'y a pas de valeurs décimales alors que les coordonnées disponibles dans les résultats respectés sont au 10 degrés près, je voulais dire 2 chiffres décimaux.

Topdombili

Je ne comprends pas ton commentaire. Vous demandez-vous à quel point les coordonnées projetées sont précises lorsque les valeurs d'origine (lat, lon) sont spécifiées avec une précision limitée? Si tel est le cas, cette réponse peut vous être utile.

whuber

rgdal reconnaît EPSG: 2056, FWIW

mdsumner

@md Merci! J'avais trouvé une référence qui indique que c'est EPSG: 9814, mais RGDAL ne l'a pas reconnu.

whuber