File:Gdd5 lgm chelsa ccsm4 envirem france 1.png

1. 3. "R" script: calculate gdd5, pet and ptopet rasters with Envirem package
   4. default: area of france
   5. with envirem
   6. source: chelsa ccsm4 lgm files
   7. version 0001b
   8. 9.10.2010

library(raster) library(rgeos) library(maps) library(viridis) library(envirem) library(base) library(stringr) library(gtools)

area_1=1 stage_1=1 stage_2=1 stage_3=1 stage_4=1 rsaaga_1=1

1. 1. france

if (area_1==1) { extent1 <- extent(-10,10,40,50) poly <- readWKT("POLYGON((

  -10 50,
  10 50,
  10 40,
  -10 40,
  -10 50

))", p4s = CRS("+proj=longlat +datum=WGS84"))

}

1. 1. europe

if (area_1==2) { extent1 <- extent(-20,80,30,70)

poly <- readWKT("POLYGON((

 -20 70,
  80 70,
  80 30,
 -20 30,
 -20 70

))", p4s = CRS("+proj=longlat +datum=WGS84"))

}

inputDir <- "./climate" outputDir <- "./snapshot" processpath2<-"./snapshot2" outdir1 <- './envirem1' dempath1<-"D:/razter1/high_longlat.tif"

name1='./envirem1/test_growingDegDays5.tif' name2='./envirem1/test_annualPET.tif'

outname1="gdd5.nc" outname2="pet.nc" outname3="ptopet.nc" outname4="precipann.nc" outdpath1="./twi1.nc"

processpath1<-outputDir indir1<-processpath2

Sys.which('gdalinfo') Sys.which('gdal_translate')

list.files(inputDir)

dir.create(outputDir) dir.create(processpath2)

1. 1. creation selection ractangle

makerasters=1 calculate_average_temp=1

if(stage_1==1) {

1. plot(poly, border = 'blue', xpd = NA, add = TRUE)
2. crop and mask input rasters

files <- list.files(inputDir, pattern = '.tif$', full.names = TRUE)

1. read in a precipitation raster and crop/mask and keep as a mask layer
   1. precipitation rasters have the term "_pr_" in the file name
   2. crop/mask using polygon

grep(files, pattern = '_prec_', value=TRUE)[1]

precipRaster <- raster(grep(files, pattern = '_prec_', value=TRUE)[1])

1. precipRaster <- crop(precipRaster, poly)

precipRaster <- mask(precipRaster, poly)

precipRaster <- crop(precipRaster, extent1)

1. precipRaster <- mask(precipRaster, extent1)

1. create a terrestrial/marine mask by setting values < 0 to NA, and values >= 0 to 1

precipRaster[precipRaster < 0] <- NA precipRaster[!is.na(precipRaster)] <- 1

plot(precipRaster, col = 'blue', legend = FALSE) title(main = 'terrestrial mask')

1. supply some additional flags to writeRaster() to employ compression to achieve smaller file sizes

tifOptions <- c("COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=6")

for (i in 1:length(files)) {

 cat(i, ' ')
 r <- raster(files[i])
 #r <- crop(r, poly)
 r <- crop(r, extent1)
 # apply terrestrial mask
 r <- mask(r, precipRaster)
 
 # if temperature, divide values by 10
 ## recognize temperature files as those containing the terms tasmin or tasmax
 if (grepl('tmin|tmax', files[i])) {

# print("t ")

   r <- r / 10
   r <- r-273.15
 }
 
 r <- aggregate(r, fact = 2)
 
 outfile <- paste0(outputDir, '/', names(r), '.tif')
 writeRaster(r, filename = outfile, format = 'GTiff', options = tifOptions, overwrite = TRUE)

}

}

