如何在 R 中对 SpatialPolygonsDataFrame 类型的文件进行子采样而不丢失 shp 文件的属性

How to subsample a file of type SpatialPolygonsDataFrame in R without losing the properties of the shp file

我是 R 编程的新手,我想用两个文件制作交互式地图,一个是 .shp,您可以从这里下载:https://www.ine.es/ss/Satellite?L=es_ES&c=Page&cid=1259952026632&p=1259952026632&pagename=ProductosYServicios%2FPYSLayout(只是 select 2021 year and go 及其下载),其中有许多多边形。然后我有一个带有商店特征数据的 csv(它包含 2 个 LON 和 LAT 字段)。

要开始做这一切,我想为 NCA 字段中的每个不同值过滤 .shp 文件(例如:巴斯克地区的 1 张地图,马德里的另一张地图,巴塞罗那的另一张地图......)。

所有这一切都不会丢失几何属性,因为如果我丢失了它们,那么我将无法以图形方式表示它们(或者也许我可以但我不知道,如果是这样,请告诉我,我将非常感激) .

他将在 el siguiente codigo 上尝试:

# Load the libraries
pacman::p_load(leaflet, leaflet.extras, mapview, rworldxtra, rgdal,raster, sf, tidyverse, readr, ggthemes)

# Load the .shp file in spdf format.
myspdf = readOGR(getwd(), layer = "SECC_CE_20210101")

#Filter
PV = myspdf %>% filter(NCA == "País Vasco") # Dont work
PV2 = myspdf[myspdf$NCA == "País Vasco"] # Dont work

当我加载 shp 文件并将其保存在变量 myspdf 中时,我可以想象这样的事情:https://ibb.co/mywDd6p

如果我执行 myspdf@data,我会在其中访问数据(我要过滤的 NCA 字段在哪里)

所以当我尝试这样过滤时:

PV = myspdf %>% filter(NCA == "País Vasco") # Dont work
PV2 = myspdf[myspdf$NCA == "País Vasco"] # Dont work

它 returns 这对我来说是这样的:https://ibb.co/VDYdByq,行完全是空的,我想获得的是相同的格式,但大约有 1700 行 x 18 列,并且几何属性也是如此。

我的另一个问题是,当我将 .shp 文件读取为 sf 时,又添加了一列几何图形,里面是存储在列表中的坐标,就像这样:https://ibb.co/M1Fn8K5,我可以很容易地对其进行过滤,但我不知道如何以图形方式(传单或地图视图...)表示它,以便您可以看到 NCA = 'Basque Country' 的多边形,有人可以给我举个例子吗?我将不胜感激

好的!我想我会完成所有工作流程!

library(sf)
library(tmap)
library(mapview)

# lets get some shops
shop <- data.frame(X = c(-4.758628, -4.758244, -4.756829, -4.759394, -4.753698,
                         -4.735330, -4.864548, -4.863816, -4.784694, -4.738924),
                   Y = c(43.42144, 43.42244, 43.42063, 43.42170, 43.41899,
                         43.41181, 43.42327, 43.42370, 43.42422, 43.40150),
                   name = LETTERS[1:10])
# Here I save them 
write.csv(shop, "shop.csv")

# because I want to show you how to import
shop <- read.csv("shop.csv")
# and convert to en sf object
shop_sf <- sf::st_as_sf(shop, coords = c("X", "Y"))
# and add a CRS
shop_sf  <- sf::st_set_crs(shop_sf, 4326)

# now I have downloaded data from your link 
# I import it in R 
spain_seccionado <- sf::st_read("España_Seccionado2021/SECC_CE_20210101.shp")
# Rq CRS is ETRS89 / UTM 30, will need to transform that

# here I am just exploring a bit the data set
names(spain_seccionado)
unique(spain_seccionado$NCA)

# I just keep Asturias, You have plenty of different way of doing that 
# this is what you tried to do here: PV = myspdf %>% filter(NCA == "País Vasco")
# but on an sp object not an sf one
Asturias <- spain_seccionado[spain_seccionado$NCA == "Principado de Asturias",]
asturias_4326 <- sf::st_transform(Asturias, 4326)
# Now both data set are in the same CRS
# a quick plot just to see if everything is correct 
plot(asturias_4326$geometry)
plot(shop_sf, col = "red", add = TRUE, pch = 5)


# An interactive map quick and dirty you will need to improve it ! 

tmap_mode("view")

llanes_shop <- tmap::tm_shape(asturias_4326) +
                    tmap::tm_borders() +
                    tmap::tm_shape(shop_sf) +
                    tmap::tm_symbols(shape = 24) +
                    tmap::tm_layout()

llanes_shop