Scripts/01_PreparePolygonFiles.R
View file @ 822fac28
......@@ -14,6 +14,25 @@ crs.laea <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS
crs.mollweide <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" # For global equal area projections
# ------------- #
#### Split EEA reference grid for a harmonized grid for rasterization ####
# Ideally we ensure that the reference grid is consistent with the
# EEA reference grid https://www.eea.europa.eu/data-and-maps/data/eea-reference-grids-2
# Since the reference grid is for Europe only available at a 10km resolution,
# we will use the original EEA grid and split it into smaller parts.
# Download European reference grids
# Source: https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/grids
for(grain in c("1km", "2km", "5km", "10km", "20km", "50km", "100km")){
print(grain)
url <- paste0("https://gisco-services.ec.europa.eu/grid/grid_",grain,"_surf.gpkg")
curl::curl_download(
url,
paste0(path_input, basename(url) )
)
}
#### Prepare European countries layer ####
# Specifically take the various European countries and
......@@ -26,14 +45,3 @@ write_sf(sub, "/home/martin/GADM_Europesubset.gpkg")
# --- #
# --- #
# Download European reference grids
# Source: https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/grids
grid_1km <- "https://gisco-services.ec.europa.eu/grid/grid_1km_surf.gpkg"
curl::curl_download(
grid_1km,
paste0(path_input, basename(grid_1km) )
)
grid <- sf::st_read( paste0(path_input, "grid_1km_surf.gpkg"))
Scripts/03_PrepareEuropeanGrid.R
View file @ 822fac28
library(sf)
library(terra)
terra::terraOptions(progress = TRUE)
terra::terraOptions(
progress = TRUE,
tempdir = "/media/martin/data/tmp/")
# Paths
path_input <- "Inputs/" # Note that inputs and other files need to be downloaded separately
......@@ -15,7 +17,8 @@ crs.laea <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS
#### Rasterize the EEA European Coastline grid ####
# Source:
# https://www.eea.europa.eu/data-and-maps/data/eea-coastline-for-analysis-1/gis-data/europe-coastline-shapefile
coast <- sf::st_read(paste0(path_input, "Europe_coastline_shapefile/Europe_coastline_poly.shp"))
coast <- sf::st_read(paste0(path_input, "Europe_coastline_shapefile/Europe_coastline_poly.shp"), quiet = TRUE) |>
sf::st_transform(crs = sf::st_crs(crs.laea))
# Load GADM dataset for global shorelines
gadm <- read_sf("/mnt/pdrive/bec_data/100_rawdata/GADM_AdminBoundaries/gadm_410.gpkg", quiet = T) |>
......@@ -28,11 +31,17 @@ gadm <- read_sf("/mnt/pdrive/bec_data/100_rawdata/GADM_AdminBoundaries/gadm_410.
# Also get the Grid from CHELSA for its extent
# This is larger than Europe and should incorporate all the areas modelled by WP3
ext <- rast(paste0(path_input, "Mask5km_ENV/Mask5km_ENV/Mask5Km.tif")) |> st_bbox()
# ext <- rast(paste0(path_input, "Mask5km_ENV/Mask5km_ENV/Mask5Km.tif")) |> st_bbox()
# Use the 10km EEA grid as reference instead
grid <- sf::st_read("Inputs/EU_ReferenceGrid/europe_10km.shp", quiet = TRUE)
ext <- sf::st_bbox(grid)
# Make a template raster at 100m for Europe
template <- rast(extent = ext,
crs = crs.laea, resolution = c(100,100))
# Make a template raster at 1000m for Europe using the surf grid as reference
template <- rast(extent = ext, crs = crs.laea, resolution = c(1000, 1000), vals = 1)
baseline <- rasterize(grid, template)
# Now disaggregate to 100m. Since they are based on the same extent this works
template <- terra::disagg(template, fact = 10, method = "near", datatype = "INT2U")
# Take at minimum the extent from the coastline layer and rasterize it to the template
ras <- terra::rasterize(coast, template, touches = TRUE)
......@@ -43,12 +52,12 @@ ras_gadm <- terra::rasterize(gadm, template, touches = TRUE)
ras2 <- sum(ras, ras_gadm, na.rm = TRUE)
ras2[ras2>=1] <- 1
terra::writeRaster(ras2, paste0(path_output,"/ReferenceGrid_Europe_100m.tif"),
terra::writeRaster(ras2, paste0(path_output,"/Gridded/ReferenceGrid_Europe_100m.tif"),
datatype = "INT1U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
rm(ras2); gc()
# ----- #
ras <- rast( paste0(path_output,"ReferenceGrid_Europe_100m.tif") )
ras <- rast( paste0(path_output,"Gridded/ReferenceGrid_Europe_100m.tif") )
# Next we mean aggregate this layer to 1km #
ras[is.na(ras)] <- 0 # First set all na data to 0 to allow fractional aggregation
ras_1km <- terra::aggregate(ras, fact = 10, fun = "mean", cores = 7)
......@@ -57,51 +66,89 @@ ras_1km <- terra::aggregate(ras, fact = 10, fun = "mean", cores = 7)
ras_1km <- round(ras_1km * 10000)
# Write the output as fractional layer
terra::writeRaster(ras_1km, paste0(path_output,"ReferenceGrid_Europe_frac_1000m.tif"),
terra::writeRaster(ras_1km, paste0(path_output,"Gridded/ReferenceGrid_Europe_frac_1000m.tif"),
datatype = "INT4U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
ras_1km_bin <- ras_1km
ras_1km_bin[ras_1km_bin>0] <- 1
terra::writeRaster(ras_1km_bin, paste0(path_output,"ReferenceGrid_Europe_bin_1000m.tif"),
terra::writeRaster(ras_1km_bin, paste0(path_output,"Gridded/ReferenceGrid_Europe_bin_1000m.tif"),
datatype = "INT2U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
rm(ras_1km_bin)
# ------------ #
# Create aggregate versions of the shared layer for different grains
# Aggregate to 1km
ras_2km <- terra::aggregate(ras_1km, fact = 2, fun = "mean", cores = 7)
# Write output layers
terra::writeRaster(ras_2km, paste0(path_output,"/Gridded/ReferenceGrid_Europe_frac_2000m.tif"),
datatype = "INT4U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
ras_2km_bin <- ras_2km
ras_2km_bin[ras_2km_bin>0] <- 1
terra::writeRaster(ras_2km_bin, paste0(path_output,"Gridded/ReferenceGrid_Europe_bin_2000m.tif"),
datatype = "INT2U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
# --- #
# Aggregate to 5km
ras_5km <- terra::aggregate(ras_1km, fact = 5, fun = "mean", cores = 7)
# Write output layers
terra::writeRaster(ras_5km, paste0(path_output,"/ReferenceGrid_Europe_frac_5000m.tif"),
terra::writeRaster(ras_5km, paste0(path_output,"/Gridded/ReferenceGrid_Europe_frac_5000m.tif"),
datatype = "INT4U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
ras_5km_bin <- ras_5km
ras_5km_bin[ras_5km_bin>0] <- 1
terra::writeRaster(ras_5km_bin, paste0(path_output,"ReferenceGrid_Europe_bin_5000m.tif"),
terra::writeRaster(ras_5km_bin, paste0(path_output,"Gridded/ReferenceGrid_Europe_bin_5000m.tif"),
datatype = "INT2U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
rm(ras_5km_bin); gc()
# --- #
# Aggregate to 10km
ras_10km <- terra::aggregate(ras_1km, fact = 10, fun = "mean", cores = 7)
# Write output layers
terra::writeRaster(ras_10km, paste0(path_output,"/ReferenceGrid_Europe_frac_10000m.tif"),
terra::writeRaster(ras_10km, paste0(path_output,"Gridded/ReferenceGrid_Europe_frac_10000m.tif"),
datatype = "INT4U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
ras_10km_bin <- ras_10km
ras_10km_bin[ras_10km_bin>0] <- 1
terra::writeRaster(ras_10km_bin, paste0(path_output,"ReferenceGrid_Europe_bin_10000m.tif"),
terra::writeRaster(ras_10km_bin, paste0(path_output,"Gridded/ReferenceGrid_Europe_bin_10000m.tif"),
datatype = "INT2U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
# --- #
# Aggregate to 20km
ras_20km <- terra::aggregate(ras_1km, fact = 20, fun = "mean", cores = 7)
# Write output layers
terra::writeRaster(ras_20km, paste0(path_output,"/Gridded/ReferenceGrid_Europe_frac_20000m.tif"),
datatype = "INT4U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
ras_20km_bin <- ras_20km
ras_20km_bin[ras_20km_bin>0] <- 1
terra::writeRaster(ras_20km_bin, paste0(path_output,"Gridded/ReferenceGrid_Europe_bin_20000m.tif"),
datatype = "INT2U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
rm(ras_20km_bin); gc()
# --- #
# Aggregate to 50km
ras_50km <- terra::aggregate(ras_1km, fact = 50, fun = "mean", cores = 7)
# Write output layers
terra::writeRaster(ras_50km, paste0(path_output,"/ReferenceGrid_Europe_frac_50000m.tif"),
terra::writeRaster(ras_50km, paste0(path_output,"/Gridded/ReferenceGrid_Europe_frac_50000m.tif"),
datatype = "INT4U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
ras_50km_bin <- ras_50km
ras_50km_bin[ras_50km_bin>0] <- 1
terra::writeRaster(ras_50km_bin, paste0(path_output,"ReferenceGrid_Europe_bin_50000m.tif"),
terra::writeRaster(ras_50km_bin, paste0(path_output,"Gridded/ReferenceGrid_Europe_bin_50000m.tif"),
datatype = "INT2U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
# --- #
# Aggregate to 100km
ras_100km <- terra::aggregate(ras_1km, fact = 100, fun = "mean", cores = 7)
# Write output layers
terra::writeRaster(ras_100km, paste0(path_output,"/Gridded/ReferenceGrid_Europe_frac_100000m.tif"),
datatype = "INT4U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
ras_100km_bin <- ras_100km
ras_100km_bin[ras_100km_bin>0] <- 1
terra::writeRaster(ras_100km_bin, paste0(path_output,"Gridded/ReferenceGrid_Europe_bin_100000m.tif"),
datatype = "INT2U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
......