Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
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
library(sf)
library(terra)
terra::terraOptions(progress = TRUE)
# Paths
path_input <- "Inputs/" # Note that inputs and other files need to be downloaded separately
path_output <- "Outputs/"
path_scripts <- "Scripts/"
# CRS proj4strings
crs.latlong <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" # Latitude longitude
crs.laea <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs"
# ---------- #
#### 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"))
# Load GADM dataset for global shorelines
gadm <- read_sf("/mnt/pdrive/bec_data/100_rawdata/GADM_AdminBoundaries/gadm_410.gpkg", quiet = T) |>
# Get all countries
# We exclude Antarctica and thus the French Southern and Antarctic Lands
dplyr::filter(NAME_0 != 'Antarctica') %>%
sf::st_cast("MULTIPOLYGON") %>%
sf::st_transform(crs = crs.laea) %>%
dplyr::mutate( dummy = 1 )
# 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()
# Make a template raster at 100m for Europe
template <- rast(extent = ext,
crs = crs.laea, resolution = c(100,100))
# Take at minimum the extent from the coastline layer and rasterize it to the template
ras <- terra::rasterize(coast, template, touches = TRUE)
# Also rasterize the GADM data and use it whenever the coastline grid indicates 0
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"),
datatype = "INT1U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
rm(ras2); gc()
# ----- #
ras <- rast( paste0(path_output,"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)
# Multiply with 10000 and round to save as Integer
ras_1km <- round(ras_1km * 10000)
# Write the output as fractional layer
terra::writeRaster(ras_1km, paste0(path_output,"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"),
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 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"),
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"),
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"),
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"),
datatype = "INT2U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)
# 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"),
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"),
datatype = "INT2U", gdal = c("COMPRESS=DEFLATE", "TFW=YES"), overwrite = TRUE)