ggmap上的反lay sf对象 - 未对齐，特定于区域

Question

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

为了参考，我用Google API

获得了AK.Base

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"))

使用下面引用的网站的Geojson文件，数据框看起来像这样。

# 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

您遇到的两个主要问题：

Answer 1

，没有重新注射，因此offset

您的Ak.base对象穿过反赛式的对象，这使得解决方案这样的解决方案this thistricky

• 修复有点令人费解，但我无法做出更简洁的方法。 首先，让我们加载所需的软件包并检查

  get_map()

  数据的范围： get_map()
• 指出值是WGS84/EPSG：4326坐标，Xmin值在WGS84/EPSG的边界范围之外：4326（-180，180，180，-90，90，90）。为了确保范围有效，您可以首先将AK.base转换为spatraster对象，并使用

  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

terra::crop()

Xmin现在在WGS84/EPSG的允许边界范围内：4326。但是，这些值应为EPSG：3857，即Google的默认CRS。要纠正这一点：

# 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], ext(r)[3],  ext(r)[4])

# Crop r
r <- crop(r, ext_crop)

# Get new extent of SpatRaster
r_ext <- ext(r)
setNames(c(r_ext[1], r_ext[2], r_ext[3], r_ext[4]), 
         c("xmin", "xmax", "ymin", "ymax"))
#       xmin       xmax       ymin       ymax 
# -179.99316 -123.83105   47.88752   73.58465 

如图所示，数据现在具有正确的CR和范围。

最终绘制您的数据。由于r是spatraster，因此您不能使用

# Convert SpatRaster extent coordinates to EPSG:3857
ext3857 <- ext(project(vect(ext(r), crs = crs(r)), "EPSG:3857"))

# Change (and not project) CRS of SpatRaster
crs(r) <- "EPSG:3857"

# Use ext3857 coordinates to modify extent values of SpatRaster
ext(r) <- c(as.integer(ext3857[c(1, 2, 3, 4)]))

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

，所以使用

ggmap()

：

tidyterra::geom_spatraster_rgb()