Comprobando si los puntos caen dentro del archivo de forma del polígono

16

Zillow tiene un conjunto de shapefiles para diferentes vecindarios de las principales ciudades de los EE. UU. Quería verificar si ciertos edificios estaban presentes en ciertos vecindarios usando 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)

Coincidoconlacadenaproj4delshapefileconmidata.frame

CRSobj<-CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 ")
[email protected] <- CRSobj

over(df, sodo)

Esto solo me da un montón de valores NA . He probado esto answer / a>

spp <- SpatialPoints(df)
[email protected] <- CRSobj
over(spp, sodo)

pero aún así solo obtienes los valores NA . ¿Alguna idea de qué otra cosa debería probar?

    
pregunta Stedy 09.02.2015 - 03:41

2 respuestas

17

El data.frame espacial no se forma correctamente. 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
    
respondido por el user32309 09.02.2015 - 04:26
6

He estado haciendo lo mismo. La respuesta de Pascal casi lo cubre, pero es posible que necesite dos pasos adicionales como se muestra 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))
    
respondido por el John Curry 15.01.2016 - 10:58

Lea otras preguntas en las etiquetas