Search code examples
rgeospatialmaskshapefileterra

Unable to mask spatial object even after transforming the CRS of the shapefile to match the CRS of the raster

I'm trying to mask a raster using a shapefile in R, specifically using the terra package. However, I'm encountering an issue where the output after masking is all NaN values. I've followed the documentation and made sure that the coordinate reference systems (CRS) of both the raster and shapefile match, but the problem persists.

By looking at the two plots, I believe the cause of the problem is that the extension of the raster is incorrect. The latitude should start at -180 up to 180. The raster and shapefiles I am using can be obtain in the following link.

# ---- NOAA global temperature data ----
# Load the terra package
if (!require(terra)) install.packages('terra', repos='https://rspatial.r-universe.dev')

# Load NOAA global temperature data
tmax_2018 <- rast("data/temperature/tmax.2018.nc")

tmax_2018_Jan1 <- tmax_2018[[1]]

> tmax_2018_Jan1
class       : SpatRaster 
dimensions  : 360, 720, 1  (nrow, ncol, nlyr)
resolution  : 0.5, 0.5  (x, y)
extent      : 0, 360, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source      : tmax.2018.nc 
varname     : tmax (Daily Maximum Temperature) 
name        : tmax_1 
unit        :   degC 
time        : 2018-01-01 UTC

enter image description here

# ---- Ecuador shapefile ----
# Load sf package
if (!require("sf")) install.packages("sf")

# Load Ecuador shapefiles
ecuador_shp <- st_read('data/ecu_shapefiles/ECU_adm1.shp')
ecuador_shp <- st_simplify(ecuador_shp, preserveTopology = TRUE, dTolerance = 100)

# Change the CRS to match the SpatRaster's CRS
ecuador_shp <- st_transform(ecuador_shp, crs(tmax_2018_Jan1))

# Transforming sf class to SpatVector class from terra
ecuador_shp <- vect(ecuador_shp)

> ecuador_shp
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 24, 9  (geometries, attributes)
 extent      : -92.00854, -75.20007, -5.017001, 1.681093  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 
 names       :  ID_0   ISO  NAME_0  ID_1  NAME_1    TYPE_1 ENGTYPE_1 NL_NAME_1 VARNAME_1
 type        : <num> <chr>   <chr> <num>   <chr>     <chr>     <chr>     <chr>     <chr>
 values      :    68   ECU Ecuador     1   Azuay Provincia  Province        NA        NA
                  68   ECU Ecuador     2 Bolivar Provincia  Province        NA        NA
                  68   ECU Ecuador     3   Cañar Provincia  Province        NA        NA

enter image description here

# ---- Masking ----
tmax_2018_Jan1_ecu <- terra::mask(tmax_2018_Jan1, ecuador_shp)
plot(tmax_2018_Jan1) # Shows an empty plot.
Solution

  • Here is (first) a simplified and condensed version of margusl's answer.

    library(terra)
tmax_2018 <- rast("tmax.2018.nc")
tmax_2018 <- rotate(tmax_2018)

ecu <- vect('ecu_shapefiles/ECU_adm1.shp')
tmax_2018_ecu <- crop(tmax_2018, ecu, mask=TRUE)

plot(tmax_2018_ecu, 1)
lines(ecu)

    If you do not want to rotate the raster data, you could shift the vector data instead.

    library(terra)
tmax_2018 <- rast("tmax.2018.nc")

ecu <- vect('ecu_shapefiles/ECU_adm1.shp')
secu <- shift(secu, 360)
tmax_2018_secu <- crop(tmax_2018, secu, mask=TRUE)

plot(tmax_2018_secu, 1)
lines(secu)

    And you could then shift the result. That could be more efficient if that is an issue.

    tmax <- shift(tmax_2018_secu, -360)