Error in Converting XY Coordinates to Lat-Long Using r proj4 Library

Question

I am trying to write a script that will allow my colleagues to easily convert XY coordinates into lat and long. The idea here is to take a spreadsheet with only X and Y available and create a new sheet with latitude and longitude values appended, and to do this without working in a GIS interface at all.

Based on the code provided here, I wrote a script that I thought would do the trick. Most of my colleagues use R, so setting the working directory, reading in the right file, and running the script should be a snap for them. I am testing this with the NYPD's year-to-date complaint data, which already has lat and long but also comes with X and Y coordinates requiring the same projection system (FIPS 3104) as the actual data my colleagues will use:

####libraries####
library(proj4)
library(tidyverse)
library(leaflet)

####working directory####
wd <- ("my working directory")
setwd(wd)

####file location####
file <- ("/NYPD Open Data -- Has State Plane Coords/NYPD_Complaint_Map__no lat lon.csv")

####load data####
xy <- read.csv(paste0(wd,file))%>%
  filter(!is.na(X_COORD_CD),
         !is.na(Y_COORD_CD))


####convert points####
  #assign projection
proj4string <- "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
  #project point data
pj <- project(xy, proj4string, inverse=TRUE)
  #create latlon
latlon <- data.frame(xy, lat=pj$y, lon=pj$x)

####check points visually####
leaflet() %>%
  addProviderTiles("CartoDB.Positron", options = providerTileOptions(noWrap = TRUE)) %>%
  setView(lng = -73.9760663, lat = 40.7916269, zoom = 12) %>%
  addCircleMarkers(data = latlon,
                   lng = ~lon, lat = ~lat,
                   radius = 5.0, color="red",stroke=FALSE,fillOpacity = .8

####save output####
write.csv(latlon, file = "test.csv")

Unfortunately, I am not a savvy enough GIS-er to debug the problem I'm having. When I run the script using NYPD's year-to-date complaint data, all of the lat-longs end up being points in Antarctica. Per NYPD's documentation for this data, the X and Y coords should be in the NAD83 New York State Plane-Long Island projection (FIPS 3104), which is the projection string I am using. I have tried swapping X and Y around and setting inverse to false, on the theory that this was a similar error to reversing lat and long when loading data into ArcGIS,but it didn't solve the issue (although it did swap the values from the lat column to the lon column and vice-versa).

I have seen questions about similar errors with pyproj, but don't really understand how to apply the solutions in R (the proj4 project command doesn't have an equivalent argument to pyproj's preserve_units, which seemed to be the problem for pyproj users). I also found this discussion about a bug in the proj4 library that sounds related, but I lack the skill/experience to get any actionable solution from it.

Answer 1

Here is a very short sample of code that shows the proj4 library is functioning and your proj4string is correct:

library(proj4)
proj4string <- "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"

# Source data
xy <- data.frame(x=1001557, y=217404)

# Transformed data
pj <- project(xy, proj4string, inverse=TRUE)
latlon <- data.frame(lat=pj$y, lon=pj$x)
print(latlon)

This outputs:

       lat       lon
1 40.76338 -73.93752

This is in NY and matches the data.

Focus on your csv file reading and writing. (R is not my strong suit, so there might be some code errors in those I am missing.)

Answer 2

After posting this I realized that it cannot correctly convert the points when I'm not telling it which columns contain the points I need it to convert.

Here's a version that seems to work just fine:

####libraries####
library(proj4)
library(tidyverse)
library(leaflet)

####working directory####
wd <- ("MY WD")
setwd(wd)

####file location####
file <- ("/NYPD Open Data -- Has State Plane Coords/NYPD_Complaint_Map__no lat lon.csv")

####load data####
original <- read.csv(paste0(wd,file))%>%
  filter(!is.na(X_COORD_CD),
         !is.na(Y_COORD_CD))%>%
  rowid_to_column()

####convert points####
  #make 2D
xy <- original %>% select(X_COORD_CD,Y_COORD_CD,rowid)
  #assign projection
proj4string <- "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
  #project point data
pj <- project(xy, proj4string, inverse=TRUE)
  #create latlon
points <- data.frame(xy, lat=pj$y, lon=pj$x)
  #merge
final <-  merge(original, points, by.x = "rowid", by.y = "rowid")

####check points -- note their location makes no sense####
leaflet() %>%
  addProviderTiles("CartoDB.Positron", options = providerTileOptions(noWrap = TRUE)) %>%
  setView(lng = -73.9760663, lat = 40.7916269, zoom = 12) %>%
  addCircleMarkers(data = final,
                   lng = ~lon, lat = ~lat,
                   radius = 5.0, color="red",stroke=FALSE,fillOpacity = .8)

####save output####
write.csv(final, file = "test.csv")

I think as long as x and y are the first two columns of the dataframe you're converting, it all works just fine. I probably could do this without the rowID and merge, just by reordering the columns instead.

Answer 3

I'd usually load my data into an sp object and use rgdal to project:

library(sp)
library(rgdal)

dfr <- read.csv('~/Downloads/NYPD_Complaint_Data_Current__Year_To_Date_.csv')
dfr <- with(dfr, dfr[!is.na(X_COORD_CD) & !is.na(Y_COORD_CD), ]
xy <- dfr[, c('X_COORD_CD', 'Y_COORD_CD')]
prj <- CRS("+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs")

spdf <- SpatialPointsDataFrame(coordinates=xy, data=dfr, proj4string=prj)
spdf4326 <- spTransform(spdf, CRS('+init=epsg:4326'))

plot(spdf4326)

The data plots in New York.

This said it looks like Richard Morgan's answer may be even simpler, but if you need to associate data with your points then sp (or sf) is the way to go.

Regarding xy or yx, R almost always uses xy because it's conventional ordering for plotting data. I've heard geographers classically used yx (lat-lon) but I think xy is more standard today, at least in R.