带有 sf 地图和替代 crs 的背景图像

Background images with sf maps and alternative crs

是否可以使用带有 sf 地图的背景图像,并使用正在使用的坐标参考系统转换图像?例如,下面的第一个图使用 NASA night lights image 和 EPSG:4326。第二张图片带有 EPSG:3035。数据点似乎已在 CRS 上正确绘制,但是,背景图像未显示。

CRS EPSG 4326: 显示正确 CRS EPSG 3035: 无背景图片 使用的代码:

# Download NASA night lights image
download.file("https://www.nasa.gov/specials/blackmarble/2016/globalmaps/
              BlackMarble_2016_01deg.jpg", 
              destfile = "BlackMarble_2016_01deg.jpg", mode = "wb")

# Load picture and render
earth <- readJPEG("BlackMarble_2016_01deg.jpg", native = TRUE)
earth <- rasterGrob(earth, interpolate = TRUE)

# Select crs as desired
crs.inuse = 4326  # WGS 84
crs.inuse = 3035  # ETRS89 / LAEA Europe

# Plot geoname localities with population >100k with ggplot    
ggplot() +
  annotation_custom(earth, xmin = -180, xmax = 180, ymin = -90, ymax = 90) +
  geom_sf(data = popn.sf,
          aes(geometry=geometry),
          alpha = 0.4,
          size = 2,
          colour = "white") +
  theme(panel.background = element_rect(fill = "#05050f", colour = "#05050f"), 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.title = element_blank(), 
        axis.text = element_blank(), 
        axis.ticks.length = unit(0, "cm"),
        legend.position = "none") +
  coord_sf(crs = st_crs(crs.inuse))

jpg 文件导入为 RasterBrick 并将其投影到您想要的最终投影的解决方法。最终图是通过 plot 函数创建的。我还绘制了海岸线,以便可以看出投影是如何与 RasterBrick.

一起工作的

您可以找到有关此技术的更多解释 here

# Libraries----
library(sf)
library(raster)
library(rnaturalearth) #To get cities and coast to reprex

# Load shapefiles----
#popm.sf
popn.sf = ne_download(50, type = "populated_places", returnclass = "sf")
popn.sf_4326 = st_transform(popn.sf, 4326)
popn.sf_3035 = st_transform(popn.sf_4326, 3035)

#coasts
coast.sf = ne_download(10,
                       category = "physical",
                       type = "coastline",
                       returnclass = "sf")
coast.sf_4326 = st_transform(coast.sf, 4326)
coast.sf_3035 = st_transform(coast.sf, 3035)


# Background----
# Download NASA night lights image
download.file(
  "https://www.nasa.gov/specials/blackmarble/2016/globalmaps/BlackMarble_2016_01deg.jpg",
  destfile = "BlackMarble_2016_01deg.jpg",
  mode = "wb"
)

#To brick
earth <- brick("BlackMarble_2016_01deg.jpg")
raster_4326 <- earth
projection(raster_4326) <-
  CRS(st_crs(popn.sf_4326)[["proj4string"]])
extent(raster_4326) <-
  c(-180, 180,-89.99, 89.99) # Sligth offset to 180,180,-90,90 to avoid errors

#Project raster
proj3035 <- st_crs(popn.sf_3035)[["proj4string"]]
raster_3035 = projectRaster(raster_4326, crs = proj3035) # Some warnings, but still working

#Extra: to plots----
png(
  "proj4326.png",
  bg = "#05050f",
  height = dim(raster_4326)[1],
  width = dim(raster_4326)[2]
)
par(mar = c(0, 0, 0, 0))
plotRGB(raster_4326, bgalpha = 0)
plot(coast.sf_4326$geometry, col = "blue", add = T)
plot(
  popn.sf_4326$geometry,
  col = adjustcolor("white", alpha.f = .4),
  pch = 20,
  cex = 3,
  add = T
)
dev.off()

png(
  "proj3035.png",
  bg = "#05050f",
  height = dim(raster_3035)[1],
  width = dim(raster_3035)[2]
)
par(mar = c(0, 0, 0, 0))
plotRGB(raster_3035, bgalpha = 0)
plot(coast.sf_3035$geometry, col = "blue", add = T)
plot(
  popn.sf_3035$geometry,
  col = adjustcolor("white", alpha.f = .4),
  pch = 20,
  cex = 3,
  add = T
)
dev.off()

原文件

使用 CRS 4326

使用 CRS 3035