Newer
Older
#'Convert Shapefile to Mask Array
#'
#'This function reads a shapefile (.shp) containing information about polygonal
#'regions. It then transfers the shapefile data into an array and subsets the
#'output based on requested region names or IDs. The accepted shapefile
#'databases are 'NUTS', 'LAU', and 'GADM', each with its own unique format.
#'However, the function can use other shapefiles databases with specifying the
#'
#'To ensure accurate comparison with the shapefile, the function loads a
#'reference dataset that provides longitude and latitude information. By
#'intersecting each subset of the shapefile with the reference coordinates, the
#'function selects only the desired regions. The final step involves creating a
#'mask array. Depending on the chosen option, the mask array is either returned
#'as the function's output or saved into a NetCDF format in the specified
#'directory.
#'
#'Note: Modules GDAL, PROJ and GEOS are required.
#'
#'@param shp_file A character string indicating the shp file path.
#'@param ref_grid A character string indicating the path to the reference
#' data. Either (1) a netCDF file or (2) a list of lon and lat to provide the
#' reference grid points. It is NULL by default.
#'@param compute_area_coverage A logical value indicating the method to find
#' the intersection of the reference grid and the shapefile. When it is TRUE,
#' the method used is the calculation of the area coverage fraction of
#' intersection. If it is FALSE, the method used is searching if the centroid
#' of the grid cell falls inside the shapefile or not. It is FALSE by default.
#'@param shp_system A character string containing the Shapefile System Database
#' Name used to subset the shapefile into regions by using parameters 'reg_ids'
#' or 'reg_names'. The accepted systems are: 'NUTS', 'LAU', and 'GADM'. When it
#' is used, you must specify either 'reg_ids' or 'reg_names'; if you don't need
#' to subset different regions, set it to NULL. It is set to "NUTS" by default
#' (optional).
#'@param shp_col_name_ids A character string indicating the column name of the
#' column in where the specified 'reg_ids' will be taken (optional).
#'@param reg_ids A character string indicating the unique ID in shapefile.
#' It is NULL by default (optional).
#'@param reg_names A named list of character string vectors indicating the
#' country and the region name. The name of the list stands for the country
#' name code and the vector character strings indicate the region name for
#' each country. It is NULL by default (optional).
#'@param reg_level An integer number from 1 to 3 indicating the 'NUTS' dataset
#' level. For other datasets this parameter is not used. One mask can only have
#' a unique level. It is set to 3 by default (optional).
#'@param lat_dim A character string indicating the latitudinal dimension. If it
#' is NULL, the latitudinal name will be searched using an internal function
#' with the following possible names: 'lat', 'latitude', 'y', 'j' and
#'@param lon_dim A character string indicating the longitudinal dimension. If it
#' is NULL, the longitudinal name will be searched using an internal function
#' with the following possible names: 'lon', 'longitude', 'x', 'i' and
#'@param region A logical value indicating if we want a dimension for the
#' regions in the resulting mask array. It is FALSE by default.
#'@param check_valid A logical value that when it is TRUE it uses the function
#' 'sf::st_make_valid' applied to the shapefile and to the coordinates.
#'@param find_min_dist A logical value indicating if we want to look for the
#' nearest coordinate between the shapefile region and the reference grid when
#' there is no intersection between the shapefile and the reference grid. It is
#' FALSE by default.
#'@param max_dist A numeric value indicating the maximum distance is accepted
#' to the closest gridpoint when there is no intersection between the shapefile
#' and the reference grid.
#'@param ncores The number of parallel processes to spawn for the use for
#' parallel computation in multiple cores.
#'@param fileout A character string of the path to save the NetCDF mask. If not
#' specified (default), only the mask array will be returned.
#'@param ... Arguments passed on to 's2_options' in function 'st_intersection'.
#' See 's2 package'.
#'@param units A character string indicating if your GIS files has a grid in degrees or meters. If it
#' is NULL, the units will be set as "meters"
#' with the following possible names: 'degrees', 'meters'
#'@param name_shp_col A character string indicating in the shape file which is the name of the column
#' with the values of the names of the different polygons. It is NULL by default.
#'@param id_shp_col A character string indicating in the shape file which is the name of the column
#' with the values of the IDs of the different polygons. It is NULL by default.
#'
#'@return A multidimensional array containing a mask array with longitude and
#'latitude dimensions. If 'region' is TRUE, there will be a dimension for
#'the region.
#'
#'@examples
#'# Exmple (1): NUTS
#'shp_file <- paste0('/esarchive/shapefiles/NUTS3/NUTS_RG_60M_2021_4326.shp/',
#'ref_grid <- list(lon = seq(10, 40, 0.5), lat = seq(40, 85, 0.5))
#'NUTS_name <- list(FI = c('Lappi', 'Kainuu'), SI = c('Pomurska', 'Podravska'))
#'mask2 <- ShapeToMask(shp_file = shp_file, ref_grid = ref_grid,
#' reg_names = NUTS_name)
#'@importFrom doParallel registerDoParallel
ShapeToMask <- function(shp_file, ref_grid, compute_area_coverage = FALSE,
shp_system = "NUTS", reg_names = NULL, reg_ids = NULL,
shp_col_name_ids = NULL, reg_level = 3,
region = FALSE, check_valid = FALSE,
find_min_dist = FALSE, max_dist = 50, ncores = NULL,
fileout = NULL, units = 'degrees', id_shape_col = NULL,
name_shape_col = NULL, ...) {
# TODO: Add saving option?
# Initial checks
# shp_file
if (is.null(shp_file)) {
stop("Parameter 'shp_file' cannot be NULL.")
}
if (!is.character(shp_file)) {
stop("Parameter 'shp_file' must be a character string.")
}
# lon_dim, lat_dim
if (!is.null(lon_dim)) {
if (!is.character(lon_dim)) {
stop("Parameter 'lon_dim' must be a character string.")
}
}
if (!is.null(lat_dim)) {
if (!is.character(lat_dim)) {
stop("Parameter 'lat_dim' must be a character string.")
}
}
if (!is.null(id_shape_col)) {
if (!is.character(id_shape_col)) {
stop("Parameter 'id_shape_col' must be a character string.")
}
}
if (!is.null(name_shape_col)) {
if (!is.character(name_shape_col)) {
stop("Parameter 'name_shape_col' must be a character string.")
}
}
if (!is.null(units)) {
if (!is.character(units)) {
stop("Parameter 'units' must be a character string.")
}
}
# ref_grid
if (is.null(ref_grid)) {
stop("Parameter 'ref_grid' cannot be NULL.")
}
if (!any(is.character(ref_grid) | is.list(ref_grid))) {
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
stop("Parameter 'ref_grid' must be either a netCDF file or a list of lon and lat.")
}
if (is.character(ref_grid)) {
if (!file.exists(ref_grid)) {
stop("ref_grid file does not exist.")
}
}
if (is.list(ref_grid)) {
if (length(ref_grid) != 2) {
stop("If 'ref_grid' is a list, it must have two items for longitude and latitude.")
}
if (is.null(lat_dim) | is.null(lon_dim)) {
lon_known_names <- c(.KnownLonNames(), 'lons')
lat_known_names <- c(.KnownLatNames(), 'lats')
lon_dim <- lon_known_names[which(lon_known_names %in% names(ref_grid))]
lat_dim <- lat_known_names[which(lat_known_names %in% names(ref_grid))]
if (identical(lon_dim, character(0)) | identical(lat_dim, character(0))) {
stop("longitude and latitude names are not recognized in 'ref_grid'. ",
"Please specify 'lon_dim' and 'lat_dim'.")
}
}
}
# shp_system, reg_ids, reg_names, shp_col_name_ids
if (!is.null(shp_system)) {
if (!is.character(shp_system)) {
stop("Parameter 'shp_system' must be a character strinig.")
}
if (all(is.null(reg_ids), is.null(reg_names))) {
stop("If 'shp_system' is used, you must provide either parameter ",
"'reg_ids' or 'reg_names'.")
} else if (!is.null(shp_col_name_ids)) {
if (is.null(reg_ids)) {
stop("If 'shp_col_name_ids' is used, parameter 'reg_ids' must be provided.")
}
if (!is.character(shp_col_name_ids)) {
stop("Parameter 'shp_col_name_ids' must be a character strinig.")
}
}
if (all(!is.null(reg_ids), !is.null(reg_names))) {
warning("Only use 'reg_ids' to get the shape region. 'reg_names' is not used.")
if (!is.null(reg_level)) {
if (!is.numeric(reg_level)) {
stop("Parameter 'reg_level' must be numeric.")
}
}
if (!is.logical(region)) {
stop("Parameter 'region' must be a logical value.")
}
if (!is.logical(check_valid)) {
stop("Parameter 'check_valid' must be a logical value.")
}
if (!is.logical(find_min_dist)) {
stop("Parameter 'find_min_dist' must be a logical value.")
}
if (!is.null(max_dist)) {
if (!is.numeric(max_dist)) {
stop("Parameter 'max_dist' must be a numeric.")
}
}
# ncores
if (!is.null(ncores)) {
if (!is.numeric(ncores)) {
stop("Parameter 'ncores' must be numeric.")
}
ncores <- round(ncores)
if (ncores < 2) {
ncores <- NULL
}
}
shp <- sf::st_read(shp_file) # class sf
if (units == 'degrees') {
shp <- sf::st_transform(shp, 4326)
} else if (units == 'meters') {
shp <- sf::st_transform(shp, 3857)
}
shp_crs <- sf::st_crs(shp)
NUTS_ID <- ADM1_PCODE <- ISO <- NULL
CNTR_CODE <- NUTS_NAME <- ADM0_EN <- ADM1_EN <- Name <- LEVL_CODE <- NULL
if (!is.null(shp_col_name_ids)) {
if (shp_col_name_ids %in% names(shp)) {
shp <- subset(shp, get(shp_col_name_ids) %in% reg_ids)
} else {
stop("Shape system name not found in shapefile names.")
}
} else if (shp_system == "NUTS") {
shp <- subset(shp, NUTS_ID %in% reg_ids)
} else if (shp_system == "ADM") {
shp <- subset(shp, ADM1_PCODE %in% reg_ids)
} else if (shp_system == "GADM") {
shp <- subset(shp, ISO %in% reg_ids)
stop("shp_system ", shp_system, " is not defined yet.")
} else if (!is.null(reg_names)) {
for (cntr_i in 1:length(reg_names)) {
if (shp_system == "NUTS") {
tmp <- subset(shp, CNTR_CODE == names(reg_names)[cntr_i])
tmp <- subset(tmp, NUTS_NAME %in% reg_names[[cntr_i]])
} else if (shp_system == "ADM") {
tmp <- subset(shp, ADM0_EN == names(reg_names)[cntr_i])
tmp <- subset(tmp, ADM1_EN %in% reg_names[[cntr_i]])
} else if (shp_system == "GADM") {
tmp <- subset(shp, Name %in% reg_names)
}
if (cntr_i == 1) {
shp_tmp <- tmp
} else {
shp_tmp <- rbind(shp_tmp, tmp)
}
}
if (shp_system == "NUTS") {
shp <- subset(shp_tmp, LEVL_CODE == reg_level)
} else if (shp_system == "ADM" | shp_system == "GADM") {
# Step (2.1): Use the reference file to get lat and lon
if (all(tools::file_ext(ref_grid) == 'nc')) {
## Method 1: ref_grid is a netCDF file
if (is.null(lat_dim) | is.null(lon_dim)) {
var_names <- easyNCDF::NcReadVarNames(ref_grid)
lat_dim <- var_names[which(var_names %in% .KnownLatNames())]
lon_dim <- var_names[which(var_names %in% .KnownLonNames())]
latlon <- easyNCDF::NcToArray(ref_grid, vars_to_read = c(lat_dim, lon_dim))
lat <- easyNCDF::NcToArray(ref_grid, vars_to_read = lat_dim)[1, ]
lon <- easyNCDF::NcToArray(ref_grid, vars_to_read = lon_dim)[1, ]
} else if (is.list(ref_grid)) {
## Method 2: ref_grid is a list of lon and lat
lat <- ref_grid[[lat_dim]]
lon <- ref_grid[[lon_dim]]
# Set degrees longitudes to -180 to 180
lon <- ifelse(lon > 180, lon - 360, lon)
order_idx <- order(lon)
lon <- lon[order_idx]
locations <- array(1:(length(lon)*length(lat)), c(lon = length(lon),
lat = length(lat)))
# Build data dataframe
lonlat_df <- data.frame(lon = rep(as.vector(lon), length(lat)),
lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE))
data_df <- lonlat_df %>%
mutate(dat = as.vector(locations))
lonlat_df_ori <- NULL
# NOTE: if target_proj = "ESRI:54030", Nord3v2 has different behavior from hub and ws!!
data_df <- st_as_sf(data_df, coords = c("lon", "lat"), crs = shp_crs)
data_df <- st_transform(data_df, crs = shp_crs)
data_df <- data_df %>%
mutate(long = st_coordinates(data_df)[, 1],
lat = st_coordinates(data_df)[, 2])
xy.sfc <- polygonize(lonlat_df = lonlat_df, data_df = data_df,
lon = lon, lat = lat, target_proj = shp_crs,
original_proj = shp_crs)
region <- TRUE
} else {
ref.df <- data.frame(data = 0,
lon = rep(lon, times = length(lat)),
lat = rep(lat, each = length(lon)))
ref.df <- modifyLonsTo180(units, ref.df, order_idx)
coord <- as.matrix(data.frame(x = ref.df$lon, y = ref.df$lat))
xy.sfg <- sf::st_multipoint(coord)
xy.sfc <- sf::st_sfc(xy.sfg)
# Assign crs of original shapefile
st_crs(xy.sfc) <- sf::st_crs(shp) # asign crs of shapefile
# Check valid
if (check_valid) {
xy.sfc <- st_make_valid(xy.sfc)
shp <- st_make_valid(shp)
}
## Step (3): Loop through each shp region to create the mask
if (region) {
cfun <- function(a, b) {
if (length(dim(a)) == 2) {
dims <- c(dim(a), 2)
dims <- c(dim(a)[1:2], dim(a)[length(dim(a))]+1)
return(array(c(a,b), dim = dims))
} else {
cfun <- function(a, b) {
a[which(b != 0)] <- b[which(b != 0)]
return(a)
}
}
if (is.null(ncores)) {
mask <- foreach(shp_i = 1:nrow(shp), .combine = 'cfun') %do%
.shapetomask(shp = shp, n = shp_i, lon = lon, lat = lat,
xy.sfc = xy.sfc, compute_area_coverage = compute_area_coverage,
max_dist = max_dist, region = region, shp_system = shp_system, ...)
} else {
registerDoParallel(ncores)
mask <- foreach(shp_i = 1:nrow(shp), .combine = 'cfun', .packages = c('sf', 'dplyr')) %dopar%
.shapetomask(shp = shp, n = shp_i, lon = lon, lat = lat,
xy.sfc = xy.sfc, compute_area_coverage = compute_area_coverage,
max_dist = max_dist, region = region, shp_system = shp_system, ...)
registerDoSEQ()
}
if (region) {
if (length(dim(mask)) == 2) dim(mask) <- c(dim(mask), 1)
names(dim(mask)) <- c(lon_dim, lat_dim, 'region')
} else {
names(dim(mask)) <- c(lon_dim, lat_dim)
attr(mask, lon_dim) <- lon
attr(mask, lat_dim) <- lat
attr(mask, "index") <- 1:nrow(shp)
## Return all the info from shp
attr(mask, "shapefile") <- attributes(shp)
# Generate attributes related with the shape file
shape_attrs = list()
if (is.null(id_shape_col) && is.null(name_shape_col)) {
if (shp_system == "NUTS") {
id_aux = shp$NUTS_ID[i]
name_aux = shp$NUTS_NAME[i]
} else if (shp_system == "ADM") {
id_aux = shp$ADM1_PCODE[i]
name_aux = shp$ADM1_NAME[i]
} else if (shp_system == "GADM") {
id_aux = shp$ISO[i]
name_aux = shp$Name[i]
} else {
id_aux = NULL
name_aux = NULL
}
} else {
id_aux = id_shape_col
name_aux = name_shape_col
id = id_aux,
name = name_aux,
shape_val = i
)
shape_attrs <- append(shape_attrs, list(entry))
}
# Step 5: Save to NetCDF
if (!is.null(fileout)) {
# Lat Lon Max Min
lat_list <- attr(mask, 'lat')
lon_list <- attr(mask, 'lon')
transformToNc (mask, lat_list, lon_list, NULL, fileout, units, shape_attrs, lat_dim, lon_dim)
.shapetomask <- function(shp, n, lon, lat, xy.sfc, compute_area_coverage = FALSE,
region = FALSE, shp_system = 'NUTS', ...) {
if (compute_area_coverage) {
mask <- xy.sfc %>%
dplyr::mutate(int = areacov(geometry, shpi)) %>%
dplyr::arrange(value) %>%
dplyr::pull(int)
dim(mask) <- c(length(lon), length(lat))
} else {
mask <- array(0, dim = c(length(lon), length(lat)))
# NOTE: Don't know it's a problem in st_intersection or st_coordinates, tmp_coords may
# not be identical as lon/lat. E.g., (29, 65) may be (29 - -3.552714e-15, 65).
shp_pol <- sf::st_intersection(x = xy.sfc, y = shpi, ...)
tmp_coords <- sf::st_coordinates(shp_pol)[, 1:2]
if (length(tmp_coords) == 0) {
dim(tmp_coords) <- NULL
} else if (is.null(dim(tmp_coords))) {
tmp_coords <- array(tmp_coords, dim = c(1, length(tmp_coords)))
}
if (!is.null(dim(tmp_coords))) {
# polygon_instersection
for (ii in 1:nrow(tmp_coords)) {
# pt_x <- which(lon == tmp_coords[ii, 1])
# pt_y <- which.min(abs(lat - tmp_coords[ii, 2]))
if (!region) {
# min(abs(lon - tmp_coords[ii, 1]))
# min(abs(lat - tmp_coords[ii, 2]))
mask[which.min(abs(lon - tmp_coords[ii, 1])),
which.min(abs(lat - tmp_coords[ii, 2]))] <- n # AQUI ES DONDE ASOCIA EL N A CADA VALOR DEL SHAPE PROBABLEMENTE RAUL
} else {
mask[which.min(abs(lon - tmp_coords[ii, 1])),
which.min(abs(lat - tmp_coords[ii, 2]))] <- 1 # SI SOLO HAY UNA REGION O NO HAY QUE DIVIDIR POR REGIONES
} else if (find_min_dist) {
x.centroid.shpi <- sf::st_coordinates(sf::st_centroid(shpi))[,1]
y.centroid.shpi <- sf::st_coordinates(sf::st_centroid(shpi))[,2]
dist <- sqrt((xy.sfc[,1] - x.centroid.shpi)**2 + (xy.sfc[,2] - y.centroid.shpi)**2)
tmp_coords <- array(xy.sfc[which(dist == min(dist, na.rm = TRUE)),], dim = c(1,2))
colnames(tmp_coords) <- c('X', 'Y')
if (max(dist) <= max_dist & (any(round(lat,2) == round(tmp_coords[1,2],2)) &
any(round(lon,2) == round(tmp_coords[1,1],2))) ) {
if (length(dim(mask)) == 2) {
mask[which.min(abs(lon - tmp_coords[, 1])),
which.min(abs(lat - tmp_coords[, 2]))] <- n
} else {
mask[which.min(abs(lon - tmp_coords[, 1])),
which.min(abs(lat - tmp_coords[, 2]))] <- 1
}
# warning(paste0('The reference grid has no intersection with region ',
# ifelse(is.character(shp_col_name_ids), shpi[[shp_col_name_ids]], paste0('n° ', n)),
# ' from the shapefile; the provided grid cell is at a distance of ', dist[which(dist == min(dist, na.rm = TRUE))],
# ' to the centroid of the region (units are: ° or meters depending on the crs of the shapefile).'))
# warning(paste0('The reference grid has no intersection with region ',
# ifelse(is.character(shp_col_name_ids), shpi[[shp_col_name_ids]], paste0('n° ', n))))
# Flip the matrix when the area_coverage is TRUE
if (compute_area_coverage && grepl("NUTS", shp_system)) {
mask <- mask[, ncol(mask):1] # Flip in the horizontal axis
# mask <- mask[nrow(mask):1, ] # Flip in vertical axis
}
# Function to create polygons from VizRobinson code
polygonize <- function(lonlat_df, data_df, lon, lat, target_proj,
original_proj) {
# Calculate polygon points from regular lat/lon
# NOTE: The original grid must be regular grid with same space
d_lon <- abs(lon[2] - lon[1]) / 2
d_lat <- abs(lat[2] - lat[1]) / 2
lon_poly <- lat_poly <- rep(NA, 4 * dim(lonlat_df)[1])
for (ii in 1:dim(lonlat_df)[1]) {
lon_poly[(ii*4-3):(ii*4)] <- c(lonlat_df$lon[ii] - d_lon, lonlat_df$lon[ii] + d_lon,
lonlat_df$lon[ii] + d_lon, lonlat_df$lon[ii] - d_lon)
lat_poly[(ii*4-3):(ii*4)] <- c(lonlat_df$lat[ii] - d_lat, lonlat_df$lat[ii] - d_lat,
lonlat_df$lat[ii] + d_lat, lonlat_df$lat[ii] + d_lat)
}
# To prevent out-of-global lon
# lon_poly[which(lon_poly >= 180)] <- 179.9
# lon_poly[which(lon_poly < -180)] <- -180
# To prevent out-of-global lat
lat_poly[which(lat_poly > 90)] <- 90
lat_poly[which(lat_poly < -90)] <- -90
lonlat_df <- data.frame(lon = lon_poly, lat = lat_poly)
# Transfer lon/lat to projection
proj_lonlat <- st_as_sf(lonlat_df, coords = c("lon", "lat"), crs = original_proj)
# NOTE: if target_proj = "ESRI:54030", on Nord3v2, st_transform has lon and lat swapped!
proj_lonlat <- st_transform(proj_lonlat, crs = target_proj)
lonlat_df_proj <- st_coordinates(proj_lonlat)
# Use id to create groups for each polygon
ids <- factor(paste0("id_", 1:dim(data_df)[1]))
values <- data.frame(id = ids,
value = data_df$dat)
positions <- data.frame(id = rep(ids, each = 4),
x = lonlat_df_proj[, 1],
y = lonlat_df_proj[, 2])
datapoly <- merge(values, positions, by = "id")
datapoly <- st_as_sf(datapoly, coords = c("x", "y"), crs = target_proj)
dplyr::group_by(.data$id) %>%
dplyr::summarise(geometry = st_combine(geometry)) %>%
# dplyr::summarise() %>% # NOTE: VERY SLOW if plot global
dplyr::mutate(value = values[order(values$id), ]$value) %>%
st_cast("POLYGON") %>%
st_convex_hull() # maintain outer polygen (no bowtie shape)
# Function to compute the coverage area ratio between the grid and the shapefile
areacov <- function(grid, shp) {
shp_int_grid <- sf::st_intersection(x = grid, y = shp)
idx <- attributes(shp_int_grid)$idx
i <- 0
fr <- array(0, length(grid))
for (idx_i in idx[, 1]) {
i <- i + 1
area_grid <- sf::st_area(grid[idx_i])
area_int <- sf::st_area(shp_int_grid[i])
fr[idx_i] <- round(as.numeric(area_int/area_grid), digits = 6)
}
return(fr)
transformToNc <- function(values, lat=NULL, lon=NULL, time=NULL, out_file, crs='degrees', shape_values = NULL, lat_dim, lon_dim) {
# Take the name of lat lon dimensions
if (is.null(lon_dim)) {
possible_lon_names <- c("lon", "longitude", "x", "i", "nav_lon")
} else {
possible_lon_names <- c(lon_dim)
}
if (is.null(lat_dim)) {
possible_lat_names <- c("lat", "latitude" , "y", "j", "nav_lat")
} else {
possible_lat_names <- c(lat_dim)
}
# Change latitude and longitude names to lat, lon
names(dim(values)) <- gsub(paste(possible_lon_names, collapse = "|"), "lon", names(dim(values)), ignore.case = TRUE)
names(dim(values)) <- gsub(paste(possible_lat_names, collapse = "|"), "lat", names(dim(values)), ignore.case = TRUE)
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
# Add coordinates
if (!is.null(lon) && length(lon) > 0) {
dim(lon) <- length(lon)
metadata <- list(lon = list(units = crs))
attr(lon, 'variables') <- metadata
names(dim(lon)) <- 'lon' # rigth_lon_name
}
if (!is.null(lat) && length(lat) > 0) {
dim(lat) <- length(lat)
metadata <- list(lat = list(units = crs))
attr(lat, 'variables') <- metadata
names(dim(lat)) <- 'lat' # rigth_lat_name
}
# Add the projection to the file
if (crs == 'degrees') {
crs = 'EPSG4326'
} else {
crs = 'EPSG3857'
}
# Format the attributes before introducing to the file
formatted_string <- paste(shape_values, collapse="")
formatted_string <- gsub("\\(", "", formatted_string)
formatted_string <- gsub("\\)", "\n", formatted_string)
formatted_string <- gsub("list", "", formatted_string)
formatted_string <- gsub('\\"', "", formatted_string)
# Add the atributes to the file
attrs <- list(global_attrs = list(
crs = crs,
shapes_values = formatted_string
))
attributes(values) <- c(attributes(values), attrs)
# Export to Nc
ArrayToNc(list(values, lon, lat), out_file)
}
modifyLonsTo180 <- function(units, df, order_idx) {
print('Check errors 1:')
# browser()
if (units == 'degrees'){
my_list <- c("lon", "longitude", "x")
data <- df[colnames(df)[1]]
name_lon <- my_list[tolower(my_list) %in% tolower(colnames(df))]
lon_data <- df[[name_lon]]
data_reorder <- df[[colnames(df)[1]]]
df[colnames(df)[1]] <- data_reorder[order_idx]
}
return (df)
}