if(stage_2==1) {

old_files <- list.files(processpath1, pattern = "*.tif", full.names = TRUE)

1. old_files

old_files_2=mixedsort(old_files)

len1=length(old_files)

print (len1)

1. old_files_2

precipfiles=grep('prec', old_files_2, value=TRUE)

print (precipfiles)

len2=length(precipfiles)

for (n in 1:len2) { infilee=precipfiles[n] #print (infilee) a=sprintf("%02d",n) #print (a) newfileename=paste0(processpath2, "/test_precip_",a,".tif") print (newfileename) file.copy(from = infilee, to = newfileename)

}

tminfiles=grep('tmin', old_files_2, value=TRUE)

print (tminfiles)

len2=length(tminfiles)

for (n in 1:len2) { infilee=tminfiles[n] #print (infilee) a=sprintf("%02d",n) #print (a) newfileename=paste0(processpath2,"/test_tmin_",a,".tif") print (newfileename) file.copy(from = infilee, to = newfileename) }

tmaxfiles=grep('tmax', old_files_2, value=TRUE)

print (tmaxfiles)

len2=length(tmaxfiles)

for (n in 1:len2) { infilee=tmaxfiles[n] #print (infilee) a=sprintf("%02d",n) #print (a) newfileename=paste0(processpath2,"/test_tmax_",a,".tif") print (newfileename) file.copy(from = infilee, to = newfileename) }

rasterTemplate <- raster('./snapshot2/test_precip_01.tif')

ETsolradRasters(rasterTemplate = rasterTemplate, year = -25000, outputDir = processpath2, overwrite = TRUE)

}

if(stage_3==1) {

old_files <- list.files(processpath1, pattern = "*.tif", full.names = TRUE)

1. old_files

old_files_2=mixedsort(old_files)

tmaxfiles=grep('tmin', old_files_2, value=TRUE)

print (tmaxfiles)

len2=length(tmaxfiles)

print("Tmean ..")

for (n in 1:len2) { a=sprintf("%02d",n) infileename1=paste0(processpath2,"/test_tmax_",a,".tif") infileename2=paste0(processpath2,"/test_tmin_",a,".tif") outfileename1=paste0(processpath2,"/test_tmean_",a,".tif") #print (n) print (outfileename1) r1=raster(infileename1) #print("RK") r2=raster(infileename2)

#plot(r2) #print("R3 ...") r3=(r1+r2)/2

#plot(r3) #print("WR ...") writeRaster(r3, filename = outfileename1, format = 'GTiff',overwrite = TRUE)

#file.copy(from = infilee, to = newfileename) }

}

if(stage_4==1) {

dir.create(outdir1)

assignNames(reset = TRUE)

assignNames(tmax = "test_tmax_##",

          tmin = "test_tmin_##",
          tmean="test_tmean_##",
          precip = "test_precip_##",
          solrad = "et_solrad_##"
          )

chelsaFiles <- list.files(processpath2, pattern = 'test', full.names = TRUE)

chelsaStack <- stack(chelsaFiles) solarFiles <- list.files(processpath2, pattern = 'solrad', full.names = TRUE) solarStack <- stack(solarFiles) print(solarFiles)

verifyFileStructure(processpath2, returnFileNames = FALSE)

verifyRasterNames(chelsaStack, solradstack = solarStack)

varnames()

inputDir <- indir1 outputDir <- outdir1

print ("Generating rasters") generateRasters(var = 'growingDegDays5',prefix = 'test_',maindir = indir1,outputDir = outdir1)

1. print("PET")

generateRasters(var = 'annualPET',prefix = 'test_',maindir = indir1,outputDir = outdir1)

1. read in the resulting raster

result <- raster(name1) plot(result, col = inferno(100)) writeRaster(result, filename = outname1, format = 'netCDF',overwrite = TRUE)

result <- raster(name2) plot(result, col = inferno(100)) writeRaster(result, filename = outname2, format = 'netCDF',overwrite = TRUE)

print(" Prec ..." ) rasterFiles <- list.files(indir1, pattern = '.tif$', full.names = TRUE)

env <- stack(rasterFiles) precip <- grep('prec', names(env), value=TRUE) precip <- stack(envprecip)

annualprecip<-sum(precip)/10

pet=raster(outname2)

ptopet=annualprecip/pet

plot(annualprecip, col = inferno(100)) plot(ptopet, col = inferno(100))

writeRaster(annualprecip, filename = outname4, format = 'netCDF',overwrite = TRUE)

writeRaster(ptopet, filename = outname3, format = 'netCDF',overwrite = TRUE)

}

if (rsaaga_1==1) {

e1 <- extent1

elev0<-raster(dempath1)

alt <- crop(elev0, e1)

plot(alt)

tempInfile <- 'saga_temp_twi_in' writeRaster(alt, filename=tempInfile, format='SAGA',overwrite = TRUE)

env1 <- rsaga.env()

outfile <- 'saga_temp_twi_out.sgrd' call <- rsaga.wetness.index(in.dem='saga_temp_twi_in.sgrd', out.wetness.index=outfile, env=env1)

twi1 <- raster('saga_temp_twi_out.sdat')

plot(alt)

plot(twi1)

writeRaster(twi1, filename = outdpath1, format = 'netCDF',overwrite = TRUE)

}