将纬度/经度转换为状态平面坐标
我有一个经纬度数据集，我想使用EPSG 2790（http://spatialreference.org/ref/epsg/2790/）或ESRI 102672（ http://spatialreference.org/ref/esri/102672/）。
这肯定是以前被问过的; 我的代码基于这里的答案（rgdal R Package中的spTransform中的“Non Finite Transformation Detected”和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
当您为数据设置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
这里有一些更多的工作代码来阐明发生了什么：
# 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")])
}
我有一个问题转向另一个方向，并在GIS Stack Exchange上发现了这个回应，这可能对未来的求职者有所帮助。 根据您的坐标系是NAD83还是NAD83（HARN），空间参考epsg代码会有所不同。 如果使用错误的系统，则可能无法将该点转换为超出任何坐标平面限制的值。
https://gis.stackexchange.com/questions/64654/choosing-the-correct-value-for-proj4string-for-shape-file-reading-in-r-maptools
在lat + long与long + lat中读取确实会对输出产生影响 - 在我的情况下（从国家飞机到WGS84），我不得不写
coordinates(data) <- ~ long+lat
您可以通过绘制一个已知的参考点来确认它是否正确转换。
链接地址: http://www.djcxy.com/p/74117.html