Outputs/globalgrid_wgs84_1km.tif 0 → 100644
View file @ 1839f8d0
File added
Scripts/00_Functions.R 0 → 100644
View file @ 1839f8d0
#' Unify extent of a rasterstack
max_extent <- function(rl){
# given list of rasters
# returns union of extent
xmin=min(sapply(rl,function(x){extent(x)@xmin}))
xmax=max(sapply(rl,function(x){extent(x)@xmax}))
ymin=min(sapply(rl,function(x){extent(x)@ymin}))
ymax=max(sapply(rl,function(x){extent(x)@ymax}))
extent(c(xmin,xmax,ymin,ymax))
}
# Saving function
writeGeoTiff <- function(file, fname,dt = "INT2S"){
require(assertthat)
if(!has_extension(fname,"tif")) fname <- paste0(fname,".tif")
writeRaster(file,fname,
format='GTiff',
datatype = dt,
NAflag = -9999,
options=c("COMPRESS=DEFLATE","PREDICTOR=2","ZLEVEL=9"),
overwrite= TRUE )
}
emptyraster <- function(x, ...) { # add name, filename,
emptyraster <- raster(nrows=nrow(x), ncols=ncol(x),
crs=x@crs,
ext=extent(x), ...)
return(emptyraster)
}
#' Align raster data by bringing it in the same geometry and extend.
alignRasters <- function(data, template, method = "bilinear",func = mean,cl = T,expand = TRUE){
lib <- c("raster")
sapply(lib, function(...) stopifnot(require(..., character.only = T)))
if(cl) beginCluster(parallel::detectCores()-1)
if(projection(data) == projection(template)){
data <- raster::crop(data, template, snap = "near")
if(class(template) == "RasterLayer"){
if(data@ncols / template@ncols >= 2){
factor <- floor(data@ncols/template@ncols)
data <- raster::aggregate(data, fact = factor, fun = func,na.rm = TRUE,
expand=expand)
}
data <- raster::resample(data, template, method = method)
}
} else {
data <- projectRaster(data, template, method = method)
}
if(cl) endCluster()
return(data)
}
split_raster <- function(infile,outpath,parts=3) {
## parts = division applied to each side of raster, i.e. parts = 2 gives 4 tiles, 3 gives 9, etc.
lib <- c("gdalUtils","parallel","reshape2")
sapply(lib, function(...) stopifnot(require(..., character.only = T)))
# Check if infile exists
stopifnot(file.exists(infile))
filename <- strsplit(basename(infile),".tif")[[1]]
# Get the dimensions of the file
dims <- as.numeric(
strsplit(gsub('Size is|\\s+', '', grep('Size is', gdalinfo(infile), value=TRUE)),
',')[[1]]
)
# Generate window frames
xy <- list()
# t is nr. of iterations per side
t <- parts - 1
for (i in 0:t) {
for (j in 0:t) {
# [-srcwin xoff yoff xsize ysize] src_dataset dst_dataset
srcwin <- paste(i * dims[1]/parts, j * dims[2]/parts, dims[1]/parts, dims[2]/parts)
xy[[paste0(i,"_",j)]] <- data.frame(infile = infile,srcwin = srcwin, file = paste0(outpath,"/",filename,"_",i,"_",j,".tif"))
}
}
df <- melt(xy)
# Then process per src_win
cat("Start splitting: ",filename)
# Create a function to split the raster using gdalUtils::gdal_translate
split <- function(input, outfile, srcwin) {
gdal_translate(input, outfile, srcwin=srcwin)
}
# Make a copy for export
df_org <- df
# Kick out files already existing
df <- df[which(!file.exists(as.character(df$file))),]
# Make cluster
cl <- makeCluster(detectCores()-1)
clusterExport(cl, c('split', 'df'))
clusterEvalQ(cl,library(gdalUtils))
system.time({
parLapply(cl, seq_len(nrow(df)), function(i) {
split(df$infile[i], df$file[i], df$srcwin[i])
})
})
stopCluster(cl)
cat("\n")
cat("Done")
return(df_org)
}
Scripts/01_PreparePolygons.R → Scripts/01_PreparePolygonFiles.R
View file @ 1839f8d0
library(sf) library(sf)
library(terra) library(tidyverse)
library(ncdf4)
library(dplyr)
library(stringr) library(stringr)
source("Scripts/00_Functions.R") # Some Convenience functions
# Paths # Paths
path_input <- "Inputs/" path_input <- "Inputs/" # Note that inputs and other files need to be downloaded separately
path_output <- "Outputs/" path_output <- "Outputs/"
path_scripts <- "Scripts/" 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"
crs.mollweide <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" # For global equal area projections
p <- read_sf("/mnt/pdrive/bec_data/100_rawdata/GADM_AdminBoundaries/gadm_410.gpkg",quiet = T) # ------------- #
#### Prepare European countries layer ####
# Get all countries on the European continent and subset # Specifically take the various European countries and
co <- p |> sf::st_drop_geometry() |> filter(CONTINENT == "Europe") |> select(SOVEREIGN) |> collect() |> distinct()
# Now select all countries with this information # Now select all countries with this information
sf::st_drop_geometry() |> filter(CONTINENT == "Europe")
sub <- p |> filter(COUNTRY %in% co$SOVEREIGN) sub <- p |> filter(COUNTRY %in% co$SOVEREIGN)
write_sf(sub, "/home/martin/GADM_Europesubset.gpkg") write_sf(sub, "/home/martin/GADM_Europesubset.gpkg")
# --- # # --- #
# Download gridded Reference from CHELSA
chelsa_reference <- str_trim(readLines(paste0(path_input,"envidatS3paths.txt")))
curl::curl_download(
chelsa_reference,
paste0(path_input, basename(chelsa_reference) )
)
# Load the object
ras <- rast( paste0(path_input, basename(chelsa_reference) ) )
# --- # # --- #
# Download European reference grids # Download European reference grids
......
Scripts/02_PrepareGlobalRasterLayers.R 0 → 100644
View file @ 1839f8d0
library(sf)
library(terra)
library(raster)
library(ncdf4)
library(gdalUtils)
library(dplyr)
library(stringr)
source("Scripts/00_Functions.R") # Some Convenience functions
# 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"
crs.mollweide <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" # For global equal area projections
# Slightly smaller x because of dateline issues. Should not be an issue for NaturaConnect
global_extent <- c(xmin = -179.999999, xmax = 180, ymin = -90, ymax = 90)
# Temporary folder for storing intermediate products
tempdir <- '/media/martin/data/tmp'
dir.create(tempdir, showWarnings = FALSE)
# -------------------------- #
#### Prepare a global layer at 1km in WGS 1984 ####
# Rationale is that we have one consistent global baseline grid that we can use for regions even beyond
# NUTS or European jurisdiction
# We use GADM as source owing that all European polygon and district layers have issues with shorelines and
# often leave out smaller islands.
# Source: https://gadm.org/ # Stored separately here
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
co <- gadm |>
dplyr::filter(NAME_0 != 'Antarctica') %>%
st_cast("MULTIPOLYGON") %>%
st_transform(crs = crs.mollweide) %>%
mutate( dummy = 1 )
# First create a global 1deg lat-long polygon
temp_latlong <- raster(xmn = global_extent[1], xmx = global_extent[2],
ymn = global_extent[3], ymx = global_extent[4],
res = 1)
temp_latlong_fname <- file.path(tempdir,"latlong.tif")
writeGeoTiff(temp_latlong, temp_latlong_fname)
# Project the grid to global mollweide
temp_moll_fname <- file.path(tempdir,"moll.tif")
# Project to get the world extent to the Mollweide projection
gdalwarp(srcfile = temp_latlong_fname,
dstfile = temp_moll_fname,
s_srs = crs.latlong,
t_srs = crs.mollweide,
r = "near",
co = c("COMPRESS=DEFLATE","PREDICTOR=2","ZLEVEL=9"),
multi = TRUE,
overwrite = TRUE, verbose=TRUE
)
global_extent_moll <- extent(raster(temp_moll_fname))
# Now to create the global template raster appropriate at full latitude - longitudes
global_raster <- raster(ext = global_extent_moll,
crs = crs.mollweide,
resolution = 1000, # Global 1km
vals = 0
)
# Clean up
file.remove(temp_latlong_fname); file.remove(temp_moll_fname)
# ------------------- #
# Now create a split of this layer
dir.create( paste0(tempdir, '/temp_splits') )
writeGeoTiff(global_raster, paste0(tempdir, '/temp_splits/layer.tif') )
ss <- split_raster(paste0(tempdir, '/temp_splits/layer.tif'),
outpath = paste0(tempdir, '/temp_splits'),
2)
for(r in 1:nrow(ss)){
print(r)
ref <- raster(as.character( ss$file )[r])
global_raster <- fasterize::fasterize(sf = co, raster = ref)
writeGeoTiff(global_raster, fname = paste0(tempdir, '/temp_splits/world_',r,'.tif')) # Save outputs and clear
rm(global_raster)
}
# Now mosaic those layers together
gdalUtils::gdalbuildvrt(gdalfile = list.files(paste0(tempdir, '/temp_splits'), "world",full.names = T),
output.vrt = paste0(tempdir, '/temp_splits/out.vrt'),
separate = FALSE)
global_raster <- raster(paste0(tempdir, '/temp_splits/out.vrt'))
writeGeoTiff(global_raster, fname = paste0(path_output, "/globalgrid_mollweide_1km.tif")) # Save outputs and clear
# Load and reproject all layers to WGS84 and save again
# 1km
gdalwarp(srcfile = paste0(path_output, "/globalgrid_mollweide_1km.tif"),
dstfile = paste0(path_output, "/globalgrid_wgs84_1km.tif"),
s_srs = s.crs,
t_srs = latlong.crs,
r = "near",
co = c("COMPRESS=DEFLATE","PREDICTOR=2","ZLEVEL=9"),
output_Raster=TRUE, multi = TRUE,
overwrite=TRUE,verbose=TRUE)
lapply(list.files(paste0(tempdir, 'temp_splits')), file.remove) # clean up
unlink(paste0(tempdir, 'temp_splits'))
# ---------- #
ras1 <- raster("../../Projects/naturemap/data/globalgrid_mollweide_1km.tif")
ras1[!is.na(ras1)] <- 1
writeGeoTiff(ras1, fname = paste0(path_output, "/globalgrid_mollweide_1km.tif")) # Save outputs and clear
ras2 <- raster("../../Projects/naturemap/data/globalgrid_wgs84_1km.tif")
ras2[!is.na(ras2)] <- 1
writeGeoTiff(ras2, fname = paste0(path_output, "/globalgrid_wgs84_1km.tif")) # Save outputs and clear
Scripts/03_PrepareEuropeanGrid.R 0 → 100644
View file @ 1839f8d0
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)