在 maptools 和 rgdal 被弃用之前，我使用此代码导入 shapefile、质心 csv 和通勤计数数据，以使用下面的代码创建流程图（为了易读，我删除了 ggplot 代码）。我正在尝试使用最近的通勤数据重新创建此流程图，以比较两个时间段之间的变化。
我读到read_sf/st_read是用来替代readOGR函数的，所以我使用st_read来导入shapefile。
要将 shapefile 连接到质心点，我需要获取多边形 ID，这是此连接的公共字段。但是，我用来获取存储为行名称的多边形 ID 的行
soton@data$id = rownames(soton@data)soton@data 错误： 没有适用于
的方法应用于“sf”类对象@
这段代码与maptools和rgdal一起工作得很好，我使用的是之前完全相同的shapefile。我非常感谢您帮助调整我的代码以使用新软件包。谢谢你。
library(tidyverse)
library(plyr)
library(ggplot2)
library(maptools) -- replaced with library(sf)
library(rgdal) -- removed this and loaded sp package
setwd("your working directory")
soton=st_read(dsn="SouthamptonMSOA",layer="SouthamptonOA")
ODmat=read.csv("ODSouthampton.csv")
centroids=read.csv("CentroidSouthampton.csv")
popdat=read.csv("SouthamptonWZ.csv")
ODmat$frx=0
ODmat$fry=0
ODmat$tox=0
ODmat$toy=0
for (i in 1:dim(ODmat)[1]){
ODmat$frx[i]=centroids$East[which(centroids$Code==ODmat$Origin[i])]
ODmat$fry[i]=centroids$North[which(centroids$Code==ODmat$Origin[i])]
ODmat$tox[i]=centroids$East[which(centroids$Code==ODmat$Destination[i])]
ODmat$toy[i]=centroids$North[which(centroids$Code==ODmat$Destination[i])]}
soton@data$id = rownames(soton@data)
soton.points = fortify(soton, region="id")
soton.df = join(soton.points, soton@data, by="id")
ODmatsub=subset(ODmat,Commuters>50)
figplot = ggplot()+
# ...
让我解释一下所有步骤：
数据准备，因为我们无法访问您的 SouthamptonMSOA 和 csv 文件：
# b <- osmdata::getbb("Southampton")
# boundaries <- osmdata::opq(b, timeout = 25*60) |>
# osmdata::add_osm_feature(key = "boundary", value = "administrative") |>
# osmdata::osmdata_sf()
#
# boundaries$osm_multipolygons |>
# sf::st_write(dsn = "data/south.gpkg")
#
# boundaries$osm_multipolygons |>
# sf::st_centroid() |>
# sf::st_coordinates() |>
# write.csv(file = "data/centroids.csv")
您可以跳过以上部分，因为您手头有数据。
让我们在图层中阅读：
sh <- sf::st_read(dsn = "data/south.gpkg")
#> Reading layer `south' from data source `data/south.gpkg' using driver `GPKG'
#> Simple feature collection with 19 features and 129 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -1.957281 ymin: 50.70574 xmax: -0.7293437 ymax: 51.38392
#> Geodetic CRS: WGS 84
仅展示内容：
sh |>
subset(select = c("name", "admin_level"))
#> Simple feature collection with 19 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -1.957281 ymin: 50.70574 xmax: -0.7293437 ymax: 51.38392
#> Geodetic CRS: WGS 84
#> First 10 features:
#> name admin_level geom
#> 1 Test Valley 8 MULTIPOLYGON (((-1.308833 5...
#> 2 Eastleigh 8 MULTIPOLYGON (((-1.312912 5...
#> 3 Southampton 6 MULTIPOLYGON (((-1.354701 5...
#> 4 New Forest 8 MULTIPOLYGON (((-1.347534 5...
#> 5 Hampshire 6 MULTIPOLYGON (((-0.9320425 ...
#> 6 Bursledon 10 MULTIPOLYGON (((-1.309554 5...
#> 7 Hamble-le-Rice 10 MULTIPOLYGON (((-1.349241 5...
#> 8 Hound 10 MULTIPOLYGON (((-1.349241 5...
#> 9 Hedge End 10 MULTIPOLYGON (((-1.320232 5...
#> 10 West End 10 MULTIPOLYGON (((-1.354701 5...
让我们读取
csvc <- read.csv("data/centroids.csv")
c
#> X.1 X Y
#> 1 1 -1.510652 51.12843
#> 2 2 -1.327202 50.93302
#> 3 3 -1.398709 50.91700
#> 4 4 -1.624719 50.85839
#> 5 5 -1.287777 51.06381
#> 6 6 -1.315738 50.88791
#> 7 7 -1.323930 50.86144
#> 8 8 -1.347316 50.87588
#> 9 9 -1.303966 50.91370
#> 10 10 -1.329964 50.93579
#> 11 11 -1.325985 50.97084
#> 12 12 -1.413366 50.96088
#> 13 13 -1.468417 50.95084
#> 14 14 -1.497882 50.98122
#> 15 15 -1.412990 50.86973
#> 16 16 -1.448976 50.89033
#> 17 17 -1.494418 50.91691
#> 18 18 -1.495827 50.85305
#> 19 19 -1.359211 50.95967
现在，我们将
csvsfc <- c|>
sf::st_as_sf(coords = c("X", "Y"), crs = "EPSG:4326")
我们只使用数据集中的 2 个多边形，其中
admin_level == 6sh |>
subset(admin_level == "6", select = c("name", "geom"))
#> Simple feature collection with 2 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -1.957281 ymin: 50.70574 xmax: -0.7293437 ymax: 51.38392
#> Geodetic CRS: WGS 84
#> name geom
#> 3 Southampton MULTIPOLYGON (((-1.354701 5...
#> 5 Hampshire MULTIPOLYGON (((-0.9320425 ...
最后我们将他们聚集在一起：
sh |>
subset(admin_level == "6", select = c("name", "geom")) |>
sf::st_join(c)
#> Simple feature collection with 19 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: -1.957281 ymin: 50.70574 xmax: -0.7293437 ymax: 51.38392
#> Geodetic CRS: WGS 84
#> First 10 features:
#> name X.1 geom
#> 3 Southampton 3 MULTIPOLYGON (((-1.354701 5...
#> 5 Hampshire 1 MULTIPOLYGON (((-0.9320425 ...
#> 5.1 Hampshire 2 MULTIPOLYGON (((-0.9320425 ...
#> 5.2 Hampshire 4 MULTIPOLYGON (((-0.9320425 ...
#> 5.3 Hampshire 5 MULTIPOLYGON (((-0.9320425 ...
#> 5.4 Hampshire 6 MULTIPOLYGON (((-0.9320425 ...
#> 5.5 Hampshire 7 MULTIPOLYGON (((-0.9320425 ...
#> 5.6 Hampshire 8 MULTIPOLYGON (((-0.9320425 ...
#> 5.7 Hampshire 9 MULTIPOLYGON (((-0.9320425 ...
#> 5.8 Hampshire 10 MULTIPOLYGON (((-0.9320425 ...
创建于 2024-03-05，使用 reprex v2.1.0