将纬度/经度转换为州平面坐标答案

【问题标题】：Convert latitude/longitude to state plane coordinates将纬度/经度转换为州平面坐标
【发布时间】：2013-09-13 10:08:10
【问题描述】：

我有一个包含经纬度的数据集，我想使用 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

【问题讨论】：

标签： r gis latitude-longitude rgdal

【解决方案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")])
}

【讨论】：

谢谢，这是我找到的第一个有用的回复。

【解决方案2】：

当您为数据设置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

【讨论】：

对于任何通过谷歌搜索伊利诺伊州找到此内容的人，请不要混淆——此代码确实有效，即使切换了经度和纬度（即 41.20 是纬度，而不是经度在伊利诺伊州！）

【解决方案3】：

我在向另一个方向转换时遇到了问题，并在 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 确实会影响输出——在我的情况下（从 State Plane 到 WGS84）我不得不写

coordinates(data) <- ~ long+lat

您可以通过绘制已知参考点来确定它是否正确转换来确认这一点。

【讨论】：

【解决方案4】：

rgdal 中的 esri 代码 (102649) 对我不起作用，我必须从 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)

【讨论】：