Comprobar si los puntos están dentro de un polígono Shapefile

Question

Zillow tiene un conjunto de archivos shape para distintos barrios de las principales ciudades estadounidenses. Quería comprobar si ciertos edificios estaban presentes en determinados barrios utilizando 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", ]

Puedo trazar sin problemas

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

Coincido con el proj4 cadena del shapefile a mi data.frame

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

over(df, sodo)

Esto sólo me da un montón de NA valores. He probado lo siguiente responder

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

pero sigue obteniendo sólo NA valores. ¿Alguna idea sobre qué más debería probar?

Answer 1

El espacio data.frame no está correctamente formado. Esto podría funcionar:

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

Answer 2

Acabo de hacer lo mismo. La respuesta de Pascal casi lo cubre, pero es posible que necesite dos pasos adicionales como a continuación.

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

Answer 3

He utilizado un enfoque similar a la respuesta aceptada en este post, pero nunca estaba realmente satisfecho con él, así que busqué alternativas y encontré el sf biblioteca.

Y usando esta librería puedes escribir código como este:

library(sf)
# Shapefile from ABS: 
# https://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/1270.0.55.004July%202016?OpenDocument
map = read_sf("data/ABS/shapes/SUA_2016_AUST.shp")

pnts_sf <- st_as_sf(pnts, coords = c('y', 'x'), crs = st_crs(map))

pnts <- pnts_sf %>% mutate(
  intersection = as.integer(st_intersects(geometry, map))
  , area = if_else(is.na(intersection), '', map$SUA_NAME16[intersection])
) 

pnts

Salida:

         geometry intersection area    
*     <POINT [°]>        <int> <chr>   
1 (138.62 -34.92)           79 Adelaide
2 (138.58 -34.93)           79 Adelaide
3 (138.52 -34.95)           79 Adelaide
4 (152.71 -27.63)           60 Brisbane
5 (153.01 -27.57)           60 Brisbane
6  (150.73 -33.9)           31 Sydney  
7 (150.99 -33.92)           31 Sydney 

He publicado este código en otro post que era una pregunta similar, aquí: Identificar el polígono que contiene el punto con el paquete R sf