3

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. Any help would be most appreciated. I have also included a dropbox link with the 2 data files I have described above: https://www.dropbox.com/sh/x0urdo6guprnw8y/tQdfzSZ384

  • 1
    The Borough data is in NAD 1983 State Plane New York Long Island zone, using US survey feet. – mkennedy Apr 3 '14 at 18:52
4

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

  • Ah, OK so I was simultaneously failing to recognize that the 2 spatial objects had different coordinate systems (longlat vs. NAD 1983), and incorrectly specifying the projection. This has been tremendously helpful. Thank you very much for your help. Very much appreciated. – thegowda Apr 3 '14 at 22:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.