Преобразование координат широты/долготы в состояние координат

Question

У меня есть набор данных с широтой и долготой, который я бы хотел преобразовать в координаты плоскости состояния для Иллинойса Востока, используя EPSG 2790 (http://spatialreference.org/ref/epsg/2790/) или, возможно, ESRI 102672 (http://spatialreference.org/ref/esri/102672/) ,

Это, безусловно, было задано раньше; мой код основан на ответах здесь ("Non Finite Transformation Detected" in spTransform in rgdal R Package и http://r-sig-geo.2731867.n2.nabble.com/Converting-State-Plane-Coordinates-td5457204.html).

Но по какой-то причине я не могу заставить его работать:

library(rgdal) 
library(sp) 
data = data.frame(long=c(41.20,40.05), lat=c(-86.14,-88.15)) 
coordinates(data) <- ~ long + lat 
proj4string(data) <- CRS("+init=epsg:4326") # latitude/longitude 
data.proj <- spTransform(data, CRS("+init=epsg:2790")) # illinois east 

Дает:

non finite transformation detected: 
    long lat    
41.20 -86.14 Inf Inf 
Error in spTransform(data, CRS("+init=epsg:2790")) : failure in points 1 
In addition: Warning message: 
In spTransform(data, CRS("+init=epsg:2790")) : 
    2 projected point(s) not finite 

Answer 1

При установке coordinates для ваших данных, вы должны установить широту до того, как долготы.

Другими словами, изменение:

coordinates(data) <- ~ long + lat

в

coordinates(data) <- ~ lat+long

И он должен работать.

library(rgdal) 
library(sp) 
data = data.frame(long=c(41.20,40.05), lat=c(-86.14,-88.15)) 
coordinates(data) <- ~ lat+long 
proj4string(data) <- CRS("+init=epsg:4326") 
data.proj <- spTransform(data, CRS("+init=epsg:2790")) 
data.proj 

Дал мне этот выход:

SpatialPoints: 
      lat  long 
[1,] 483979.0 505572.6 
[2,] 315643.7 375568.0 
Coordinate Reference System (CRS) arguments: +init=epsg:2790 +proj=tmerc 
+lat_0=36.66666666666666 +lon_0=-88.33333333333333 +k=0.9999749999999999 +x_0=300000 
+y_0=0 +ellps=GRS80 +units=m +no_defs 

Answer 2

Вот несколько больше рабочего кода, который разъясняет, что происходит:

# convert a state-plane coordinate to lat/long 
data = data.frame(x=400000,y=0) 
coordinates(data) <- ~ x+y 
proj4string(data) <- CRS("+init=epsg:2804") 
latlong = data.frame(spTransform(data, CRS("+init=epsg:4326"))) 
setnames(latlong,c("long","lat")) 
latlong 

дает:

long  lat 
1 -77 37.66667 

и:

# convert a lat/long to state-plane 
data = latlong 
coordinates(data) <- ~ long+lat 
proj4string(data) <- CRS("+init=epsg:4326") 
xy = data.frame(spTransform(data, CRS("+init=epsg:2804"))) 
setnames(xy,c("y","x")) 
xy 

дает:

> xy 
     y    x 
1 4e+05 -2.690839e-08 

А вот функцию:

# this is for Maryland 
lat_long_to_xy = function(lat,long) { 
    library(rgdal) 
    library(sp) 
    data = data.frame(long=long, lat=lat) 
    coordinates(data) <- ~ long+lat 
    proj4string(data) <- CRS("+init=epsg:4326") 
    xy = data.frame(spTransform(data, CRS("+init=epsg:2804"))) 
    setnames(xy,c("y","x")) 
    return(xy[,c("x","y")]) 
} 

Answer 3

у меня была проблема преобразования в другом направлении и нашел ответ на ГИС Stack бирже, которая может оказаться полезной для будущего убежища. В зависимости от того, будет ли ваша система координат NAD83 или NAD83 (HARN), пространственный справочный epsg-код будет отличаться. Если вы используете неправильную систему, возможно, она не сможет преобразовать точку в значения за пределы любой координатной плоскости.

https://gis.stackexchange.com/questions/64654/choosing-the-correct-value-for-proj4string-for-shape-file-reading-in-r-maptools

Чтения в латах + длинный по сравнению с длиной + латы делает разницу в output- в моем случае (переход от State Plane к WGS84) я должен был написать

coordinates(data) <- ~ long+lat 

Вы можете подтвердите это, построив известную опорную точку, чтобы определить, правильно ли она была преобразована.

Answer 4

Руководство ESRI код (102649) в rgdal не работает для меня, я должен был вручную кодировать его со страницы proj4js, чтобы перейти от состояния самолета (0202 Arizona Central) в WGS84:

d<- data.frame(lon=XCord, lat=YCord) 
coordinates(d) <- c("lon", "lat") 
proj4string(d) <- CRS("+proj=tmerc +lat_0=31 +lon_0=-111.9166666666667 +k=0.9999 +x_0=213360 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs") 
CRS.new <- CRS("+init=epsg:4326") # WGS 84 
d.ch102649 <- spTransform(d, CRS.new)