library(easyNCDF)
library(sf)
#library(sp)
#library(rgeos)
### 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.grid1 <- '/esarchive/recon/ecmwf/era5land/monthly_mean/tas_f1h/tas_201006.nc'
## 1x1 deg
ref.grid2 <- '/esarchive/exp/ecmwf/system5c3s/monthly_mean/tas_f6h/tas_20170201.nc' 
## list of lon and lat
ref.grid_FI <- list(lon = seq(10, 40, 0.5), lat = seq(40, 85, 0.5)) 

#============================
# NUTS
#============================
NUTS.id <- paste0("FI1D", c(1:3, 5, 7:9))
NUTS.name <- list(FI = c('Lappi', 'Kainuu'), SI = c('Pomurska', 'Podravska'))
shp.file <- paste0('/esarchive/shapefiles/NUTS3/NUTS_RG_60M_2021_4326.shp/NUTS_RG_60M_2021_4326.shp')
mask1 <- shp_mask(shp.file, ref.grid1, reg.ids = NUTS.id)
mask2 <- shp_mask(shp.file, ref.grid2, reg.names = NUTS.name)
mask3 <- shp_mask(shp.file, ref.grid_FI, reg.names = NUTS.id)

str(mask1)
str(mask2)
unique(sort(c(mask3)))
which(mask3 != 0, arr.ind=T) 

#============================
# 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.grid1, 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.grid1, reg.names = list("United Republic of Tanzania" = c("Dar-es-salaam", "Mara")), shp.system = "ADM")

identical(mask, mask2)
#[1] TRUE