r - 在R中的shapefile中检查点向量属于哪个区域

Question

我有一个带有区域/多边形和特征的 shapefile。然后我有一个纬度/经度点的数据框。

我想在给定点所在的区域索引的小标题中添加一列。

一些可重现的代码：

### SETUP

library("tidyverse")

library("rgdal")
library("maptools")
library("rgeos")

library("ggmap")

unzipURL <- function(url, exdir){
  if(!dir.exists(exdir)){
    temp <- tempfile()
    download.file(url,temp)
    unzip(temp, exdir = exdir)
    unlink(temp)
  }
  return(exdir)
}

### DATA INGESTION

# Municipalities in Milan

municipi.mi.url <- "https://geoportale.comune.milano.it/Download/area_download/SIT/Municipi/Municipi_RDN2008.zip"

municipi.shp.path <- unzipURL(municipi.mi.url, exdir = file.path(".", "Municipi", criterion = is_git_root))

municipi.mi <- readOGR(dsn = municipi.shp.path, layer = "Municipi")
municipi.mi <- spTransform(municipi.mi,"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

# ARPA CONTROL STATIONS

centraline.mi.url <- "http://dati.comune.milano.it/dataset/d6960c75-0a02-4fda-a85f-3b1c4aa725d6/resource/59b3476e-9081-4601-a95a-a64925da3af5/download/qaria_stazione_shp.zip"

centraline.shp.path <- unzipURL(centraline.mi.url, exdir = file.path(".", "Centraline", criterion = is_git_root))

centraline.mi <- readOGR(dsn = centraline.shp.path, layer = "qaria_stazione")

### DATA PREP

municipi.mi@data$id = rownames(municipi.mi@data)
municipi.mi.points = fortify(municipi.mi, region="id")
municipi.mi.df = inner_join(municipi.mi.points, municipi.mi@data, by="id")

municipi.mi.df <- municipi.mi.df %>% mutate(MUNICIPIO = as.factor(MUNICIPIO))

# PLOT

ggplot(municipi.mi.df, mapping = aes(long, lat)) + 
  geom_polygon(mapping = aes(group=group, fill = MUNICIPIO)) +
  geom_path(mapping = aes(group=group), color="white") +
  geom_point(data=centraline.mi@data, mapping=aes(x=lon, y=lat), color = "white", cex=3) + 
  coord_map() +
  scale_fill_brewer(palette = "RdYlBu", type= "div") +
  theme(plot.background=element_blank(), 
        panel.background = element_blank(), 
        axis.title=element_blank(), 
        axis.ticks=element_blank(), 
        axis.text=element_blank()) +
  ggtitle("Municipalità di Milano",
          subtitle="Dislocazione centraline ARPA controllo della qualità dell'aria")

这给出了下面的情节

如前所述，我需要添加一列，centraline.mi@data其中region包含给定点所在的多边形的索引（在这个特定的示例中，我可以观察到，但这对于我拥有的所有层来说并不一定如此）

标签： rshapefileggmapsprgdal