I have a spatial polygon object and a spatial points object. The latter was created from xy latlong data vectors (called latitude and longitude, respectively) in a dataframe, while the former was simply read into R directly using rgdal. My code is as follows:
boros <- readOGR(dsn = ".", "nybb")
rats <- read.csv("nycrats_missing_latlong_removed_4.2.14.csv", header = TRUE)
coordinates(rats) <- ~longitude + latitude
At this point neither spatial object is projected. If I project these objects as follows:
proj4string(boros) <- CRS("+proj=lcc")
proj4string(rats) <- CRS("+proj=lcc")
Both objects are now projected, and both will successfully map with the plot() function as follows:
plot(boros)
plot(rats)
However when I try to plot them together:
plot(boros)
plot(rats, add = TRUE)
I get the first plot only, without the rats object superimposed over boros. However, and this is the big problem, I get NO error message, so I have been unable to determine where the disconnection is between these two spatial objects being able to speak to each other. Both commands run smoothly without error or warning, yet I am left with just the single plot. And when I check the projections of each object with proj4string() I get the same projection returned for each object.
I have now spent many, many hours over several days trying various ways of creating two spatial objects whose CRS and projections match such that they can be mapped together using plot(). Incidentally, one approach I took was to create a shapefile in ArcGIS for the rats object, which worked fine to create the file. But I am still left with the same inability of the two spatial objects to work together in R. I have tried many different projections for both objects, spTransform on both objects, nothing seems to work.
I have also included a dropbox link with the 2 data files I have described above: https://www.dropbox.com/sh/x0urdo6guprnw8y/tQdfzSZ384
1 Answer 1
If you draw the axes (argument axes=TRUE in your plot statements), you can see the different coordinate systems:
library("rgdal")
boros <- readOGR(dsn=".", "nybb")
rats <- read.csv("nycrats_missing_latlong_removed_4.2.14.csv", header=TRUE)
coordinates(rats) <- ~longitude + latitude
op <- par(mfrow=c(1,2))
plot(rats, axes=TRUE)
plot(boros, axes=TRUE)
par(op)
plot with axes argument
The boros data set is using NAD 1983 State Plane New York Long Island coordinate reference system (according to @mkennedy's comment above). Using spTransform we obtain the following result:
crs <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
proj4string(rats) <- crs
proj4string(boros) <- "+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 +ellps=GRS80 +towgs84=0,0,0"
plot(spTransform(boros, CRS(crs)), axes=TRUE)
plot(rats, add=TRUE, col="#FF000050", pch=19, cex=0.3)
plot