Converting Latitude / Longitude Coordinates to Plane State

Question

I have a latitude and longitude dataset that I would like to convert to state plane coordinates for Illinois East using EPSG 2790 ( http://spatialreference.org/ref/epsg/2790/ ) or possibly ESRI 102672 ( http : //spatialreference.org/ref/esri/102672/ ).

It has definitely been asked before; my code is based on the answers here ( "Infinite transformation detected" in spTransform in the rgdal R package and http://r-sig-geo.2731867.n2.nabble.com/Converting-State-Plane-Coordinates-td5457204.html ).

But for some reason I can't get it to work:

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 

gives:

 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

When you set the coordinates for your data, you must set the latitude before longitude.

In other words, change:

coordinates(data) <- ~ long + lat

to

coordinates(data) <- ~ lat+long

And it should work.

 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 

Gave me this conclusion:

 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

Here is another working code that explains what happens:

 # 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 

gives:

  long lat 1 -77 37.66667 

and

 # 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 

gives:

 > xy yx 1 4e+05 -2.690839e-08 

And here is the function:

 # 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

I had a problem converting in the other direction and found this answer on the GIS Stack Exchange, which may be useful to future searchers. Depending on whether your coordinate system is NAD83 or NAD83 (HARN), the spatial reference epsg code will be different. If you use the wrong system, it may not be able to convert the point to values ​​outside of any coordinate plane.

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

Reading in lat + long versus long + lat makes a difference in output - in my case (switching from State Plane to WGS84) I had to write

 coordinates(data) <- ~ long+lat 

You can confirm this by building a known anchor point to determine if it was converted correctly.

Answer 4

The code esri (102649) in rgdal did not work for me, I had to manually encode it from the proj4js page in order to switch from the state plane (0202 Arizona Central) to 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)