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

Question

我是R编程的新手，我想从两个文件中制作一个交互式地图，其中一个是可以从这里下载的.shp：https://www.ine.es/ss/Satellite?L=es_ES&c=页面&cid=1259952026632&p＝1259952026632&pagename=ProductosYServicios%2FPYSLayout(只需选择2021年及其下载(，其中有许多多边形。然后我有一个带有存储特征数据的csv(它包含2个LON和LAT字段(。

为了开始做这一切，我想为NCA字段中的每个不同值过滤.shp文件(例如：巴斯克国家的一张地图，马德里的一张，巴塞罗那的另一张…(。

所有这些都没有失去几何属性，因为如果我失去了它们，我就无法用图形表示它们(或者我可以，但我不知道，如果是，请告诉我，我会非常感激(。

他可能有一个共同的名字：

# 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

它返回给我的是：https://ibb.co/VDYdByq，行完全为空，我想获得的是相同的格式，但大约有1700行x18列，以及几何属性。

我的另一个问题是，当我将.shp文件读取为sf时，又添加了一列几何图形，其中是存储在列表中的坐标，如下所示：https://ibb.co/M1Fn8K5，我可以很容易地过滤它，但我不知道如何用图形表示(传单或地图视图…(，这样你就可以看到NCA="巴斯克国家"的多边形，有人能给我举个例子吗？我将非常感谢

Answer 1

好！我想我会做所有的工作流程！

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