使用sf（在R中）进行多面几何分组

Question

我有一个自定义的shapefile。当我一次性绘制它时，它可以正常工作。但是我想按照某些变量进行分组，以绘制特定形状的区域。例如：

+-------------+--------+-------+
|   county    | region | sales |
+-------------+--------+-------+
| washoe      |      1 |     5 |
| carson city |      1 |    10 |
| clark       |      2 |    15 |
| harmon      |      2 |    20 |
+-------------+--------+-------+

如果我运行：

leaflet() %>% addTiles %>% addPolygons(data=df)

这将根据这四个位置绘制四个独立的多边形。

但假设我想基于区域创建多边形。因此，理论上，输出将只包含两个多边形。一个是蒙哥马利和史蒂文斯的组合形状，另一个是莫西斯和哈蒙县的组合形状，并且我还想总结销售总额。

x <- df %>% group_by(region) %>% summarise(totals=sum(ones))

我会收到以下错误信息：

Error in CPL_geos_union(st_geometry(x), by_feature) : 
  Evaluation error: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point -119.97497995771261 39.521600907266631 at -119.97497995771261 39.521600907266631

关于几何变量的信息：

df$geometry

Geometry set for 2106 features 
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: -115.8968 ymin: 35.00184 xmax: -114.0428 ymax: 36.85366
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

我需要将它转换成其他格式吗？

- user9302275

---

我无法使用您在评论中链接的数据，使用问题中的代码重复geos错误（请编辑您的问题并在其中链接数据）。因此，这可能只是特定版本的geos存在问题，或者这不是问题中“df”所需的数据。首先，“df”必须具有名为“ones”的列才能起作用。那到底是怎么回事？还有，这是https://gis.stackexchange.com/questions/345011/grouping-by-with-multipolygon-geometry-with-sf-using-r的副本吗？ - Spacedman

Answer 1

使用OP的数据进行编辑：

shp <- st_read("file/path/sample_shape.shp")

shp_union <- shp %>% 
  mutate(turf = as.character(turf),
         turf = ifelse(is.na(turf), "NA", turf)) %>% 
  group_by(turf) %>% 
  summarize(geometry = st_union(geometry),
            total_sales = sum(sales))

原始几何形状

ggplot() + geom_sf(data = shp, aes(fill = sales))

新几何图形

ggplot() + geom_sf(data = shp_union, aes(fill = total_sales))

- Eugene Chong

---

谢谢你的回答，但我的shapefile仍然无法工作。创建“group”组件的目的是什么？我需要这样做吗，还是只是为了你的示例而定制的？ - user9302275

那只是我加入到“sf”包内置的“nc”数据中的任意分组，因为我没有你的数据。你能否发布一小段你的数据？ - Eugene Chong

所以我上面使用的数据是虚假的。这里是真实数据的链接。原理相同，只不过是区域-草皮-销售额，而不是县-区域-销售额。https://github.com/gooponyagrinch/custom_shape - user9302275

哦，还有一件事。读取完文件后，请添加以下内容：testdata$areas <- 1:2106你会发现总共有2106行。这些是它最初读取的多边形，在你第一次运行时会在地图上看到。 - user9302275

它通过分组将几何图形合并。比如你有每个国家的多边形。如果你按大陆分组，然后执行 st_union()，你将得到6个多边形，每个有一个有人居住的大陆。 - Eugene Chong