r - 将纬度/经度转换为州平面坐标

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

下面是一些更有效的代码，可以阐明正在发生的事情:

# 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")])
}