Transform coordinates from Long Island Zone to Lat. Long format with R

Question

I have the following data set with coordinates in the State Plane Coordinate System from 1983. All coordinates are located in the Long Island Zone (3104).

dput(example)
structure(c(1008031L, 1000852L, 1001869L, 1005306L, 986887L, 
998031L, 1018703L, 1014319L, 1016186L, 1006977L, 1006891L, 1000883L, 
1001403L, 999812L, 1010077L, 1015918L, 984241L, 1013735L, 986848L, 
998243L, 1007312L, 1005663L, 992415L, 999771L, 1006787L, 987215L, 
990271L, 1015773L, 999342L, 1007245L, 1007098L, 996980L, 1006886L, 
999643L, 1008769L, 1016489L, 1004212L, 986848L, 1001512L, 1002584L, 
1001753L, 1004625L, 990725L, 1013435L, 1010795L, 1007509L, 1009419L, 
NA, 1009731L, 999007L, 999007L, 1000195L, 985863L, 990064L, 1008192L, 
1008306L, NA, 1003280L, 1006541L, 1001264L, 1003844L, 1008345L, 
987951L, 999104L, 1009013L, 998201L, 984182L, 1004940L, 1004513L, 
999659L, 1018204L, 1005918L, 1008158L, 999629L, 982208L, 1008008L, 
983985L, 1003591L, 992033L, 1012144L, 1008285L, 1004196L, 999937L, 
1007579L, 1001610L, 1013897L, 985504L, 1003588L, 1000088L, 1002230L, 
999304L, 1001393L, 997666L, 999148L, 997501L, 1004670L, 994699L, 
1005950L, 994821L, 998160L, 233036L, 228179L, 190702L, 186668L, 
173599L, 234924L, 241414L, 182198L, 178657L, 178140L, 242280L, 
236356L, 235184L, 238138L, 181374L, 245648L, 149582L, 211309L, 
212883L, 176387L, 243183L, 237170L, 149315L, 188471L, 242047L, 
215403L, 203844L, 240835L, 233575L, 234932L, 166665L, 174885L, 
193881L, 228852L, 244547L, 247336L, 178750L, 212883L, 232231L, 
248715L, 182080L, 242885L, 204176L, 251857L, 183147L, 245160L, 
235573L, NA, 243613L, 229814L, 229814L, 229856L, 212233L, 225331L, 
245037L, 245316L, NA, 229886L, 243541L, 232832L, 250988L, 235949L, 
220453L, 192913L, 242619L, 173610L, 150037L, 169914L, 180635L, 
229932L, 239783L, 190990L, 244973L, 243379L, 170319L, 246638L, 
205857L, 242274L, 215119L, 236944L, 245256L, 183865L, 238365L, 
183413L, 241367L, 238753L, 216029L, 249617L, 230093L, 176647L, 
227192L, 200533L, 177016L, 187285L, 170971L, 233870L, 176744L, 
179753L, 177866L, 227234L), .Dim = c(100L, 2L), .Dimnames = list(
    NULL, c("xcoord", "ycoord")))

I need to have the coordinates in lat./long. format or wgs84.

Can someone tell me how I do this? Thanks for your help.

Answer 1

Alternative to rgdal is proj4 package. It has transform function. I guessed the projection being used was http://www.spatialreference.org/ref/epsg/32118/

library(proj4)
dat <- structure(c(1008031L, 1000852L, 1001869L, 1005306L, 986887L, 
998031L, 1018703L, 1014319L, 1016186L, 1006977L, 1006891L, 1000883L, 
1001403L, 999812L, 1010077L, 1015918L, 984241L, 1013735L, 986848L, 
998243L, 1007312L, 1005663L, 992415L, 999771L, 1006787L, 987215L, 
990271L, 1015773L, 999342L, 1007245L, 1007098L, 996980L, 1006886L, 
999643L, 1008769L, 1016489L, 1004212L, 986848L, 1001512L, 1002584L, 
1001753L, 1004625L, 990725L, 1013435L, 1010795L, 1007509L, 1009419L, 
NA, 1009731L, 999007L, 999007L, 1000195L, 985863L, 990064L, 1008192L, 
1008306L, NA, 1003280L, 1006541L, 1001264L, 1003844L, 1008345L, 
987951L, 999104L, 1009013L, 998201L, 984182L, 1004940L, 1004513L, 
999659L, 1018204L, 1005918L, 1008158L, 999629L, 982208L, 1008008L, 
983985L, 1003591L, 992033L, 1012144L, 1008285L, 1004196L, 999937L, 
1007579L, 1001610L, 1013897L, 985504L, 1003588L, 1000088L, 1002230L, 
999304L, 1001393L, 997666L, 999148L, 997501L, 1004670L, 994699L, 
1005950L, 994821L, 998160L, 233036L, 228179L, 190702L, 186668L, 
173599L, 234924L, 241414L, 182198L, 178657L, 178140L, 242280L, 
236356L, 235184L, 238138L, 181374L, 245648L, 149582L, 211309L, 
212883L, 176387L, 243183L, 237170L, 149315L, 188471L, 242047L, 
215403L, 203844L, 240835L, 233575L, 234932L, 166665L, 174885L, 
193881L, 228852L, 244547L, 247336L, 178750L, 212883L, 232231L, 
248715L, 182080L, 242885L, 204176L, 251857L, 183147L, 245160L, 
235573L, NA, 243613L, 229814L, 229814L, 229856L, 212233L, 225331L, 
245037L, 245316L, NA, 229886L, 243541L, 232832L, 250988L, 235949L, 
220453L, 192913L, 242619L, 173610L, 150037L, 169914L, 180635L, 
229932L, 239783L, 190990L, 244973L, 243379L, 170319L, 246638L, 
205857L, 242274L, 215119L, 236944L, 245256L, 183865L, 238365L, 
183413L, 241367L, 238753L, 216029L, 249617L, 230093L, 176647L, 
227192L, 200533L, 177016L, 187285L, 170971L, 233870L, 176744L, 
179753L, 177866L, 227234L), .Dim = c(100L, 2L), .Dimnames = list(
    NULL, c("xcoord", "ycoord")))

lonlat <- project(dat, 
    '+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ', 
    inverse=T)

Answer 2

You can do this with rgdal and sp.

library(rgdal)
library(sp)
# I made a guess at the PROJ4 string that your data was in. This is using SRID: 2831
long.island.proj4<-CRS("+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
# Remove the missing values in your data, and convert to SpatialPoints
sp.points<-SpatialPoints(example[complete.cases(example),],proj4string=long.island.proj4)
# Project to lat/long
spTransform(sp.points,CRS('+proj=longlat'))