将蒙版应用于R中的空间栅格？

Question

我想使用代表区域边界的 shapefile 来掩盖我的背景地图。为此，我使用 read_osm 将空间栅格读入 Rstudio

library(sp)
library(tmaptools)
HB_map <- spData::nz %>% 
  filter(Name=="Hawke's Bay") %>% 
  tmaptools::read_osm(type = "stamen-terrain")

然后我已经导入了我的 shapefile

libary(sf)
Regional_boundary <- sf::st_read("regional_boundary.shp")
sf::st_crs(Regional_boundary)= 2193 
Regional_boundary_sf_poly <- sf::st_transform(Regional_boundary, 27200) %>% 
  sf::st_cast(to="POLYGON")

我有许多 GIS 数据集，因此我重新投影了我的栅格，使其与我的 GIS 数据处于同一投影中（我不确定这一点是否正确）

test_map <- projectRaster(HB_map, crs ="+init=epsg:27200")

然后我检查数据预测是否一致

crs(Regional_boundary_sf_poly)

[1] "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,- 4.5993 +units=m +no_defs"

crs(test_map)

+init=epsg:27200 +proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +datum=nzgd49 +units=m +no_defs +ellps=intl +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993

现在我敷面膜：

Library(raster)
test_map_mask <- raster::mask(test,Regional_boundary_sf_poly,inverse=FALSE,updatevalue=NA, updadateNA=FALSE)

并检查结果

tmap::qtm(test_map_mask)

除了地图颜色不再像原来的“雄蕊地形”，而是不同深浅的橙色“之外，这一切似乎都有效。如何调整设置以使地图看起来像原来的但带有蒙版？

感谢您的帮助。 亲切的问候， 西蒙

Answer 1

感谢 dieghrnan 的回答，我现在可以获得一些示例数据，下面的内容对我来说运行良好

library(sf)
library(cartography)
library(raster) 
library(spData) 

nz <- spData::nz
HB <- cartography::getTiles(nz, type = "Stamen.Terrain")
HB_mask <- raster::mask(HB, nz)
plotRGB(HB_mask)

如果你想要一个投影：

test <- projectRaster(HB_mask, crs ="+init=epsg:27200")
# or perhaps
# test <- projectRaster(HB_mask, method="ngb", crs ="+init=epsg:27200")
plotRGB(test)

旧答案：

您在哪一步松开颜色？我的猜测是在 projectRaster？我无法让 rJava 工作，所以我无法运行它，但我会尝试

colortable(test_map_mask) <- colortable(HB_map)
tmap::qtm(test_map_mask)

Answer 2

新答案

尝试使用cartography::getTiles 获取瓷砖，它似乎有效。请注意，您应该安装了 sf v0.9-0 和 cartography v2.4-0，它们最近都在 CRAN 上进行了更新。在getTileshere 上提供完整的磁贴列表。

library(sf)
library(cartography)
library(dplyr)

nz <- spData::nz %>% st_transform(27200)

#Get tile
HB_map <- cartography::getTiles(nz, type = "Stamen.Terrain")

class(HB_map)
#> [1] "RasterBrick"
#> attr(,"package")
#> [1] "raster"

st_crs(HB_map)$proj4string
#> [1] "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993 +units=m +no_defs "
st_crs(nz)$proj4string
#> [1] "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993 +units=m +no_defs "

library(raster)
test_map_mask <-
  raster::mask(
    HB_map,
    nz,
  )

#tmap
tmap::qtm(test_map_mask)
#> Warning: replacing previous import 'sf::st_make_valid' by
#> 'lwgeom::st_make_valid' when loading 'tmap'
#> Numeric values of test_map_mask interpreted as RGB values with max.value = 255. Run tm_shape(test_map_mask) + tm_raster() to visualize the data.

^{由reprex package (v0.3.0) 于 2020 年 4 月 4 日创建}

旧答案

您可能想要检查cartography 中的getPngLayer 和pngLayer，基本上，您可以创建一个png 文件作为用sf 遮罩的光栅并绘制它。

如果你的问题是绘图，pngLayer 基本上是raster::plotRGB 的包装。

额外的好处是，您还可以使用 getTiles 和 cartography 获得 Stamen Terrain，但请确保安装最新版本（今天在 CRAN 上发布）。

代表：

library(sf)
library(cartography)
library(dplyr)

HB_map <- spData::nz %>% getTiles(type = "Stamen.Terrain")
class(HB_map)
#> [1] "RasterBrick"
#> attr(,"package")
#> [1] "raster"

Regional_boundary <- spData::nz  
st_crs(HB_map)$proj4string
#> [1] "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "
st_crs(Regional_boundary)$proj4string
#> [1] "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "

library(raster)

test_map_mask <-
  raster::mask(
    HB_map,
    Regional_boundary,
  )
#> Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
#> that

par(mar=c(0,0,0,0))
tilesLayer(test_map_mask)
plot(st_geometry(Regional_boundary), border = "red", lwd=2, add = TRUE)

^{由reprex package (v0.3.0) 于 2020 年 4 月 4 日创建}

小插图：https://dieghernan.github.io/cartographyvignette/

源码：https://github.com/riatelab/cartography/blob/master/R/tilesLayer.R