如何将 geom_sf 生成的地图放在 ggmap 生成的栅格之上 [英] How to put a geom_sf produced map on top of a ggmap produced raster
本文介绍了如何将 geom_sf 生成的地图放在 ggmap 生成的栅格之上的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!
问题描述
我尝试了以下代码:
库(ggplot2)图书馆(ggmap)图书馆(SF)nc <- st_read(system.file("shape/nc.shp", package = "sf"))str(nc)类‘sf’和‘data.frame’:100 obs.共 15 个变量:美元面积:数量 0.114 0.061 0.143 0.07 0.153 0.097 0.062 0.091 0.118 0.124 ...$ 周长:num 1.44 1.23 1.63 2.97 2.21 ...$ CNTY_ : 编号 1825 1827 1828 1831 1832 ...$ CNTY_ID : 编号 1825 1827 1828 1831 1832 ...$ NAME :具有 100 个级别Alamance"、Alexander"、..的因子:5 3 86 27 66 46 15 37 93 85 ...$ FIPS :因子 w/100 级别 "37001","37003",..: 5 3 86 27 66 46 15 37 93 85 ...$ FIPSNO : 编号 37009 37005 37171 37053 37131 ...$ CRESS_ID : int 5 3 86 27 66 46 15 37 93 85 ...$ BIR74:编号 1091 487 3188 508 1421 ...$ SID74 : 数量 1 0 5 1 9 7 0 0 4 1 ...$ NWBIR74:编号 10 10 208 123 1066 ...$ BIR79:编号 1364 542 3616 830 1606 ...$ SID79 : 编号 0 3 6 2 3 5 2 2 2 5 ...$ NWBIR79:编号 19 12 260 145 1197 ...$几何:长度为100的sfc_MULTIPOLYGON;第一个列表元素:1 的列表..$ :1 的列表.. ..$ : num [1:27, 1:2] -81.5 -81.5 -81.6 -81.6 -81.7 .....- attr(*, "class")= chr "XY" "MULTIPOLYGON" "sfg"- attr(*, "sf_column")= chr "几何"- attr(*, "agr")= 因子 w/3 个级别 "constant","aggregate",..: NA NA NA NA NA NA NA NA NA .....- attr(*, "names")= chr "AREA" "PERIMETER" "CNTY_" "CNTY_ID" ...地图 <- get_map("北卡罗来纳州", maptype = "satellite", zoom = 6, source = "google")a <- unlist(attr(map,"bb")[1, ])bb <- st_bbox(nc)ggplot() +annotation_raster(map, xmin = a[2], xmax = a[4], ymin = a[1], ymax = a[3]) +xlim(c(bb[1], bb[3])) + ylim(c(bb[2], bb[4])) +geom_sf(data = nc, aes(fill = AREA))
两层没有正确重叠;我尝试使用 coord_sf()
更改投影,但没有成功.
有什么建议吗?谢谢
解决方案
我自己也在努力解决这个问题,并在 (v0.2.0) 于 2018 年 6 月 13 日创建.
I tried the following code:
library(ggplot2)
library(ggmap)
library(sf)
nc <- st_read(system.file("shape/nc.shp", package = "sf"))
str(nc)
Classes ‘sf’ and 'data.frame': 100 obs. of 15 variables:
$ AREA : num 0.114 0.061 0.143 0.07 0.153 0.097 0.062 0.091 0.118 0.124 ...
$ PERIMETER: num 1.44 1.23 1.63 2.97 2.21 ...
$ CNTY_ : num 1825 1827 1828 1831 1832 ...
$ CNTY_ID : num 1825 1827 1828 1831 1832 ...
$ NAME : Factor w/ 100 levels "Alamance","Alexander",..: 5 3 86 27 66 46 15 37 93 85 ...
$ FIPS : Factor w/ 100 levels "37001","37003",..: 5 3 86 27 66 46 15 37 93 85 ...
$ FIPSNO : num 37009 37005 37171 37053 37131 ...
$ CRESS_ID : int 5 3 86 27 66 46 15 37 93 85 ...
$ BIR74 : num 1091 487 3188 508 1421 ...
$ SID74 : num 1 0 5 1 9 7 0 0 4 1 ...
$ NWBIR74 : num 10 10 208 123 1066 ...
$ BIR79 : num 1364 542 3616 830 1606 ...
$ SID79 : num 0 3 6 2 3 5 2 2 2 5 ...
$ NWBIR79 : num 19 12 260 145 1197 ...
$ geometry :sfc_MULTIPOLYGON of length 100; first list element: List of 1
..$ :List of 1
.. ..$ : num [1:27, 1:2] -81.5 -81.5 -81.6 -81.6 -81.7 ...
..- attr(*, "class")= chr "XY" "MULTIPOLYGON" "sfg"
- attr(*, "sf_column")= chr "geometry"
- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
..- attr(*, "names")= chr "AREA" "PERIMETER" "CNTY_" "CNTY_ID" ...
map <- get_map("north carolina", maptype = "satellite", zoom = 6, source = "google")
a <- unlist(attr(map,"bb")[1, ])
bb <- st_bbox(nc)
ggplot() +
annotation_raster(map, xmin = a[2], xmax = a[4], ymin = a[1], ymax = a[3]) +
xlim(c(bb[1], bb[3])) + ylim(c(bb[2], bb[4])) +
geom_sf(data = nc, aes(fill = AREA))
The two layers do not overlap properly; I tried changing projection with coord_sf()
but I did not have success.
any suggestion? thanks
解决方案
I've struggled with this myself, and with the help of this post I've come up with a solution. The bounding box of the ggmap object is in WGS84 (EPSG:4326), but the actual raster is in EPSG:3857. You have to hack the bounding box of the ggmap object to be in the same CRS as the underlying data:
library(ggplot2)
library(ggmap)
library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.3.0, proj.4 5.1.0
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
# Transform nc to EPSG 3857 (Pseudo-Mercator, what Google uses)
nc_3857 <- st_transform(nc, 3857)
map <- get_map("north carolina", maptype = "satellite", zoom = 6, 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 = nc_3857, aes(fill = AREA), inherit.aes = FALSE)
Created on 2018-06-13 by the reprex package (v0.2.0).
这篇关于如何将 geom_sf 生成的地图放在 ggmap 生成的栅格之上的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!
查看全文