Applying a mask to a spatial raster in R?

Question

I would like to mask my background map using a shapefile representing a regional boundary. To do this I have read a spatial raster into Rstudio using read_osm

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

I have then imported my 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")

I have a number of GIS data sets so I re-project my raster so it is in the same projection as my GIS data (i'm not sure if this bit is right)

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

I then check the data projections are consistent

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

Now I apply my mask:

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

And check the results

tmap::qtm(test_map_mask)

It all seems to work except the map colours no longer look like the original "stamen-terrain" but are instead different shades of orange". How do adjust the settings to make the map look like the original but with the mask?

Thanks for you help. Kind regards, Simon

Answer 1

NEW ANSWER

Try to get the tiles with cartography::getTiles, it seems to work. Note that you should have sf v0.9-0 and cartography v2.4-0 installed, both recently updated on CRAN. A complete list of tiles serves available on 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.

^{Created on 2020-04-04 by the reprex package (v0.3.0)}

OLD ANSWER

You may want to check getPngLayer and pngLayer from cartography, basically with this you create a png file as a raster masked with a sf and plot it.

If your problem is the plotting, pngLayer Is basically a wrapper of raster::plotRGB.

Extra ball is that you can get also Stamen Terrain using getTiles , also on cartography, but make sure you install the latest version (released today on CRAN).

Reprex:

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)

^{Created on 2020-04-04 by the reprex package (v0.3.0)}

Vignette: https://dieghernan.github.io/cartographyvignette/

Source code: https://github.com/riatelab/cartography/blob/master/R/tilesLayer.R

Answer 2

Thanks to the answer by dieghrnan I can now get some example data, and the below now runs fine for me

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)

If you want a projection:

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

old answer:

At what step do you loose the colors? my guess is at projectRaster? I cannot get rJava to work so I cannot run this, but I would try

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