diff --git a/R/get_biome_classes.R b/R/get_biome_classes.R index cedbc07..aec8dba 100644 --- a/R/get_biome_classes.R +++ b/R/get_biome_classes.R @@ -20,7 +20,7 @@ get_biome_classes <- function(dataset) { nc_in <- ncdf4::nc_open(this_file) biome_attributes <- ncdf4::ncatt_get(nc_in, "biome") # format the table - if (dataset == "Beyer2020" | dataset == "Example") { + if (dataset == "Beyer2020" || dataset == "Example") { indeces_biomes <- which(substr(names(biome_attributes), 1, 5) == "biome") biome_categories <- data.frame( id = 0:(length(indeces_biomes) - 1), diff --git a/data-raw/notes_on_ecoregions.R b/data-raw/notes_on_ecoregions.R new file mode 100644 index 0000000..c362856 --- /dev/null +++ b/data-raw/notes_on_ecoregions.R @@ -0,0 +1,53 @@ +path <- "/media/andrea/Elements/resolve2017/" +# create the path if it does not exist +if (!dir.exists(path)) { + dir.create(path) +} +# download the resolve dataset +# download.file("https://storage.googleapis.com/teow2016/Ecoregions2017.zip", destfile = file.path(path, "Ecoregions2017.zip")) + +# read the zipped file +ecoregions <- terra::vect(paste0("/vsizip/",file.path(path, "Ecoregions2017.zip"))) +worldclim_pres <- pastclim::region_slice( + time_ce = 1985, + bio_variables = c("bio01"), + dataset = "WorldClim_2.1_10m" +) +ecoregions_rast <- terra::rast(terra::ext(worldclim_pres), + resolution=terra::res(worldclim_pres), + crs = terra::crs(worldclim_pres)) + +# rasterise the biomes at the resolutions of Worldclim +ecoregions_rast <- terra::rasterize(ecoregions,ecoregions_rast, field = "BIOME_NAME") +# store the levels +ecoregions_meta <- levels(ecoregions_rast) +# expand whilst there are mismatches with worldclim +while(terra::global((!is.na(worldclim_pres) & is.na(ecoregions_rast)),"sum")>0){ + ecoregions_rast <- terra::focal(ecoregions_rast, 3, "modal", na.policy="only") +} + +# now we mask by Worldclim +ecoregions_rast <- terra::mask(ecoregions_rast, worldclim_pres) + +# reassign levels +levels(ecoregions_rast) <- ecoregions_meta +time(ecoregions_rast, "years") <- time(worldclim_pres) +# now save this +ncdf_filename <- file.path(path, "ecoregions_1985.nc") +terra::writeCDF(ecoregions_rast, varname = "biome", + longname = "biome from RESOLVE Ecoregions 2017", + filename = ncdf_filename, prec = "integer", + overwrite = TRUE) +# fix a few things in the netcdf file +nc_in <- ncdf4::nc_open(ncdf_filename, write = TRUE) +ncdf4::ncatt_put(nc_in, varid = "time", attname = "units", attval = "years since 1950-01-01 00:00:00.0") +ncdf4::ncatt_put(nc_in, varid = "time", attname = "long_name", attval = "years BP") +ncdf4::ncatt_put(nc_in, varid = "time", attname = "axis", attval = "T") +format_line <- function(x){paste0(paste(x, collapse=" "), "; ")} +ncdf4::ncatt_put(nc_in, varid = "biome", attname = "biomes", + attval = paste(apply(ecoregions_meta[[1]], 1, format_line),collapse = "")) +ncdf4::nc_close(nc_in) + + +#TODO we need to save this as integers, and with info on the levels +