admin管理员组文章数量:1316011
I would like to overlay an sf object (mapping geometry with anization's data) over the google map showing terrain features. This was successful when using other US states, but not Alaska. What do I need to do/change in order to overlay the map more precisely?
ak1<-ggmap(AK.base)+
geom_sf(data=ak_region_download,mapping=aes(geometry=geometry,fill=n,color=''),inherit.aes = F,alpha=.5)+
scale_fill_continuous(na.value=NA,low = "darkblue", high = "orange",
name = 'number',
breaks=c(1,5,10,20,30,40))+
scale_color_discrete('')+
theme_bw()
ak1
For reference I achieved AK.base with a google API
AK.base <- get_map(location = c(lon=-152,lat=63.5), zoom = 4, scale = "auto",
maptype = c("terrain"),
messaging = FALSE, urlonly = FALSE, filename = "ggmapTemp",
crop = FALSE, color = c("color"), source = c("google"))
And the dataframe looks something like this, using the geojson file from the website referenced below.
# ::alaska-borough-and-census-area-boundaries/about
#
# Download shapefile from above link, place ZIP in folder called DataRaw within
# project directory, unzip
ak_region_download <-
here("DataRaw/Alaska_Borough_and_Census_Area_Boundaries.shp") |>
sf::st_read()
ak_region_download$n <- 1:30
I would like to overlay an sf object (mapping geometry with anization's data) over the google map showing terrain features. This was successful when using other US states, but not Alaska. What do I need to do/change in order to overlay the map more precisely?
ak1<-ggmap(AK.base)+
geom_sf(data=ak_region_download,mapping=aes(geometry=geometry,fill=n,color=''),inherit.aes = F,alpha=.5)+
scale_fill_continuous(na.value=NA,low = "darkblue", high = "orange",
name = 'number',
breaks=c(1,5,10,20,30,40))+
scale_color_discrete('')+
theme_bw()
ak1
For reference I achieved AK.base with a google API
AK.base <- get_map(location = c(lon=-152,lat=63.5), zoom = 4, scale = "auto",
maptype = c("terrain"),
messaging = FALSE, urlonly = FALSE, filename = "ggmapTemp",
crop = FALSE, color = c("color"), source = c("google"))
And the dataframe looks something like this, using the geojson file from the website referenced below.
# https://gis.data.alaska.gov/datasets/DCCED::alaska-borough-and-census-area-boundaries/about
#
# Download shapefile from above link, place ZIP in folder called DataRaw within
# project directory, unzip
ak_region_download <-
here("DataRaw/Alaska_Borough_and_Census_Area_Boundaries.shp") |>
sf::st_read()
ak_region_download$n <- 1:30
Share
Improve this question
asked Jan 29 at 21:55
M3LbaM3Lba
1657 bronze badges
1 Answer
Reset to default 2There's two main issues you're encountering:
- the
get_map()
data from Google are in EPSG:3857, but the CRS they get loaded into R with is WGS84/EPSG:4326 without reprojecting, hence the offset - your AK.base object crosses the anti-meridian, which makes solutions like this one trickier to implement
The fix is a bit convoluted, but I could not work out a more succinct approach.
First, let's load the required packages and examine the extent of the get_map()
data:
library(sf)
library(terra)
library(ggmap)
library(tidyterra)
AK.base <- get_map(
location = c(lon = -152, lat = 63.5),
zoom = 4,
scale = "auto",
maptype = c("terrain"),
messaging = FALSE,
urlonly = FALSE,
filename = "ggmapTemp",
crop = FALSE,
color = c("color"),
source = c("google"))
# View extent of get_map() object
ak_bbox <- setNames(unlist(attr(AK.base, "bb")),
c("ymin", "xmin", "ymax", "xmax"))
ak_bbox[c("xmin", "xmax", "ymin", "ymax")]
# xmin xmax ymin ymax
# -180.08105 -123.83105 47.88752 73.58465
Note that the values are WGS84/EPSG:4326 coordinates and the the xmin value is outside the bounding extent of WGS84/EPSG:4326 (-180, 180, -90, 90). To ensure the extent is valid, you can first convert AK.base to a SpatRaster object and crop it using terra::crop()
:
# Convert get_map() object to SpatRaster
r <- rast(AK.base)
# Create extent object that does not cross anti-meridian
ext_crop <- ext(-180, ext(r)[2:4])
# Crop r
r <- crop(r, ext_crop)
# Get new extent of SpatRaster
r_ext <- ext(r)
setNames(c(r_ext[1:4]),
c("xmin", "xmax", "ymin", "ymax"))
# xmin xmax ymin ymax
# -179.99316 -123.83105 47.88752 73.58465
The xmin is now within the permitted bounding extent of WGS84/EPSG:4326. However, these values should be EPSG:3857, which is Google's default CRS. To correct this:
# Convert r extent coordinates to EPSG:3857
ext3857 <- ext(project(vect(ext(r), crs = crs(r)), "EPSG:3857"))
# Change (and not project) CRS of r
crs(r) <- "EPSG:3857"
# Use ext3857 coordinates to modify extent values of r
ext(r) <- ext3857
r
# class : SpatRaster
# dimensions : 1280, 1278, 3 (nrow, ncol, nlyr)
# resolution : 4891.97, 4891.97 (x, y)
# extent : -20036747, -13784809, 6088163, 12349884 (xmin, xmax, ymin, ymax)
# coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857)
# source(s) : memory
# colors RGB : 1, 2, 3
# names : red, green, blue
# min values : 12, 8, 8
# max values : 254, 254, 254
As shown, the data now have the correct CRS and extent.
Finally, plot your data. As r is a SpatRaster, you can't use ggmap()
so use tidyterra::geom_spatraster_rgb()
instead:
# Project ak_region_download to same CRS as SpatRaster
ak_region_download <- st_transform(ak_region_download, st_crs(r))
# Plot
ggplot() +
geom_spatraster_rgb(data = r, maxcell = Inf) +
geom_sf(data = ak_region_download,
aes(fill = n, color = NA),
inherit.aes = F,
alpha = .5) +
scale_fill_continuous(na.value = NA, low = "darkblue", high = "orange",
name = "number",
breaks = c(1,5,10,20,30,40)) +
coord_sf(xlim = ext(r)[c(1,2)],
expand = FALSE) +
scale_color_discrete("")+
theme_bw()
If you examine the result, you will see that the Canada/US border does not align with your sf data. I checked the Google data against:
ak_sf <- tigris::states() |>
filter(NAME == "Alaska") |>
st_transform(st_crs(r))
and the same issue occurred. I can't say definitively whether this a projection issue or an issue with the data, but if you want to avoid the misaligned borders, you can use ggmap::get_googlemap()
. It allows you to modify the default output using the style =
parameter. The example below sets the country borders to "off":
AK.base <- get_googlemap(center = c(lon = -152, lat = 63.5),
source = "google",
maptype = "terrain",
scale = 2,
zoom = 4,
style = c(feature = "administrative.country",
element = "geometry.stroke",
visibility = "off")
)
Applying the rest of the code results in:
本文标签: rOverlay sf object on ggmapnot alignedspecific to regionStack Overflow
版权声明:本文标题:r - Overlay sf object on ggmap - not aligned, specific to region - Stack Overflow 内容由网友自发贡献,该文观点仅代表作者本人, 转载请联系作者并注明出处:http://www.betaflare.com/web/1741994994a2409863.html, 本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容,一经查实,本站将立刻删除。
发表评论