使用 ggplot/ggmap 在 R 中绘制凸包

Question

我有兴趣在 Ames 数据集中绘制邻域的多边形/凸包表示，类似于 O'Reilly 版本的 Tidy-modelling with R 中使用的方法。

同一段落的第 4 章似乎有两个不同的图像，具体取决于您使用的来源，第一个图像来自 here 的 O'Reilly 版本的书，第二个图像来自 tmwr.org

我个人更喜欢第一个，因为它更容易看到社区的重叠。 Tmwr.org 有一个 GitHub 存储库，但是那里没有捕获凸包代码。

从他们的仓库中，我可以看到他们从this网站下载了一个形状文件，并将其用于后台的道路测绘，然后覆盖多边形。

我正在尝试创建多边形并使用谷歌地图作为背景，但是我在添加谷歌地图图层方面没有取得太大成功。

在代码的底部，我可以绘制谷歌地图图层，但是当我引入

ames_sf

时，我收到以下错误：

Error in `geom_sf()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 4th layer.
Caused by error:
! object 'lon' not found

正如您从我用 ggplot Vs 生成的凸包中看到的那样。谷歌地图图像他们使用不同的坐标值，这可能是问题所在。我很困惑，因为第 46 行

ames_sf

是用

crs = 3857

定义的，然后在第 90 行

ggmap(map)

的底部被传递到

coord_sf

，其中

st_crs

也设置为 3857。

下面是我到目前为止拼凑在一起的代码（其中一些取自here），我正在寻找一些在谷歌地图背景中分层的帮助。任何可以分享的帮助将不胜感激。

# LOAD THE REQUIRED LIBRARIES
pacman::p_load(sf, ggplot2, dplyr, RColorBrewer, AmesHousing)

# ASSIGN A PALETTE OF 8 COLORS FROM THE 'DARK2' PALETTE TO A VARIABLE
col_vals <- brewer.pal(8, "Dark2")

# CREATE A NAMED VECTOR 'ames_cols' WITH NEIGHBORHOOD NAMES MAPPED TO SPECIFIC COLORS FROM 'col_vals'
# COPIED FROM THE TMWR.ORG REPO
ames_cols <- c(
  North_Ames = col_vals[3],
  College_Creek = col_vals[4],
  Old_Town = col_vals[8],
  Edwards = col_vals[5],
  Somerset = col_vals[8],
  Northridge_Heights = col_vals[4],
  Gilbert = col_vals[2],
  Sawyer = col_vals[7],
  Northwest_Ames = col_vals[6],
  Sawyer_West = col_vals[8],
  Mitchell = col_vals[3],
  Brookside = col_vals[1],
  Crawford = col_vals[4],
  Iowa_DOT_and_Rail_Road = col_vals[2],
  Timberland = col_vals[2],
  Northridge = col_vals[1],
  Stone_Brook = col_vals[5],
  South_and_West_of_Iowa_State_University = col_vals[3],
  Clear_Creek = col_vals[1],
  Meadow_Village = col_vals[1],
  Briardale = col_vals[7],
  Bloomington_Heights = col_vals[1],
  Veenker = col_vals[2],
  Northpark_Villa = col_vals[2],
  Blueste = col_vals[5],
  Greens = col_vals[6],
  Green_Hills = col_vals[8],
  Landmark = col_vals[3],
  Hayden_Lake = "red" # 'Hayden_Lake' SPECIFICALLY ASSIGNED THE COLOR RED
)

# LOAD AND PREPARE THE AMES HOUSING DATA
ames_data <- make_ames() %>%
  janitor::clean_names()

# CONVERT AMES DATA INTO AN SF (SIMPLE FEATURES) OBJECT FOR SPATIAL ANALYSIS AND
# TRANSFORM EPSG 3857 (Pseudo-Mercator, what Google uses)
ames_sf <- st_as_sf(ames_data, coords = c("longitude", "latitude"), crs = 3857)

# CREATE POLYGONS FOR EACH NEIGHBORHOOD BY AGGREGATING POINTS INTO A CONVEX HULL
ames_sf <- ames_sf %>%
  group_by(neighborhood) %>%
  summarize() %>%
  st_convex_hull()

# PLOT THE CONVEX HULLS USING GGPLOT
#ggplot(data = ames_sf) +
#  geom_sf(aes(fill = neighborhood), color = "white", size = 0.5) +
#  scale_fill_manual(values = ames_cols) +
#  theme_minimal() +
#  labs(title = 'Neighborhoods in Ames') +
#  theme(legend.position = "none")

map <- get_map("Ames, Iowa", maptype = "roadmap", zoom = "auto", scale = "auto", source = "google")

# Define a function to fix the bbox to be in EPSG:3857
ggmap_bbox <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector,
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(
    unlist(attr(map, "bb")),
    c("ymin", "xmin", "ymax", "xmax")
  )

  # Coonvert the bbox to an sf polygon, transform it to 3857,
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))

  # Overwrite the bbox of the ggmap object with the transformed coordinates
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}

# Use the function:
map <- ggmap_bbox(map)

ggmap(map) +
  coord_sf(crs = st_crs(3857)) + # force the ggplot2 map to be in 3857
  geom_sf(data = ames_sf, aes(fill = neighborhood), color = "white", size = 0.5, alpha = 0.5)

使用 ggmaps 我可以创建多边形并将其打印到坐标网格，但我很难添加到地图背景中。

Answer 1

当您比较问题中的这两个图时，您应该注意到的第一件事（或至少第二件事）是纬度/经度值的差异。即使不太确定纬度/经度的正常范围应该是多少，将美国的住房数据放在赤道（纬度 ~0）上也是一个非常明确的指标，表明最好后退几步。

这里的罪魁祸首是crs处理，使用

st_as_sf(ames_data, ... ,crs = 3857)

，您实际上并不是在转换坐标，而是将经度/纬度值定义为3857伪墨卡托的东距和北距，这里的

crs

应该与输入坐标的实际crs相匹配。 为了及早发现此类问题，请在每个步骤之后（如果需要）尽可能多地可视化地理空间数据。而且使用类似

mapview()

之类的东西要方便得多

以下示例使用

ggmap

和 OpenStreetMap 数据的 Carto Light 渲染，而不是

ggspatial

和 Google 磁贴。

library(AmesHousing)
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(ggplot2)
library(dplyr)
library(mapview)
library(ggspatial)

ames_sf <- make_ames() %>% 
  select(Neighborhood, Latitude, Longitude) %>% 
  janitor::clean_names() %>% 
  # crs in st_as_sf does not transform but defines the projection of source coordinates,
  # when dealing with lat/lon, it's almost always WGS84 aka EPSG:4326
  st_as_sf(coords = c("longitude", "latitude"), crs = "WGS84") %>% 
  group_by(neighborhood) %>% 
  summarise() %>% 
  st_convex_hull()
print(ames_sf, n = 3)
#> Simple feature collection with 28 features and 1 field
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: -93.69315 ymin: 41.9865 xmax: -93.57743 ymax: 42.06339
#> Geodetic CRS:  WGS 84
#> # A tibble: 28 × 2
#>   neighborhood                                                          geometry
#> * <fct>                                                            <POLYGON [°]>
#> 1 North_Ames    ((-93.60507 42.03456, -93.60795 42.03456, -93.62474 42.03461, -…
#> 2 College_Creek ((-93.6843 42.01339, -93.68792 42.01404, -93.69205 42.01621, -9…
#> 3 Old_Town      ((-93.62186 42.0258, -93.62374 42.02648, -93.62416 42.02997, -9…
#> # ℹ 25 more rows

# visualize with mapview to check that coordinates / location are OK
mapview(ames_sf, zcol = "neighborhood")

ggplot(ames_sf) +
  annotation_map_tile(type = "cartolight", zoomin = 0 ) +
  geom_sf(aes(fill = neighborhood), alpha = .2, show.legend = FALSE)
#> Zoom: 13
#> Fetching 16 missing tiles
#> ... |                                                                              |======================================================================| 100%
#> ...complete!

^{创建于 2023-11-06，使用 reprex v2.0.2}