Применение маски к пространственному растру в R?

Question

Я хотел бы замаскировать мою фоновую карту, используя шейп-файл, представляющий региональную границу. Для этого я прочитал пространственный растр в Rstudio, используя read_osm

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

Затем я импортировал свой шейп-файл

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

У меня есть несколько наборов данных ГИС, поэтому я перепроектирую мой растр, поэтому он находится в той же проекции, что и мои данные ГИС (я не уверен, верен ли этот бит)

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 + unit = m + no_defs "

crs(test_map)

+ init = epsg: 27200 + proj = nzmg + lat_0 = -41 + lon_0 = 173 + x_0 = 2510000 + y_0 = 6023150 + датум = nzgd49 + единиц = 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

НОВЫЙ ОТВЕТ

Попробуйте получить плитки с cartography::getTiles, похоже, работает. Обратите внимание, что у вас должны быть установлены sf v0.9-0 и cartography v2.4-0, оба недавно обновлены на CRAN . Полный список плиток доступен на getTiles здесь .

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.

^{Создано в 2020-04-04 представьте пакет (v0.3.0)}

СТАРЫЙ ОТВЕТ

Вы можете проверить getPngLayer и pngLayer из cartography, в основном с этим вы создайте png-файл как растр, замаскированный sf, и нанесите его на график.

Если ваша проблема заключается в построении графика, pngLayer в основном является оболочкой raster::plotRGB.

Extra ball Вы также можете получить Stamen Terrain, используя getTiles, также на cartography, но не забудьте установить последнюю версию (выпущенную сегодня на 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)

^{Создано в 2020-04-04 пакетом представитель (v0.3.0)}

Виньетка: https://dieghernan.github.io/cartographyvignette/

Исходный код: https://github.com/riatelab/cartography/blob/master/R/tilesLayer.R

Answer 2

Благодаря ответу 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? Я не могу заставить работать r Java, поэтому я не могу запустить это, но я бы попробовал

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