18

Zillow has a set of shapefiles for different neighborhoods of major US cities. I wanted to check if certain buildings were present in certain neighborhoods using R:

library(rgeos)
library(sp)
library(rgdal)

df <- data.frame(Latitude =c(47.591351, 47.62212,47.595152),
                 Longitude = c(-122.332271,-122.353985,-122.331639),
                 names = c("Safeco Field", "Key Arena", "Century Link"))
coordinates(df) <- ~ Latitude + Longitude

wa.map <- readOGR("ZillowNeighborhoods-WA.shp", layer="ZillowNeighborhoods-WA")

sodo <- wa.map[wa.map$CITY == "Seattle"  & wa.map$NAME == "Industrial District", ]

I can plot without any issues

plot(sodo)
points(df$Latitude ~ df$Longitude, col = "red", cex = 1)

enter image description here

I match the proj4 string from the shapefile to my data.frame

CRSobj <- CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 ")
df@proj4string <- CRSobj

over(df, sodo)

This just gives me a bunch of NA values. I have tried this answer

spp <- SpatialPoints(df)
spp@proj4string <- CRSobj
over(spp, sodo)

but still get only NA values. Any ideas what else I should try?

19

The spatial data.frame is not correctly formed. This might work:

library(rgeos)
library(sp)
library(rgdal)

wa.map <- readOGR("ZillowNeighborhoods-WA.shp", layer="ZillowNeighborhoods-WA")

sodo <- wa.map[wa.map$CITY == "Seattle"  & wa.map$NAME == "Industrial District", ]

# Don't use df as name, it is an R function
# Better to set longitudes as the first column and latitudes as the second
dat <- data.frame(Longitude = c(-122.332271,-122.353985,-122.331639),
                  Latitude =c(47.591351, 47.62212,47.595152),
                  names = c("Safeco Field", "Key Arena", "Century Link"))
# Assignment modified according
coordinates(dat) <- ~ Longitude + Latitude
# Set the projection of the SpatialPointsDataFrame using the projection of the shapefile
proj4string(dat) <- proj4string(sodo)

over(dat, sodo)
#  STATE COUNTY    CITY                NAME REGIONID
#1    WA   King Seattle Industrial District   271892
#2  <NA>   <NA>    <NA>                <NA>       NA
#3  <NA>   <NA>    <NA>                <NA>       NA

over(sodo, dat)
#           names
#122 Safeco Field
7

I've just been doing the same thing. Pascal's answer almost covers it but you may need two extra steps as below.

#After you create your list of latlongs you must set the proj4string to longlat
proj4string(dat) <- CRS("+proj=longlat")

#Before you re-set the proj4string to the one from sodo you must actually convert #your coordinates to the new projection
dat <- spTransform(dat, proj4string(sodo))
  • It is not clear to me in which cases these additional steps are needed. For me the solution of the answer by user32309 was good enough. – djhurio Nov 27 '17 at 22:00
  • 3
    It depends on what format your data is in. If its in projection A and you just want to declare that you use proj4string and you should be good. But if its in projection B and you need to actuallly do a conversion to projection A you need to use spTransform. – John Curry Nov 28 '17 at 22:44

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.