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

convertfipslatitude longituderxy

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.

Best Answer

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.)

Related Question