r - st_centroid() 不在 Lambert-93 shapefile 上的多边形内
问题描述
我有这个问题。我有一个来自IGN 研究所的 shapefile 。我想根据一个特征聚合 shapefile,然后我想绘制生成的多边形的质心。
library(sf)
library(dplyr)
# load a shapefile and plot
sh = st_read("~/Downloads/CONTOURS-IRIS_2-1__SHP__FRA_2018-06-08/CONTOURS-IRIS/1_DONNEES_LIVRAISON_2018-06-00105/CONTOURS-IRIS_2-1_SHP_LAMB93_FXX-2017/CONTOURS-IRIS.shp",stringsAsFactors=FALSE)
# paris
p = sh %>%
mutate(INSEE_N = as.integer(INSEE_COM)-75100) %>%
filter((INSEE_N < 21) & (INSEE_N > 0))
# aggregate and plot centroid
p_arr = p %>%
group_by(INSEE_COM) %>%
summarise(n = n()) %>%
select(n)
plot(p_arr)
plot(st_geometry(st_centroid(p_arr)), pch = 3, col = 'red', add = TRUE)
我很惊讶这些质心不包含在它们各自的多边形中。如您所见,我没有触及 CRS 规范,我只是阅读了 shapefile。我在这里做错了什么?
编辑
顺便说一句,在 WGS84 中我根本看不到质心:
p_wgs = st_transform(p,4326)
p_arr_wgs = p_wgs %>%
group_by(INSEE_COM) %>%
summarise(n = n()) %>%
select(n)
plot(p_arr_wgs,axes=TRUE)
plot(st_geometry(st_centroid(p_arr_wgs)), pch = 3, col = 'red', add = TRUE)
Warning messages:
1: In st_centroid.sf(p_arr_wgs) :
st_centroid assumes attributes are constant over geometries of x
2: In st_centroid.sfc(st_geometry(x), of_largest_polygon = of_largest_polygon) :
st_centroid does not give correct centroids for longitude/latitude data
我看到了警告。但为什么不呢?数据看起来不太远?我希望至少能看到一些东西
> st_geometry(st_centroid(p_arr_wgs))
Geometry set for 20 features
geometry type: POINT
dimension: XY
bbox: xmin: 2.261954 ymin: 48.82841 xmax: 2.421378 ymax: 48.89257
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
First 5 geometries:
POINT (2.336403 48.86258)
POINT (2.342896 48.86828)
POINT (2.360002 48.86287)
POINT (2.357606 48.85435)
解决方案
不确定 ypour 代码有什么问题,但下面对我来说工作正常。
library(sf)
library(dplyr)
library(mapview)
paris <- st_read("./SO-answers/CONTOURS-IRIS_2-1__SHP__FRA_2018-06-08/CONTOURS-IRIS/1_DONNEES_LIVRAISON_2018-06-00105/CONTOURS-IRIS_2-1_SHP_LAMB93_FXX-2017/CONTOURS-IRIS.shp",
stringsAsFactors = FALSE)
p_arr <- paris %>%
mutate(INSEE_N = as.integer(INSEE_COM)-75100) %>%
filter((INSEE_N < 21) & (INSEE_N > 0)) %>%
group_by(INSEE_COM) %>%
summarise()
mapview( list( p_arr, st_centroid( p_arr ) ) )
推荐阅读
- r - 如何有条件地在 r 中拆分数据帧?
- nativescript - 如何在 iOS 平台上更新或刷新 RadListView 项?
- android - 如何通过 id 获取 Firebase 用户?
- python - foo=list、foo=[] 和 foo=list() 之间的区别?
- java - 如何使用 Java 编写余弦定理?
- intellij-idea - 我复制的项目出现在新项目旁边的括号中
- javascript - 更改图表的位置
- laravel - 让 Laravel “Hello World” 运行的最简单方法是什么?
- python-3.x - 如何获取小部件项的值以显示在 TextArea 中并将值转换为字典
- java - 错误:无法执行目标 org.apache.maven.plugins:maven-deploy-plugin:2.8.2:deploy 405,