r - 带有 sf 地图和替代 crs 的背景图像
问题描述
是否可以将背景图像与 sf 地图一起使用,并且还可以使用正在使用的坐标参考系统对图像进行转换?例如,下面的第一个图是使用带有 EPSG:4326的NASA 夜灯图像。第二张图片是 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
为 aRasterBrick
并将其投影到所需的最终投影的解决方法。最终的图是通过plot
函数创建的。我还绘制了海岸线,因此可以看到投影如何与RasterBrick
.
您可以在此处找到有关此技术的更多说明。
# 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()
推荐阅读
- ansible - 如何在ansible中将今天的日期和时间作为参数传递
- c++ - 未找到 assimp-vc140-mt.dll
- ruby-on-rails - 如何使用 rack-cors 在 Access-Control-Expose-Headers 标头中声明 Content-Range
- visual-studio - 如何使用 PackageReferences 将 Microsoft.ReportingServices.ReportViewerControl.WebForms NuGet 包安装到非 SDK 样式项目?
- python - 从 Wikipedia 文章中删除参考书目
- go - 引擎盖下发生了什么,所以地图的这种并发使用是活泼的
- java - 在流终端操作中使用通配符
- pycharm - 失败时如何使用 pycharm 将代码推送到服务器?
- c++ - 如何将 QVector3D 与 QMatrix3x3 相乘?
- dsl - 如何在语法工具包中的 Lexer 识别元素之上包装复合 PsiElement