File:Lgm grassfrac mpiesm chelsa rf 1.png

2. 1. Climate data random forest downscaler
4. CMIP5 MPIESM LGM output and CHELSA MPIESM data
5. 1. 27.5.2021 v 0000.0002

library(raster) library(rgdal) library(ncdf4) library(tictoc) library(randomForest)

1. library(gbm)
2. library(caret)

install_libraries=FALSE

if(install_libraries==TRUE)

{

options(repos = "https://ftp.acc.umu.se/mirror/CRAN/")

install.packages("raster")
install.packages("rgdal")
install.packages("ncdf4")
install.packages("tictoc")
install.packages("randomForest")

}

writeout<-function(oras, outn, varnamex, varunitx, longnamex) {

crs(oras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(oras, filename=outn, overwrite=TRUE, format="CDF", varname=varnamex, varunit=varunitx, longname=longnamex, xname="lon", yname="lat")

}

1. 1. snow

load_cmip5_variable <- function(dataname1, invarname1, posit, numyears, month1) {

print("Loading cmip5 data ...")

nppin1 <- nc_open(dataname1)

vext1<-c(0,360,-90,90)

lok1=posit*12+month1 mara=numyears*12

stacksnow1<-stack()

for(n in 1:mara) { # print (".") snow00 <- ncvar_get( nppin1, varid=invarname1,start=c(1,1,lok1), count=c(-1,-1,1) )

snow01=t(snow00) snow02<-apply(snow01,2,rev)

snow0=raster(snow02)

extent(snow0)<-vext1 names(snow0)<-invarname1

snow2=rotate(snow0)

crs(snow2) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

stacksnow1 <- stack( stacksnow1 , snow2 )

lok1=lok1+12

}

   rasnow0<-mean(stacksnow1)
   
    print (rasnow0)
  
   rasnow1=rasnow0
   
    rasnow1[is.na(rasnow1)] <- 0

return(rasnow1)

1. 1. loadipslnpp

 }

create_small_stack<-function(small1, stack1) { print("Create small stack ...")

mara1=nlayers(stack1)

print(mara1)

sasa1<-stack()

sasa1<-stack(sasa1, small1)

for(n in 1:mara1) { print(n) r1<-subset(stack1,n) s1<-resample(r1,small1) plot(s1) sasa1<-stack(sasa1, s1)

}

return(sasa1) }

create_big_stack_variables<-function(inlist1, invarlist1, rext1, xsiz1, ysiz1) { print("Loading big raster stack ...")

 	mara1=length(inlist1)
 	
 	print (mara1)
 	 
 	instak1<-stack()
 	sabluna1<-raster(ncol=xsiz1, nrow=ysiz1)

for(n in 1:mara1) {

rasname1<-inlist1[n] rasvarname1<-invarlist1[n] inras0<-raster(rasname1) inras1<-crop(inras0, rext1) extent(sabluna1)<-extent(inras1) inras2<-resample(inras1, sabluna1) plot(inras2) names(inras2)<-rasvarname1 instak1<-stack(instak1,inras2)

}

    ## note remove NA
#    dstak1[is.na(dstak1)] <- 0
    
    
 	 return(instak1)
 	  	

}

1. 3. program init

infilname1<-"./indata1/grassFrac_Lmon_MPI-ESM-P_lgm_r1i1p1_185001-194912.nc" invarname1<-"grassFrac"

1. infilenamelist2<-c("./indata1/CHELSA_PMIP_MPI-ESM-P_BIO_12.tif",
2. "./indata1/CHELSA_PMIP_MPI-ESM-P_BIO_01.tif",
3. "./indata1/CHELSA_PMIP_MPI-ESM-P_BIO_10.tif",
4. "./indata1/CHELSA_PMIP_MPI-ESM-P_BIO_18.tif")

1. invarnamelist2<-c("PrecipAnn",
2. "TempAnn",
3. "PrecipWarm",
4. "TempWarm")

infilenamelist2<-c("./indata1/CHELSA_PMIP_MPI-ESM-P_BIO_12.tif", "./indata1/CHELSA_PMIP_MPI-ESM-P_BIO_01.tif")

invarnamelist2<-c("PrecipAnn",

   "TempAnn")
   

targetlist2<-c("Target",invarnamelist2)

posit=0 numyears=8 month1=7 rext1<-c(-180,-120,50,80) # Beringia

1. rext1<-c(-15,60,30,70) # europe0

xsize1=750*3 ysize1=400*3

1. 1. program run

instak1<-create_big_stack_variables(infilenamelist2, invarnamelist2, rext1, xsize1, ysize1) smallraw0<-load_cmip5_variable(infilname1, invarname1, posit, numyears, month1)

small1<-crop(smallraw0,rext1)

1. remove na

small1[is.na(small1)] <- 0

1. plot(instak1)

plot(small1)

sasta1<-create_small_stack(small1, instak1)

1. names(sasta1)<-c("Target","PrecipAnn", "TempAnn")

names(sasta1)<-targetlist2

smallframe1=as.data.frame(sasta1) bigframe1=as.data.frame(instak1)

head(smallframe1)

head(bigframe1)

print("Fit and predict ...")

tic()

1. raf1 <- randomForest(Target ~ PrecipAnn+TempAnn, data = smallframe1)

raf1 <- randomForest(Target ~ ., data = smallframe1)

red1 <- predict(raf1, bigframe1)

1. raf1 <- gbm(smallframe1$Target ~ ., data = smallframe1)

1. red1 <- gbm(raf1, bigframe1)

toc()

1. print (str(red1))

1. print(head(red1))

odak1<-as.matrix(red1)

1. print (head(odak1))

1. print (length(odak1))

oras1<- raster(nrow=ysize1, ncol=xsize1)

oras1 <- setValues(oras1,odak1)

extent(oras1)<-rext1

1. plot(oras1)

print("Write output ...")

writeout(oras1,"./grass.nc","grassFrac (LGM MPI-ESM)", "MPI-ESM", "grassFrac, Last Glacial Maximum")

print(".")