3
votes

I tried to use this instruction to convert a set of x-y coordinates in NAD83 State Plane Coordinates to regular Lan/Lat coordinates in degrees. I could reproduce for the example given in this post but it fails to give right answer to my set! Following is what I have tried and what I obtained [wrong answer].

library(rgdal)
nad83_coords <- data.frame(x=c(577430), y=c(2323270)) # My coordinates in NAD83
coordinates(nad83_coords) <- c('x', 'y')
proj4string(nad83_coords)=CRS("+init=esri:102272") # West Illinois
coordinates_deg <- spTransform(nad83_coords,CRS("+init=epsg:3436")) # West Illinois

> x             y
1 1894451.52045 7622263.00755

The correct solution must be:

lat= 38.2525     deg
lon=-90.07364722 deg

I have a code to do this job in MATLAB, but I am really interested to know if "rgdal" or any other library in R can do this. This way, my world would be more beautiful! Thank you for reading.

1
% INPUTS: % northing northing from State Plane reference point [USft] % easting easting from State Plane reference point [USft] % zone is the State Plane Coordinate System of 1983 zones % OUTPUTS: % lat (latitude) is positive north % lon (longitude) is positive WEST % lat & lon are in degrees (internal claculatiion in radians)Rotail
EPSG 3436 is also a NAD 83 projection. When I use EPSG 4326 (WGS84) I get -96.57821 x 42.86485. Those are off your correct solution, but I'm curious where those values are coming from.Badger
When I tried EPSG:3426 in only spTransform I get: 577430, 2323270. And when I try EPSG:3426 in both spTransform and proj4string, I get: 577430, 2323270 for x and y, respectively. These don't even look right!Rotail
This look weird that you use EPSG 4326 [somewhere in Africa] and get reasonable answer. I'm still working on reproducing your result.Rotail
3426 is also NAD 83, if you want lat long coordinates you need a WGS84 projection. You can use the ESRI: proj code you provided and the EPSG: 4326. The pin is at 0,0 (off the west coast of Africa), this is because WGS84 is a global projection, it is not specific for any single area.Badger

1 Answers

2
votes

In order to achieve Latitude and Longitude coordinates you need to go from you projected coordinate system (NAD 83) to a geographic coordinate system (WGS 84). In the case of your data you are using a projection in feet, so your West Illinois projection is correct. However, your mistake in the orignal post and highlighted below was that you spTransform from NAD83 to NAD83 giving you erroneous data. More information on the difference between projected and geographic coordinate systems can be found here. Instead you need to use a WGS84 projection in your transformation, as follows: ** as noted by Jim, the conversion from feet to meters is now incorporated.

library(rgdal)
   nad83_coords <- data.frame(x=c(577430), y=c(2323270)) # My coordinates in NAD83
nad83_coords <- nad83_coords *.3048 ## Feet to meters
    coordinates(nad83_coords) <- c('x', 'y')
    proj4string(nad83_coords)=CRS("+init=esri:102272") # West Illinois
    ## Erroneous code, NAD83 to NAD83
    ## coordinates_deg <- spTransform(nad83_coords,CRS("+init=epsg:3436"))
    ## Corrected with WGS84 to yield Lat/Long
    coordinates_deg <- spTransform(nad83_coords,CRS("+init=epsg:4326"))
    coordinates_deg
    SpatialPoints:
                x        y
    [1,] -96.57822 42.86484
    Coordinate Reference System (CRS) arguments: +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

In the event you did not want to convert we can use the following projection EPSG: 3531

library(rgdal)
   nad83_coords <- data.frame(x=c(577430), y=c(2323270)) # My coordinates in NAD83
    coordinates(nad83_coords) <- c('x', 'y')
    proj4string(nad83_coords)=CRS("+init=EPSG:3531") # West Illinois
    ## Erroneous code, NAD83 to NAD83
    ## coordinates_deg <- spTransform(nad83_coords,CRS("+init=epsg:3436"))
    ## Corrected with WGS84 to yield Lat/Long
    coordinates_deg <- spTransform(nad83_coords,CRS("+init=epsg:4326"))
    coordinates_deg
    SpatialPoints:
                x        y
    [1,] -96.57821 42.86485
    Coordinate Reference System (CRS) arguments: +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

I see what you're saying about where the coordinates should be, looking at the spatial reference site, your coordinates are within the window of that projection, but they aren't coming out as expected. Still digging.