Commit 19950b36 authored by Brage.Forland's avatar Brage.Forland

Added missing file for WorldClimate spatial clipping.

parent 7cf86d98
###
## This script clips the worldclim dat aset to the extent of Norway. It assumes that the
## WorlClim files are already downloaded and upacked in the location
## given by inputRasterLocation below.
##
## Input files must be downloaded from:
##
## http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_10m_bio.zip
## http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_5m_bio.zip
## http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_2.5m_bio.zip
## http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_30s_bio.zip
##
library(rgdal)
library(sf)
library(fasterize)
library(raster)
# Clip WORLDClim data to the extent of a given region. Save a GTIFF of the selected layers.
#
# clipRegion A Simple Features Layer containing the clipping region
# inputRasterLocation Directory containing the already downloaded WORLDClim data set
# outputLocation Directory where the clipped output files will be saved
# climateVariables Which layers to keep from the WORLDClim data set
# resolution WORLDClim resolution, one of c('10m', '5m', '2.5m', '30s')
#
if(Sys.getenv("ARTSAPP_HOME")!="") {
ARTSAPP_HOME=Sys.getenv("ARTSAPP_HOME")
} else {
ARTSAPP_HOME=normalizePath(getwd())
}
clipRegionLoc <- paste0(ARTSAPP_HOME, "/data/GADMNorway")
clipRegionLayerName <- "NOR_adm0"
#clipRegionLoc <- paste0(ARTSAPP_HOME, "/data/norway")
#clipRegionLayerName <- "borders"
inputRasterLocation = paste0(ARTSAPP_HOME,'/WCD.tmp')
outputLocation <- paste0(ARTSAPP_HOME, "/data/WORLDCLIMData");
clipWORLDCLIMData <- function(inputRasterLocation, outputLocation, clipRegion,
resolution, climateVariables = c("01","04","12","15"),insideMaskOnly=TRUE ) {
climateLayerDescriptions <- list(
"01" = "Annual Mean Temperature (°C)",
"04" = "Annual Temperature Range (°C)",
"12" = "Annual Precipitation (mm)",
"15" = "Coefficient of Precipitation Variation (%)"
)
climateLayerNames <- paste('wc2.0_bio_', rep(resolution,length(climateVariables)),"_", climateVariables, sep="")
climateLayerFiles <- paste(inputRasterLocation, "/", climateLayerNames,".tif", sep="")
names(climateLayerFiles) <- climateVariables
if(! all(file.exists(climateLayerFiles))){
stop(paste("Could not find raster input files:",paste(climateLayerFiles, collapse=", ")))
}
if(!dir.exists(outputLocation)) {
dir.create(outputLocation, recursive=TRUE)
}
# Templategrid, raster of valid observation sites.
rasterTemplate <- raster(climateLayerFiles[1])
if(st_crs(rasterTemplate) != st_crs(clipRegion)) {
message(paste("Changing projection of clipRegion", st_crs(clipRegion)$proj4string, "->", st_crs(rasterTemplate)$proj4string));
clipRegion <- st_transform(clipRegion,st_crs(rasterTemplate))
}
rasterTemplate <- crop(rasterTemplate, clipRegion, snap="out")
# Cut away outside valid region,
# rasterMask <- fasterize(clipRegion, rasterTemplate)w
rasterMask <- raster::rasterize(clipRegion, rasterTemplate)
# but only keep points where climate data is available.
rasterMask <- raster::mask(rasterMask, rasterTemplate)
# The CRS is the same, this ensures that the string representation of the CRS are also identical.
crs(rasterMask) <- crs(rasterTemplate)
readr::write_tsv(data.frame(cbind(unlist(climateLayerNames), unlist(climateLayerDescriptions))),
paste0(outputLocation,"/layers.txt"),
col_names=FALSE)
raster::writeRaster(rasterMask, paste0(outputLocation,"/template.tif"),format="GTiff", datatype="INT1U", overwrite=TRUE)
outputFiles <- lapply(X = climateLayerFiles, FUN = function(climateLayerFile, rasterMask, outputLocation, insideMaskOnly) {
rasterLayer <- raster(climateLayerFile)
rasterLayer <- crop(rasterLayer,rasterMask)
if(insideMaskOnly) {
rasterLayer <- raster::mask(rasterLayer, rasterMask)
}
outputFileName <- basename(climateLayerFile)
writeRaster(rasterLayer,paste(outputLocation,outputFileName,sep="/"), format="GTiff", overwrite=TRUE)
outputFileName
}, rasterMask=rasterMask, outputLocation=outputLocation, insideMaskOnly=insideMaskOnly )
outputFiles
}
clipRegion <- st_read(clipRegionLoc, clipRegionLayerName)
for(resolution in c("10m","5m","2.5m","30s")) {
if(file.exists(sprintf("%s/wc2.0_bio_%s_15.tif", inputRasterLocation, resolution))) {
outputSubLocation <- paste(outputLocation, resolution, sep="/")
outputSubLocationUnmasked <- paste(outputLocation, "unmasked", resolution, sep="/")
message(paste("Clipping worldclim 2.0 data: ", resolution, "to", outputSubLocation))
clipWORLDCLIMData(inputRasterLocation, outputSubLocation, clipRegion, resolution)
clipWORLDCLIMData(inputRasterLocation, outputSubLocationUnmasked, clipRegion, resolution, insideMaskOnly=FALSE)
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment