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
109
110
111
112
113
114
115
116
117
118
119
120
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