library(easyNCDF)
library(sp)
library(rgeos)
library(rgdal)
### Choose one to source the function
#source("/path/to/file/shp_mask.R")
source("/esarchive/scratch/aho/tmp/shp_mask.R")

# reference grid file
## 0.1x0.1 deg
ref.file1 <- '/esarchive/recon/ecmwf/era5land/monthly_mean/tas_f1h/tas_201006.nc'
## 1x1 deg
ref.file2 <- '/esarchive/exp/ecmwf/system5c3s/monthly_mean/tas_f6h/tas_20170201.nc' 


#============================
# NUTS
#============================
NUTS.id <- paste0("FI1D", c(1:3, 5, 7:9))
NUTS.name <- list(FI = c('Lappi', 'Kainuu'), SI = c('Pomurska', 'Podravska'))

mask1 <- shp_mask(shp.file, ref.file1, reg.ids = NUTS.id)
mask2 <- shp_mask(shp.file, ref.file2, reg.names = NUTS.name)

str(mask1)
str(mask2)


#============================
# ADM
#============================
shp.file <- '/esarchive/scratch/cdelgado/focus_outputs/Shapefiles/tza_admbnda_adm1/tza_admbnda_adm1_20181019.shp'

# Use reg.ids to specify the shape regions
mask <- shp_mask(shp.file, ref.file1, reg.ids = c("TZ07", "TZ20"), shp.system = "ADM")

which(mask == 1, arr.ind = T)
which(mask == 2, arr.ind = T)
str(attr(mask, 'index'))

# Use reg.names to specify the shape regions
mask2 <- shp_mask(shp.file, ref.file1, reg.names = list("United Republic of Tanzania" = c("Dar-es-salaam", "Mara")), shp.system = "ADM")

identical(mask, mask2)
#[1] TRUE