Newer
Older
#library(easyNCDF)
#library(sp)
#library(rgeos)
#library(rgdal)
########################
# module load GDAL PROJ GEOS
########################
# shp.file: The shp file
# ref.file: The netCDF file to provide the reference grid points
# NUTS.id: The unique ID in NUTS
# NUTS.name: A list of country and the region name
# NUTS.level: Can only have the same level in one mask
# lat_dim: 'latitude' for example
# lon_dim: 'longitude' for example
# savefile: If NULL, return an array
### Some example inputs ###
#shp.file <- paste0('/esarchive/shapefiles/NUTS3/NUTS_RG_60M_2021_4326.shp/NUTS_RG_60M_2021_4326.shp')
#ref.file <- '/esarchive/exp/ecmwf/system5c3s/monthly_mean/tas_f6h/tas_20170201.nc'
#ref.file <- '/esarchive/exp/ecmwf/s2s-monthly_ensfor/weekly_mean/tas_f6h/tas_20191212.nc'
#ref.file <- '/esarchive/recon/ecmwf/era5land/monthly_mean/tas_f1h/tas_201006.nc'
#NUTS.id <- paste0("FI1D", c(1:3, 5, 7:9))
#NUTS.name <- list(FI = c('Lappi', 'Kainuu')) #NOTE: NUTS.name may be repetitive; use level to filter
#NUTS.name <- list(FI = c('Lappi', 'Kainuu'), SI = c('Pomurska', 'Podravska'))
# mask1 <- shp_mask(shp.file, ref.file, NUTS.id = NUTS.id)
# mask2 <- shp_mask(shp.file, ref.file, NUTS.name = NUTS.name)
#NOTE: One region is one number; need to have the option to combine them?
#TODO: Suppress warnings?
#TODO: Substitute packages
#----------NEW--------------
shp_mask <- function(shp.file, ref.file, shp.system = "NUTS", ids = NULL, names = NULL,
level = 3, lat_dim = NULL, lon_dim = NULL, savefile = NULL) {
#--------NEW_END------------
####################################################
# Step 1: Load the shapefile
shp <- rgdal::readOGR(shp.file)
if (all(is.null(ids), is.null(names))) {
stop("Either provide parameter 'ids' or 'names'.")
} else if (!is.null(ids)) {
#----------NEW-------------
## Method 1: Directly use IDs
if (shp.system == "NUTS") {
shp <- subset(shp, NUTS_ID %in% ids)
} else if (shp.system == "ADM") {
shp <- subset(shp, ADM1_PCODE %in% ids)
} else {
stop("shp.system ", shp.system, " is not defined yet.")
}
#-------NEW_END-----------
if (!is.null(names)) {
warning("Only use 'ids' to get the shape region. 'names' is not used.")
}
## Method 2: Use country code & region name
for (cntr_i in 1:length(names)) {
#---------NEW------------
if (shp.system == "NUTS") {
tmp <- subset(shp, CNTR_CODE == names(names)[cntr_i])
tmp <- subset(tmp, NUTS_NAME %in% names[[cntr_i]])
} else if (shp.system == "ADM") {
tmp <- subset(shp, ADM0_EN == names(names)[cntr_i])
tmp <- subset(tmp, ADM1_EN %in% names[[cntr_i]])
}
#-------NEW_END-----------
if (cntr_i == 1) {
shp_tmp <- tmp
} else {
shp_tmp <- rbind(shp_tmp, tmp)
}
}
shp <- subset(shp_tmp, LEVL_CODE == NUTS.level)
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
}
# plot(shp)
####################################################
# Step 2: Use the reference file to get lat and lon
if (is.null(lat_dim) | is.null(lon_dim)) {
var_names <- easyNCDF::NcReadVarNames(ref.file)
lat_dim <- var_names[which(var_names %in% s2dv:::.KnownLatNames())]
lon_dim <- var_names[which(var_names %in% s2dv:::.KnownLonNames())]
}
latlon <- NcToArray(ref.file, vars_to_read = c(lat_dim, lon_dim))
lat <- NcToArray(ref.file, vars_to_read = lat_dim)[1, ]
lon <- NcToArray(ref.file, vars_to_read = lon_dim)[1, ]
## Create data frame & sp class for ref grid
ref.df <- data.frame(data = 0,
lon = rep(lon, times = length(lat)),
lat = rep(lat, each = length(lon)))
coord <- data.frame(x = ref.df$lon, y = ref.df$lat)
ref.sp <- SpatialPointsDataFrame(coord,
data = data.frame(ref.df$data))
proj4string(ref.sp) <- sp::proj4string(shp)
####################################################
# Step 3: Create mask
## Create mask array with 0; 1, 2, etc. will be filled in for each shp region
mask <- array(0, dim = c(length(lon), length(lat)))
names(dim(mask)) <- c(lon_dim, lat_dim)
## Loop through each shp region
for (shp_i in 1:nrow(shp)) {
shp_pol <- gIntersection(ref.sp, shp[shp_i, ])
if (!is.null(shp_pol)) {
for (ii in 1:nrow(shp_pol@coords)) {
mask[which(lon == shp_pol@coords[ii, 1]), which(lat == shp_pol@coords[ii, 2])] <- shp_i
}
} else {
warning("shp ID '", shp$NUTS_ID[shp_i], "' cannot be identified in this reference grid.")
}
}
####################################################
# Step 4: Add attributes
attr(mask, lon_dim) <- lon
attr(mask, lat_dim) <- lat
#---------NEW------------
if (shp.system == "NUTS") {
attr(mask, "index") <- as.list(shp$NUTS_ID)
names(attr(mask, "index")) <- 1:nrow(shp)
} else if (shp.system == "ADM") {
attr(mask, "index") <- as.list(shp$ADM1_PCODE)
names(attr(mask, "index")) <- 1:llength(shp)
}
#---------NEW_END------------
## Return all the info from shp
attr(mask, "shapefile") <- attributes(shp)
####################################################
# Step 5: Save the file or return the array
if (is.null(savefile)) {
return(mask)
} else {
#TODO
ArrayToNc()
}
}