File:Europe rainfall trace21ka 21000BP annual.svg

1. 1. Trace21ka Ccsm3 Glac1D Etopo1 Rainfall "variable" downscaler
2. zulu basis rain shadows for Europe
3. 1. Language "R" 3.6.1 2019
   2. tested under Windows 10 Rscript
   3. v 0006

install_libraries=FALSE

if(install_libraries==TRUE) {

install.packages("raster")
install.packages("rgdal")
install.packages("sp")
install.packages("spatialEco")
install.packages("ncdf4")
install.packages("dissever")
install.packages("viridis")
install.packages("dplyr")
install.packages("lattice")
install.packages("RColorBrewer")
install.packages("rgeos")
install.packages("sp")
install.packages("reshape2")
install.packages("data.table")
install.packages("stringr")
install.packages("rlist")
install.packages("pipeR")
install.packages("maptools")
install.packages("gdata", dependencies=TRUE)
install.packages("abind")
install.packages("Cairo")
install.packages("pals")
install.packages("REdaS")
install.packages("easyNCDF")
install.packages("numbers")
install.packages("rasterVis")
install.packages("OceanView")
install.packages("rainfarmr")

}

library(raster) library(rgdal) library(ncdf4) library(lattice) library(maptools) library(rgeos) library(spatialEco) library(dissever) library(rainfarmr)

library(RColorBrewer) library(viridis) library(pals) library(data.table) library(stringr) library(rlist) library(pipeR) library(rasterVis)

1. library(OceanView)

library(sp) library(reshape2)

library(dplyr) library(REdaS) library(easyNCDF) library(numbers)

1. library(gdata)

library(abind)

library(Cairo)

1. 1. NOTE SET THESE

1. 1. set orografy creating off

1. argreading=0 ## WARNING NOT IMPLEMENTED default 0
2. kalk_distance=0 # for speed 0, if u not neet to calc distance
3. kalk_tables=0 # suppress erroris from tables, id kalk_distance=0

1. 1. must be 1, if you create output foor certain year, first

1. 1. Control variables

arg_reading =0

orography_preprocess=0

1. 2. output area type
   3. 0: normal local output
   4. 1: global output, according to Glac!D
   5. NOT OK 2: custom orography, local NOK

1. 1. acquire trace21ka variable data, default 1

capture_variable=0

1. 1. capture **MONTHLY** wind and rain

capture_wind_rain_variables=1

1. 1. acquire glac1d elevation data,
   2. default 1 ( Tarasov glac1d)
   3. 2 Peltier 1ce 6G-D 122k

capture_elevation=0

1. 1. preprocess rain ds variable rasters

capture_rain_orography=0

dsobject=1 # 0 current time 1 ice age

calculate_distance=1 # 1 slow calculate distance raster

1. 1. monthly downscaling

rain_monthly_downscaling=1

1. 1. normal one year etopo1 single area downscaling

normal_ds=0

1. 1. rainfarm test utility

run_rainfarm=9

stage2=0 ## 0: no "stage 2", 1: "stage 2"

1. 1. either global_ds or normal_ds must set 1
   2. global area downscaling to glac1d

global_ds=0

1. ! afterburner accurate "kustomorofilee1.nc" downscaling, default 0
   1. need accurate orog file (GTOPO30 slice, SRTM etc.)

accurate_ds=0

1. kustomorofilee1="./predata/europe30.nc"
   1. plotting is under alpha stage!

plot_result=1

1. 1. downscaling method from orography
   2. 0 delta, 1 spatialeco, 2 dissever, 3 temperature delta lapse rate 6,5 c/ 1 km
   3. method 2 slow

method_ds_oro=1

1. 1. downscaling method for temperature

method_ds_var=3

kustomorofilee1="../atrace1/predata/dordogne.nc"

1. 1. year to render
2. fage=14559

fage=21000

1. 1. number of years to average, since "fage"

numyears=1

1. 1. month of year to render

fmonth=1

1. 1. sea level, below current: eq 120 means height of -120 to current
   2. not needed set by default, calculated from age w/ thisscript

sealevel=-110

1. 1. area to render, lon1, lon2 lat1 lat2 rectangle

1. 1. europe

lon1=-15.0 lon2=40.0 lat1=30.0 lat2=70.0

1. 1. NOTE : use longitude 200 instead of -160
2. beringia
3. lon1=160
4. lon2=240
5. lat1=50.0
6. lat2=80.0

1. 1. global
2. klopaali1=1

1. lon1=-180.0
2. lon2=180.0
3. lat1=-90.0
4. lat2=90.0

1. 1. dordogne 1

1. lon1=-2.0
2. lon2=4.0
3. lat1=42.0
4. lat2=48.0

1. 1. dordogne 2
2. lon1=0.0
3. lon2=2.0
4. lat1=44.0
5. lat2=45.5

1. 1. dordogne 3
2. lon1=-1.0
3. lon2=2.0
4. lat1=44.0
5. lat2=45.5

filename1="KKK" tslocation1=-999 elevlocation1=-999

variaabeli="RAIN" unitti="mm" lonkkunimi="Rain" lonvar1="lon" latvar1="lat"

1. 1. input data dirs
   2. etopo1 dir

predata="./predata/"

1. 1. Trace21ka TS dir

predata2="./predata2/"

datastore1="d:/datav3/"

dir.create("./data_processing") dir.create("./plotz") dir.create("./indata_wind/")

dir.create(predata)

1. outpath11<-"./data_processing/area.nc"
2. outpath12<-"./data_processing/area_neworog.nc"
   1. glac1D HDC file

inpath_tarasov1<-"../atrace1/predata/TOPicemsk.GLACD26kN9894GE90227A6005GGrBgic.nc" inpath_peltier1<-"../atrace1/predata/IceT.I6F_C.131QB_VM5a_1deg.nc"

1. outpath22<-"./data_processing/oroin.nc"
   1. ! grid-registreted etopo1
2. inpath11<-"./predata/ETOPO1_Ice_g_gdal.grd"
   1. cell reg maybe inaccurate, but func w/dissever
3. inpath_etopo1<-"../trace1/predata/ETOPO1_Ice_c_gdal.grd"

inpath_etopo1<-"D:/datavarasto/etopo1.nc"

1. inpath_etopo1<-"./predata/ETOPO1_Ice_g_geotiff.tif"

inpath_etopo2<- "./predata/etopo1_360.nc" data_processing="./data_processing/"

windpath1="./indata_wind/windirs1.nc"

1. 1. wind, rain vars 12 montsh

windpath2="./indata_wind/windirs2.nc" rainpath2="./indata_rain/rain12.nc"

1. 1. out rain downscaled rain 12

outrainpath2="./indata_out/rain12_ds1.nc" outrainpath3="./indata_out/raina_ds1.nc"

rainfarm_name1="./indata_out/rainfarm.nc"

invarfname1<-"./data_processing/tempin.nc" inorofname2<-"./data_processing/area.nc" inorofname3<-"./data_processing/area_neworog.nc"

outorofname1<-"./data_processing/area_glacial.nc"

outvarfname1<-"./data_result/result.nc" outvarfname_accurate1<-"./data_result/accurate.nc"

inorofname360_1<-"./data_processing/oroin.nc" invarfname360_1<-"./data_processing/rainin.nc" inicefname360_1<-"./data_processing/maskin.nc"

inorofname180_1<-"./data_processing/oroin_180.nc" invarfname180_1<-"./data_processing/rainin_180.nc" inicefname180_1<-"./data_processing/maskin_180.nc"

plottauzdir1="./plotz/" plotfname0="result0.svg" plotfname1="./plotz/result1.svg"

resultdir1<-"./data_result/"

dir.create(data_processing) dir.create(plottauzdir1) dir.create(resultdir1)

1. 1. csv related variables, not yet used

smallrst="./data/small/" desticsv="./data/" destirst="./data/rst/" smalltemperaturepath<-"./data/small/small_temperature.nc" elevpath<-"./data/rst/elevation.nc" smallelevpath<-"./data/small/small_elevation.nc" distancepath<-"./data/rst/distance.nc" smalldistancepath<-"./data/small/small_distance.nc" latpath<-"./data/rst/latitude.nc" smalllatpath<-"./data/small/small_latitude.nc" lonpath<-"./data/rst/longitude.nc" smalllonpath<-"./data/small/small_longitude.nc"

indir1<-"D:/datavarasto/"

iname1<-"CHELSA_bio10_12.tif"

iname2="CHELSA_PMIP_CCSM4_BIO_12.tif" iname3="etopo1.nc"

if (dsobject==0) { iname4="../atrace1/data_processing/area.nc" }

if (dsobject==1) { iname4="../atrace1/data_processing/area_glacial.nc" }

1. unlink("indata_rain")
2. unlink("indata_map")
3. unlink("indata_out")
4. unlink("indata_small")

dir.create("indata_rain") dir.create("indata_map") dir.create("indata_out") dir.create("indata_small")

indirname1<-paste0(indir1, iname1) indirname2<-paste0(indir1, iname2) indirname3<-paste0(indir1, iname3) indirname4<-iname4

indirname5<-"../atrace3/data_processing/rainin_180.nc"

outname1<-"./indata_rain/chelsa_current_annual_rain.nc" outname2<-"./indata_rain/chelsa_lgm_ccsm4_annual_rain.nc"

outname3<-"./indata_map/etopo1.nc" outname4<-"./indata_map/lons.nc" outname5<-"./indata_map/lats.nc" outname6<-"./indata_map/distance.nc" outname7<-"./indata_map/slope.nc" outname8<-"./indata_map/aspect.nc" outname9<-"./indata_map/hillshade.nc" outname10<-"./indata_map/tpi.nc" outname11<-"./indata_map/westness.nc" outname12<-"./indata_map/blurelev.nc" outname13<-"./indata_map/windir.nc" outname14<-"./indata_map/etopo_big.nc"

smallrainame1<-"./indata_small/smallrain.c"

1. 1. france
2. reso1min=2.5
3. coarsereso1deg=1

1. france
2. lon1=-5
3. lat1=42
4. lon2=8
5. lat2=50

1. europe

hillshade_slope=2 hillshade_aspect=260

reso1min=6 coarsereso1degx=3.5 coarsereso1degy=3.5

lon1=-15 lat1=30 lon2=40 lat2=70

slon1=-20 slat1=25 slon2=45 slat2=75

1. rainfarm doms lat, lon

superlon1=-20 superlon2=70 superlat1=30 superlat2=90

1. superlon1=-20
2. superlon2=40
3. superlat1=30
4. superlat2=90

1. 1. icerainame1<-"../atrace3/data_processing/rainin_180.nc"

icemaskname1<-"../atrace3/data_processing/oroin_180.nc"

smallrainame1<-"./indata_small/smallrain.nc"

refrainame1<-"./indata_rain/chelsa_current_annual_rain.nc" refrainame2<-"./indata_rain/chelsa_lgm_ccsm4_annual_rain.nc" outname1<-"./indata_out/raintest.nc" outname2<-"./indata_out/raintest2.nc" origoname1<-"./indata_out/rainsmall.nc"

errorname1<-"./indata_out/errordif.nc"

iname1<-"./indata_map/etopo1.nc" iname2<-"./indata_map/distance.nc" iname3<-"./indata_map/aspect.nc" iname4<-"./indata_map/tpi.nc" iname5<-"./indata_map/lons.nc" iname6<-"./indata_map/hillshade.nc" iname7<-"./indata_map/slope.nc" iname8<-"./indata_map/windir.nc" iname9<-"./indata_map/blurelev.nc" iname10<-"./indata_map/lats.nc"

ixname1<-"./indata_map/etopo_big.nc"

1. TPI for different neighborhood size:

tpiw <- function(x, w=15) { m <- matrix(1/(w^2-1), nc=w, nr=w) m[ceiling(0.5 * length(m))] <- 0 f <- focal(x, m) x - f }

acquire_rain_orography_data<-function() {

xsiz1<-as.integer(((lon2-lon1)*60)/reso1min) ysiz1<-as.integer(((lat2-lat1)*60)/reso1min)

smallxsiz1<-as.integer((slon2-slon1)/coarsereso1degx) smallysiz1<-as.integer((slat2-slat1)/coarsereso1degy)

print (xsiz1) print (ysiz1)

print (smallxsiz1)

	print (smallysiz1) 
	         

ext1 <- extent(lon1,lon2,lat1,lat2) ext2 <- extent(slon1,slon2,slat1,slat2)

r1<-raster(indirname1) kropped10<-crop(r1,ext1) kropped11<-crop(r1,ext2)

r2<-raster(indirname2)

kropped20<-crop(r2,ext1)

r3<-raster(indirname3) kropped31<-crop(r3,ext1)

r4<-raster(indirname4) kropped41<-crop(r4,ext1)

1. 1. glacial rain

if(dsobject==1) { r5<-raster(indirname5) kropped51<-crop(r5,ext1) kropped50<-kropped51 }

kropped30=kropped31

kropped30[kropped30 < 1] <- 0

kropped40=kropped41 kropped40[kropped40 < 1] <- 0

xy <- matrix(rnorm(xsiz1*ysiz1),xsiz1,ysiz1) #image(xy)

rast0 <- raster(xy) extent(rast0) <- extent(lon1,lon2,lat1,lat2) crs(rast0) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" #plot(rast0)

xy2 <- matrix(rnorm(smallxsiz1*smallysiz1 ), smallxsiz1,smallysiz1)

small1<- raster(xy2) extent(small1) <- extent(slon1,slon2,slat1,slat2) crs(small1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" #plot(rast0)

if(dsobject==0) { smallrain1=resample(kropped11,small1) } if(dsobject==1) { smallrain1=resample(kropped51,small1) smallrain1=smallrain1*3600*24*365

}

rout1=resample(kropped10,rast0) rout2=resample(kropped20, rast0)

## curent oro ##rout3=resample(kropped30, rast0) ## glacial oro

print("Glacial orography") rout3=resample(kropped40, rast0)

#plot(kropped1)

rout2=rout2/10

crs(smallrain1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(smallrain1, filename=smallrainame1, overwrite=TRUE, format="CDF", varname="Rain", varunit="mm", longname="Rainfall_small", xname="lon", yname="lat")

crs(rout1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(rout1, filename=outname1, overwrite=TRUE, format="CDF", varname="Rain", varunit="mm", longname="Rainfall_current", xname="lon", yname="lat")

crs(rout2) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(rout2, filename=outname2, overwrite=TRUE, format="CDF", varname="Rain", varunit="mm", longname="Rainfall_LGM", xname="lon", yname="lat")

crs(rout3) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(rout3, filename=outname3, overwrite=TRUE, format="CDF", varname="Elev", varunit="m", longname="Elevation", xname="lon", yname="lat")

latras <- lonras <- rout3 xy <- coordinates(rout3) lonras[] <- xy[, 1] latras[] <- xy[, 2]

rasa1<-rout3

rasa1[rasa1 > 1] <- NA rasa1[rasa1 <0] <- 1

windir<- rout3

windir<-windir*0 windir<-windir+1 windir<-windir*4.55

#plot(rasa1) if(calculate_distance==1) { print( "calc distance slow ...") distrasa1 <- distance(rasa1) }

crs(kropped31) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(kropped31, filename=outname14, overwrite=TRUE, format="CDF", varname="Elevation2", varunit="mm", longname="Elevation2", xname="lon", yname="lat")

crs(windir) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(windir, filename=outname13, overwrite=TRUE, format="CDF", varname="Windir", varunit="radins", longname="Windir", xname="lon", yname="lat")

crs(lonras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(lonras, filename=outname4, overwrite=TRUE, format="CDF", varname="Lonx", varunit="mm", longname="Lonx", xname="lon", yname="lat")

crs(latras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(latras, filename=outname5, overwrite=TRUE, format="CDF", varname="Laty", varunit="mm", longname="Laty", xname="lon", yname="lat")

if(calculate_distance==1) { crs(distrasa1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(distrasa1, filename=outname6, overwrite=TRUE, format="CDF", varname="Dist", varunit="m", longname="Distance", xname="lon", yname="lat") }

slope <- terrain(rout3, opt='slope') aspect <- terrain(rout3, opt='aspect') hill <- hillShade(slope, aspect, hillshade_slope, hillshade_aspect) #plot(hill, col=grey(0:100/100), legend=FALSE, main='France')

tpi5 <- tpiw(rout3, w=15)

crs(slope) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(slope, filename=outname7, overwrite=TRUE, format="CDF", varname="Slope", varunit="mm", longname="Slope", xname="lon", yname="lat")

crs(aspect) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(aspect, filename=outname8, overwrite=TRUE, format="CDF", varname="Aspect", varunit="mm", longname="Aspect", xname="lon", yname="lat")

crs(hill) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(hill, filename=outname9, overwrite=TRUE, format="CDF", varname="Hillshade", varunit="m", longname="Hillshade", xname="lon", yname="lat")

crs(tpi5) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(tpi5, filename=outname10, overwrite=TRUE, format="CDF", varname="tp15", varunit="m", longname="tpi5", xname="lon", yname="lat")

etopo=rout3 etopo[etopo < 1] <- 0 etopo[is.na(etopo[])] <- 0

blurelev1 <- focal(rout3, w=matrix(1,nrow=7,ncol=7), fun=mean, pad=TRUE, na.rm = TRUE)

crs(blurelev1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(blurelev1, filename=outname12, overwrite=TRUE, format="CDF", varname="Blurelev", varunit="m", longname="Blurelev", xname="lon", yname="lat") }

1. 1. CODE BEGINS

exte1<- extent(lon1,lon2,lat1,lat2)

1. 1. global doordinates: pacific or europe centered
   2. rasnum 1: -180 - 180
   3. rasnum 2: 0 - 360

rasnum1=1

if(lon1>180)

  {

rasnum1<-2;

  }
   
   if(lon2>180)
   {

rasnum1<-2; }

sealevel_from_age<-function(inage) { #agetable1<-c(24000, 18000, 15000, 14000, 12000,11000,7000,0) #leveltable1<-c(-120, -120, -110, -80, -65, -60,-5,0 )

agetable1<-c(24000,22000,20000,18000,16000, 15000,14000,13000,12000,11000, 10000,8000,7000,6000 ,0) leveltable1<-c(-125,-130,-125,-120,-115, -110,-80,-75,-65,-60, -45,-15,-10,-5 ,0 )

1. plot(agetable1, leveltable1, main = "Sea level curve")

#points(approx(agetable1, leveltable1), col = 2, pch = "*")

# 	selev1=points(approx(agetable1, leveltable1), col = 2, pch = "*")

seali<<-approx(agetable1, leveltable1, xout=inage)

    sealevel<-as.numeric(seali)

   print("Calculated sea level ")
   print(seali)
   print(sealevel)

}

default_settings_on<-function() { print("Default settings on.") arg_reading <<-1 orography_preprocess<<-1 capture_variablee<<-1 capture_elevation<<-1 normal_ds<<-1 global_ds<<-0 accurate_ds<<-0 plot_result<<-1 }

delete_directories<-function() { print("Delete dirs") getwd() unlink(data_processing, recursive = TRUE) unlink(plottauzdir1, recursive = TRUE) unlink(resultdir1, recursive = TRUE)

}

read_argus<-function() {

if (arg_reading >0) { args = commandArgs(trailingOnly=TRUE)

if (length(args)==0) { ## stop("Not args. "=FALSE) print("usage like: rscript --vanilla tracs3.r 10000") # stop("Abort.") }

if (length(args)>0) {

print("Reading of argus done.")

eka<-args[1] toka<-args[2] fage<<-as.numeric(eka) print(fage)

delete_directories() default_settings_on() #fmonth=as.numeric(toka) } }

} ## argus

1. 1. 1. MAIN PROGE

preprocess_orography<-function() {

print("Processing orography, wait ...") rin1<-raster(inpath_etopo1)

rin1

x1 <- crop(rin1, extent(-180, 0, -90, 90)) x2 <- crop(rin1, extent(0, 180, -90, 90)) extent(x1) <- c(180, 360, -90, 90) rin2 <- merge(x1, x2)

rin2

    rasnum1=1
    
    if(lon1>180)
    {

rasnum1=2; }

      if(lon2>180)
    {

rasnum1=2; }

   if(rasnum1==1)

{ rout1<- crop(rin1, exte1) }

else { rout1<- crop(rin2, exte1) }

rout1

1. plot(rout1)

   heightdelta=sealevel*-1
   

rout2=rout1+heightdelta

rout2[rout2 < 0] <- 0

plot(rout2, col=viridis(100))

1. inorofname2<-"./data_processing/area.nc"
2. inorofname3<-"./data_processing/area_neworog.nc"

crs(rout2) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(rout2, inorofname3, overwrite=TRUE, format="CDF", varname="Elev", varunit="meters",

 longname="Elevation", xname="lon", yname="lat")

crs(rin2) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(rin2, inpath_etopo2, overwrite=TRUE, format="CDF", varname="Elev", varunit="meters",

 longname="Elevation", xname="lon", yname="lat")

writeRaster(rout1, inorofname2, overwrite=TRUE, format="CDF", varname="Elev", varunit="meters",

 longname="Elevation", xname="lon", yname="lat")

1. Pazka

# rout2 #plot(rout2)

}

1. 1. loadvar 3d

load_trace_var10<-function(ivars1, lons1, lats1,variaabeli, unitti, fage,numyears)

{

print("Input variable processing ...")

  	varlocation1<<-((fage-22000)/10*-1)+1
  	
  	
  	print (varlocation1)
  	
  	

if (numyears==1) { vaari1 <-ivars1[,,varlocation1] print (vaari1) plot(vaari1)

1. exit(-1)

}

if (numyears>1) { ## vaari10<-ivars1[,,varlocation1]

vaari10=0 for (nn in 0:(numyears-1)) { vaari10 <-vaari10+ivars1[,,varlocation1+nn] }

vaari1=vaari10/numyears }

   if (variaabeli=="TS")
   {

vaari1<-vaari1-273.15

    }
    

vaari1<- apply(t(vaari1),2,rev)

rrvar1<-raster (vaari1)

rrvar1@extent<-extent(0, 360, -90, 90)

  print("Write ..")

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

writeRaster(rrvar1, invarfname360_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

   rrvar1_180<-rotate(rrvar1)

   rrvar1_180@extent<-extent(-180, 180, -90, 90)

plot(rrvar1_180, col=rainbow(120))

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

writeRaster(rrvar1_180, invarfname180_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

#stop("Joo")

}

tracelocation10 <- function(dirra1, variaabeli, fage, numyears, fmonth) {

   aage2=fage-numyears
   aage1=fage
   

 filename1="trace.01-36.22000BP.clm2.RAIN.22000BP_decavg_400BCE.nc"

1. lista0<-list.files(path=dirra1,pattern="trace*", ignore.case=TRUE)

fagek1=as.integer(aage1) fagek2=as.integer(aage2)

if(fagek2>fagek1) { huba=fagek1 fagek1=fagek2 fagek2=huba

}

filename1<-paste0(predata2, filename1)

print(filename1)

   putin1 <- nc_open(filename1)
    	
    	

print(class(putin1))

    #     stop("k")
   #   ivars0<-0

ivars1 <- ncvar_get(putin1, variaabeli) lones1<- ncvar_get(putin1, lonvar1)

   latis1<- ncvar_get(putin1, latvar1)	

elevlocation1<<-round( (fage-26000)/100)*-1

1. print (lones1)

#  print (latis1)
#  print (ivars1[,,1])
  

1. class(ivars1)
2. print(dim (ivars1))

  #print (ivars1[0,,])
  
  

1. ra0=ivars1[,,4]

1. print (ra0)

1. ra1<-as.raster(ra0)

1. plot (ra1)

  #	varlocation1<<-((fage-age1)*12*-1)+fmonth

 load_trace_var10(ivars1, lons1, lats1,variaabeli, unitti, fage ,numyears)

 }

load_glac_elev<-function() {

number2<-elevlocation1

input2 <- nc_open(inpath_tarasov1)

class(input2)

#input2

print("INPUT GLAC1D ELEV")

elev <- ncvar_get(input2, "HDC") lats2 <- ncvar_get(input2, "YLATGLOBP5") lons2 <- ncvar_get(input2, "XLONGLOB1")

print("INPUT GLAC ELEV DONE")

elev0<-elev[,,number2] elev1<-apply(t(elev0),2,rev) elev1[elev1 <0.0] <- 0.0

rr2<-raster(elev1)

rr2@extent<-extent(0, 360, -90, 90)

2. plot(rr2, col=rainbow(64))

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

writeRaster(rr2, inorofname360_1, overwrite=TRUE, format="CDF", varname="Elev", varunit="Celcius", longname="Elevation", xname="lon", yname="lat")

   rr2_180<-rotate(rr2)
   
   rr2_180@extent<-extent(-180, 180, -90, 90)

plot(rr2_180, col=rainbow(120))

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

writeRaster(rr2_180, inorofname180_1, overwrite=TRUE, format="CDF", varname="Elev", varunit="Celcius", longname="Elevation", xname="lon", yname="lat")

}

1. 1. load_var_nc

load_glac_icemask<-function() {

number2<-elevlocation1

1. inpath22

input2 <- nc_open(inpath_tarasov1)

1. class(input2)

#input2

print("INPUT GLAC1D ICEMASK")

elev <- ncvar_get(input2, "ICEM") lats2 <- ncvar_get(input2, "YLATGLOBP5") lons2 <- ncvar_get(input2, "XLONGLOB1")

print("INPUT GLAC ICEMASK DONE")

elev0<-elev[,,number2] elev1<-apply(t(elev0),2,rev) elev1[elev1 <0.0] <- 0.0

rr2<-raster(elev1)

rr2@extent<-extent(0, 360, -90, 90)

2. plot(rr2, col=rainbow(64))

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

writeRaster(rr2, inicefname360_1, overwrite=TRUE, format="CDF", varname="Icem", varunit="percent", longname="Icemask", xname="lon", yname="lat")

   rr2_180<-rotate(rr2)
   
   rr2_180@extent<-extent(-180, 180, -90, 90)

plot(rr2_180, col=rainbow(120))

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

writeRaster(rr2_180, inicefname180_1, overwrite=TRUE, format="CDF", varname="Icem", varunit="percent", longname="Icemask", xname="lon", yname="lat")

}

1. 1. load_var_nc

load_peltier_elev<-function() { print("Input Peltier 6G 122 elevation data.")

input2 <- nc_open(inpath_peltier1)

ele <- ncvar_get(input2, "stgit") lats2 <- ncvar_get(input2, "Lat") lons2 <- ncvar_get(input2, "Lon") time2 <- ncvar_get(input2, "Time")

   possi2<- -999
   leenu2<-length(time2)
   hage<-(fage/1000)
     

for (n2 in 1:(leenu2) ) { taa2<-time2[n2]

if(hage==taa2) { possi2=n2 break }

if(hage>taa2) { possi2=n2-1 break }

}

elev0<-ele[,,possi2]

elev1<-apply(t(elev0),2,rev) elev1[elev1 <0.0] <- 0.0

rr2<-raster(elev1)

rr2@extent<-extent(0, 360, -90, 90)

plot(rr2, col=rainbow(64))

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

writeRaster(rr2, inorofname360_1, overwrite=TRUE, format="CDF", varname="Elev", varunit="Celcius", longname="Elevation", xname="lon", yname="lat")

   rr2_180<-rotate(rr2)
   
   rr2_180@extent<-extent(-180, 180, -90, 90)

plot(rr2_180, col=rainbow(120))

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

writeRaster(rr2_180, inorofname180_1, overwrite=TRUE, format="CDF", varname="Elev", varunit="Celcius", longname="Elevation", xname="lon", yname="lat")

## peltier ice6g elevation

}

nok_process_rasters<-function() { # WARNING!!!! vprocedure is NOT OK

## NOTE DISTANCE CALCULATION IS OFF SAKE FOR SPEED ## if you use distance, you must it below in script

re_in <- raster("./data_processing/oroin.nc") rt_in<-raster("./data_processing/tempin.nc")

re_in[re_in <0.0] <- 0.0

plot(re_in) plot(rt_in)

#extent(re_in)<-(0,360,-90,90)

big_sabluna<- raster(nrow=180, ncol=360)

#stop("TEST BRK")

extent(big_sabluna)<-c(0,360,-90,90) res(big_sabluna)<-c(1.0, 1.0)

#small_sabluna<-raster(nrow=48, ncol=96) #extent(small_sabluna)<-c(0, 360.000, -90, 90 ) #res(small_sabluna)<-c(3.75, 3,78947) #small_sabluna<-raster(nrow=50, ncol=100) #small_sabluna<-raster(nrow=5, ncol=10) small_sabluna<-raster(nrow=45, ncol=90)

extent(small_sabluna)<-c(0, 360.000, -90, 90 )

1. extent(small_sabluna)<-c(-1.875, 358.125, -89.01354, 89.01354 )
2. res(small_sabluna)<-c(3.75, 3.708898)

re_big <- resample(re_in,big_sabluna, method='bilinear') extent(re_big)<-c(0,360,-90,90) rt_big <- resample(rt_in,big_sabluna, method='bilinear')

1. extent(re_big)<-c(0,360,-90,90)

re_small<-resample(re_in,small_sabluna, method='bilinear')

extent(re_small)<-c(0, 360, -90, 90 )

1. res(re_small)<-c(3.75, 3.75)

rt_small <- resample(rt_in,small_sabluna, method='bilinear')

extent(rt_small)<-c(0,360,-90,90)

1. res(rt_small)<-c(3.75, 3.75)

rlon_big <- rlat_big<- re_big xy <- coordinates(re_big) rlon_big[] <- xy[, 1]

1. JN warning

rlat_big[] <- xy[, 2]

rlat_small<-resample(rlat_big,small_sabluna, method='bilinear') rlon_small<-resample(rlon_big,small_sabluna, method='bilinear')

1. plot(rlat_small)

1. re_big<-rotate(re_big)

1. extent (re_big)

1. plot(re_big)

1. plot (rb_out)
2. extent(rb_out)

1. extent(rt_in)
2. rt_in

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

writeRaster(rt_small, smalltemperaturepath, overwrite=TRUE, format="CDF", varname="tempa", varunit="Units",

 longname="Temp", xname="lon", yname="lat")

1. plot (rt_in)

plot (rt_small)

rt_in

rt_small

1. re_big<-rotate(re_big)
2. re_small<-rotate(re_small)
3. rt_big<-rotate(rt_big)
4. rt_small<-rotate(rt_small)
5. rlon_big<-rotate(rlon_big)
6. rlon_small<-rotate(rlon_small)
7. rlat_big<-rotate(rlat_big)
8. rlat_small<-rotate(rlat_small)

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

writeRaster(re_big, elevpath, overwrite=TRUE, format="CDF", varname="elev", varunit="Units",

 longname="Elevation", xname="lon", yname="lat")

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

writeRaster(re_small,smallelevpath, overwrite=TRUE, format="CDF", varname="elev", varunit="Units",

 longname="Elevation", xname="lon", yname="lat")

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

writeRaster(rlat_big, latpath, overwrite=TRUE, format="CDF", varname="lata", varunit="Units",

 longname="Lata", xname="lon", yname="lat")

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

writeRaster(rlat_small,smalllatpath, overwrite=TRUE, format="CDF", varname="lata", varunit="Units",

 longname="Lata", xname="lon", yname="lat")

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

writeRaster(rlon_big, lonpath, overwrite=TRUE, format="CDF", varname="lata", varunit="Units",

 longname="Longa", xname="lon", yname="lat")

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

writeRaster(rlon_small, smalllonpath, overwrite=TRUE, format="CDF", varname="lata", varunit="Units",

 longname="Longa", xname="lon", yname="lat")

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

writeRaster(rt_big, temperaturepath, overwrite=TRUE, format="CDF", varname="tempa", varunit="Units",

 longname="Tempa", xname="lon", yname="lat")

r_1=re_big r_1[r_1 > 0] <- NA r_1[r_1 < 1] <- 1

if(kalk_distance==1) { rdistance_big <- distance(r_1)

# meters to degrees

rdistance_big<-rdistance_big/111120

rdistance_small<-resample(rdistance_big,small_sabluna, method='bilinear')

rdistance_big<-rotate(rdistance_big) rdistance_small<-rotate(rdistance_small)

#plot(rdistance_big)

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

writeRaster(rdistance_small, smalldistancepath, overwrite=TRUE, format="CDF", varname="distance", varunit="Units", longname="Distance", xname="lon", yname="lat")

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

writeRaster(rdistance_big, distancepath, overwrite=TRUE, format="CDF", varname="distance", varunit="Units", longname="Distance", xname="lon", yname="lat") }

} ## process_rasters

nok_calculate_tables<-function() {

# WARNING procedure is NOT OK

print ("Tables.")

climdir="./data/" smalldir="./data/small/"

elevras_small<-raster("./data/small/small_elevation.nc") distras_small<-raster("./data/small/small_distance.nc") tempras_small<-raster("./data/small/small_temperature.nc")

#elevras_small #distras_small #tempras_small

temp_df<-as.data.frame(tempras_small, xy = TRUE) # jn warn

temp.df[is.na(temp.df)] <- 0.0 elev_df<-as.data.frame(elevras_small, xy = TRUE) dist_df<-as.data.frame(distras_small, xy = TRUE)

#temp_df

plot(tempras_small, col=rainbow(100)) #plot(elevras_small) #plot(distras_small)

small_stack <-stack(elevras_small, tempras_small, distras_small)

#plot(small_stack)

small_df<-data.frame( rasterToPoints( small_stack ) )

#plot(small_df)

## JN warning !!!!

small.df[is.na(small.df)] <- 0.0

#small_df

#colnames(small_df) <- c(library(reshape2) colnames(small_df) <- c("Lon","Lat","Elev", "AvgTemp","distCoast" )

#temp_df

#plot(tempras_small)

colnames(small_df) <- c("Lon","Lat","Elev", "AvgTemp","distCoast" )

small_df["StationID"] <- seq(1,nrow(small_df))+10

#small_df2<- small_df[c("Lat","Lon","Elev", "AvgTemp","distCoast" )]

#StationName,StationID,Lat,Lon,Elev,AvgTemp,distCoast

#colnames(small_df2 )

#head(small_df2,5)

small_df2 <- small_df[c("StationID","Lon","Lat","Elev", "AvgTemp","distCoast" )]

#small_df2

#len(small_df2.AvgTemp) #len(small_df2.Lon)

colnames(small_df2 )

head(small_df2,5)

sample_df=sample_n(small_df2, 600)

## JN warning !!!!

sample_df[is.na(sample_df)] <- 0.0

write.csv(small_df2,"./data/full.csv") write.csv(sample_df,"./data/sample.csv")

##rekt<-dplyr::filter(mag <5.5 & mag > 4.5)

#dplyr::filter(small_df2), lat>

d1<-filter(small_df, Lon >-15.0) d2<-filter(d1, Lon <40.0) d3<-filter(d2, Lat <70.0)

europedf<-filter(d3, Lat >30.0)

write.csv(europedf,"./data/europe.csv")

#plot(europedf)

}

1. 1. Downscaling utilities

downscaler <- function (coarse_rastera, fine_rastera, method) { ## methods: 0 delta, 1 spatialeco, 2 dissever, 3 temperature lapse 6.5 C/1 km lm

   print ("Downscaler()")			

coarse_raster<-coarse_rastera fine_raster<-fine_rastera p1<-fine_raster p2<-fine_raster

plot(fine_raster) plot(coarse_raster, col=viridis(200))

coarseoro<- resample(p1, coarse_raster) coarseoro_big<-resample(coarseoro, p1) orodelta<-(p1-coarseoro_big)

baset1 <- resample(coarse_raster, p1)

raster_stack <- stack(p1,p2)

min_iter <- 5 # Minimum number of iterations max_iter <- 10 # Maximum number of iterations p_train <- 0.1 # Subsampling of the initial data

	 if(method>9999)
	 {

method=2 }

## dissever run

   if(method==2)

{ oma_juttu <- dissever(coarse = coarse_raster, fine = raster_stack, method = "lm", p = p_train, min_iter = min_iter,max_iter = max_iter, verbose=1) orotemp<-oma_juttu$map }

## spatialeco downscale if(method==1) { oma_juttu2 <- raster.downscale(p1, coarse_raster) orotemp<-oma_juttu2$downscale }

1. 1. delta regression 1,1

if(method==0) {

orotemp<-orodelta

   	}

1. 1. delta regression by lapse rate

if(method==3) { orotemp<-orodelta*0.0065*-1

   	}

#biassi=4

#tempiso<-baset1+oma_juttu$map+biassi

coarseorotemp<- resample(orotemp, coarse_raster) coarseorotemp_big<-resample(coarseorotemp, p1)

orotempdelta<-orotemp-coarseorotemp_big

outtemp<-baset1+orotempdelta

1. plot(outtemp, col=rev(rainbow(256)) )

outtempr<-rotate(outtemp)

plot(outtempr)

     return(outtemp)

}

downscale_orography<-function(method) {

  print("Orography ..")

1. 1. downscale orography

   oroko1<<-inorofname180_1;
   
   rasnumi1<<-rasnum1    
   
   "rasnum1"
   print(rasnumi1)
    
   
    if(rasnumi1==2)
    {

print ("ORO 0-360") oroko1<<-inorofname360_1; }

coarse_raster1<-raster(oroko1) fine_raster1<-raster(inorofname3)

  kover_raster1<-resample(coarse_raster1, fine_raster1)
  kover_koeff1<-kover_raster1
  
  kover_koeff1[kover_koeff1>0]<- 1

    plot(fine_raster1)
   plot(coarse_raster1)

   names(fine_raster1)<-"Oroa"

## because europe lies between -15 ... 40, must rotate

1. 1. if(rasnumi1==1)

##    {

1. 1. "rotate oro"
   2. coarse_raster1<-rotate(coarse_raster1)
   3. }

    coarse_raster1[coarse_raster1 < 0] <- 0

plot(fine_raster1, col=viridis(100)) plot(coarse_raster1, col=viridis(100))

outras1<-downscaler(coarse_raster1, fine_raster1,method)

   outras1[outras1<1]<- 0

"Oro ds debug 1"

1. names(outras1)<-"Oro"

1. sabluna1<-raster(inorofname3)

 "sabluna1"
#  plot(sabluna1)
  
#  sabluna1[sabluna1<0]<- 0

1. sabluna1[sabluna1>0]<- 1

1. outras1<-outras1*sabluna1

#  outras1<-outras1*kover_koeff1
  
#  plot(outras1)
  
 

1. 1. elevenaame1<-"./data_processing/area_neworog.nc"
2. elevenaame2<-"./data_processing/area_glacial.nc"

1. plot(outras1, col=rev(plasma(256)))

writeRaster(outras1, outorofname1, overwrite=TRUE, format="CDF", varname="Oro", varunit="meters", longname="Orogr", xname="lon", yname="lat")

}

downscale_temperature<-function(method) {

1. 1. downscale temperature

   print("DS of variable...")
   
    rasnumi1<<-rasnum1   

varoko1<<-invarfname360_1;

if(rasnumi1==1) {

1. coarse_raster2<-rotate(coarse_raster2)

        varoko1<<-invarfname180_1;

}

coarse_raster2<-raster(varoko1)

fine_raster2<-raster(outorofname1)

print("TS rsaters")

#coarse_raster2

#fine_raster2

1. if(klopaali1==1)
2. {
3. fine_raster2=rotate(fine_raster2)
4. }

plot(fine_raster2)

   plot(coarse_raster2)

#fine_raster2<-outras1

1. fine_raster2<-raster(inorofname2)
   1. if(rasnumi1==1)

##    {

1. 1. coarse_raster2<-rotate(coarse_raster2)
   2. }

plot(fine_raster2, col=viridis(100)) plot(coarse_raster2, col=viridis(100) )

outras2<-downscaler(coarse_raster2, fine_raster2,method)

1. names(outras1)<-"TS"

plot(outras2, col=rev(rainbow(256)))

1. writeRaster(outras2, filename='./outdata/result.nc', format="CDF", overwrite=TRUE)

writeRaster(outras2, outvarfname1, overwrite=TRUE, format="CDF", varname="TS", varunit="Celcius", longname="Temperature", xname="lon", yname="lat")

}

normal_ts_downscale<-function()

     {

## normal ## methods 0 delta, 1 spatialeco 2 dissever 3 tempdelta

print("Normal etopo1 ds. Orography ...")

  if(capture_elevation==2)
  {	

# force method to 0 because of peltier topography downscale_orography(0)

  }
  else
  {	

downscale_orography(method_ds_oro)

  }

print("Temperature ...") downscale_temperature(method_ds_var)

1. 1. spatialeco
   2. downscale_orography(1)
   3. downscale_temperature(1)

1. 1. delta

1. 1. downscale_orography(0)
2. downscale_temperature(3)

}

world_ts_downscale<-function()

     {

rasnum1<<- 2 klopaali1<<-1

inorofname1<<-"./data_processing/oroin.nc" invarfname1<<-"./data_processing/tempin.nc"

1. 1. inorofname1<-"./data/rst/elevation.nc"
   2. invarfname1<-"./data/small/small_temperature.nc"

## normal "WORLD DS" # downscale_orography(2) ## WORLD DS #outorofname1<<-"./data/rst/elevation.nc" outorofname1<<-"./data_processing/oroin.nc" # outorofname2<<-"./data_processing/oroin.nc"

downscale_temperature(3) }

custom_ts_downscale_1<-function(orofilee1, tempfilee1)

     {

print("Custom ds type 1.") rasnum1<<- 1

    coarse_raster=raster(tempfilee1)
    fine_raster=raster(orofilee1)
    
   # ex1<-coar
    

    method=3

outras<-downscaler(coarse_raster, fine_raster,method)

plot(outras, col=rev(rainbow(256)))

writeRaster(outras, "./data_result/accurate.nc", overwrite=TRUE, format="CDF", varname="TS", varunit="Celcius", longname="Temperature", xname="lon", yname="lat")

}

plottaus1<-function(zvarnaame1, outfilenaame1) {

1. elevenaame0<-"./data_processing/area_neworog.nc"

elevenaame1<-"../atrace1/data_processing/area_glacial.nc"

   while (dev.cur()>1) dev.off()

grayscale_colors <- gray.colors(100,start = 0.0,end = 1.0,gamma = 2.2,alpha = NULL)

relev1<-raster(elevenaame1)

razteri1<-raster(zvarnaame1)

1. plot(relev1)

aage<<-fage

ram1<-relev1

maintitle1<-paste0("Annual rainfall, ",toString(aage), " BP") subtitle1<-"Trace21ka CCSM3 downscaled" xlab1<-"Lon" ylab1<-"Lat"

print (maintitle1)

limx1<- -15 limx2<-45

   limy1<-35
   limy2<-65	
      
   limitx1<-c(limx1, limx2)
   limity1<-c(limy1, limy2)
   
   templevs1<-c(0,200,400,600,800,1000,1200,1400)
   ltemp1<-length(templevs1)

tempmin1<-templevs1[0]

   tempmax1<-templevs1[ltemp1-1]

print(tempmin1) print(tempmax1)

x1 <- rasterToContour(razteri1, levels=templevs1)

1. 1. lev1 <- seq(tempmin1,tempmax1, by=5)

 lev1=seq(0,1400,200)

class(x1) ek1<-rasterToContour(ram1)

svg(filename="out.svg", width=12, height=10, pointsize=20)

plot( razteri1, main=maintitle1, sub=subtitle1, xlab=xlab1, ylab=ylab1, xlim=limitx1, ylim=limity1, breaks=templevs1, col=rev(viridis(ltemp1)) )

#  plot(x1,xlim=c( -15, 50 ), ylim=c( 30, 60 ), alpha=0.5,add=TRUE)

plot(ram1, xlim=limitx1, ylim=limity1 , zlim=c(0,1), breaks=c(-1,0,1), col=viridis(2) , legend=FALSE, add=TRUE)

1. plot(ek1, add=TRUE, legend=FALSE)

dev.off() }

downscale_raster <- function (coarse_rastera, fine_rastera, method) { ## methods: 0 delta, 1 spatialeco, 2 dissever, 3 temperature lapse 6.5 C/1 km lm

   print ("Downscaler()")			

coarse_raster<-coarse_rastera fine_raster<-fine_rastera p1<-fine_raster p2<-fine_raster

1. plot(fine_raster)
2. plot(coarse_raster, col=viridis(200))

coarseoro<- resample(p1, coarse_raster) coarseoro_big<-resample(coarseoro, p1) orodelta<-(p1-coarseoro_big)

baset1 <- resample(coarse_raster, p1)

raster_stack <- stack(p1,p2)

min_iter <- 5 # Minimum number of iterations max_iter <- 20 # Maximum number of iterations p_train <- 1.0 # Subsampling of the initial data

	 if(method>9999)
	 {

method=2 }

## dissever run

   if(method==2)

{ oma_juttu <- dissever(coarse = coarse_raster, fine = raster_stack, method = "glm", p = p_train, min_iter = min_iter,max_iter = max_iter, verbose=1) orotemp<-oma_juttu$map }

## spatialeco downscale if(method==1) { oma_juttu2 <- raster.downscale(p1, coarse_raster) orotemp<-oma_juttu2$downscale }

1. 1. delta regression 1,1

if(method==0) {

orotemp<-orodelta

   	}

1. 1. delta regression by lapse rate

if(method==3) { orotemp<-orodelta*0.0065*-1

   	}

#biassi=4

#tempiso<-baset1+oma_juttu$map+biassi

coarseorotemp<- resample(orotemp, coarse_raster) coarseorotemp_big<-resample(coarseorotemp, p1)

orotempdelta<-orotemp-coarseorotemp_big

outtemp<-baset1+orotempdelta

1. plot(outtemp, col=rev(rainbow(256)) )

1. outtempr<-rotate(outtemp)

#plot(outtempr)

     return(outtemp)

}

downscale_dissever <- function (coarse_rastera, fine_stack, dismethod, samplerate) {

   print ("Dissever()")		
       names(fine_stack)
       
       
   	

coarse_raster<-coarse_rastera

   p1<-fine_stack$Elevation

1. plot(p1)

1. return(0)

coarseoro<- resample(p1, coarse_raster) coarseoro_big<-resample(coarseoro, p1) orodelta<-(p1-coarseoro_big)

baset1 <- resample(coarse_raster, p1)

raster_stack <- fine_stack

min_iter <- 5 # Minimum number of iterations max_iter <- 10 # Maximum number of iterations p_train <- samplerate # Subsampling of the initial data

oma_juttu <- dissever(coarse = coarse_raster, fine = raster_stack, method = dismethod, p = p_train, min_iter = min_iter,max_iter = max_iter, verbose=1) orotemp<-oma_juttu$map

#tempiso<-baset1+oma_juttu$map+biassi

coarseorotemp<- resample(orotemp, coarse_raster) coarseorotemp_big<-resample(coarseorotemp, p1)

orotempdelta<-orotemp-coarseorotemp_big

outtemp<-baset1+orotempdelta

1. plot(outtemp, col=rev(rainbow(256)) )

1. outtempr<-rotate(outtemp)

#plot(outtempr)

     return(outtemp)

}

geterrormean<-function(rasteri) { error1<-(rasteri-refrain0)

1. plot(error1, col=viridis(6))

errorpros1<-(error1/refrain0)*100 errorpros2<-abs(errorpros1) #errorsd1=cellStats(errorpros1, stat='sd', na.rm=TRUE, asSample=TRUE) #errormean1=cellStats(errorpros1, stat='mean', na.rm=TRUE, asSample=TRUE) errormean2=cellStats(errorpros2, stat='mean', na.rm=TRUE, asSample=TRUE) return(errormean2) }

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")

}

rain_raster_calcutest1 <- function(stak, kerroin) { smallrain0<-raster(smallrainy)

etopo0=stak$Elevation.1

rain<-resample(smallrain0, etopo0) smalletopo<-resample(etopo0,smallrain0) big_smalletopo<-resample(smalletopo, etopo0)

   deltaetopo<-etopo0-big_smalletopo

smallrainzero<-(smallrain0-(smalletopo*kerroin))

big_smallrainzero<-resample(smallrainzero, etopo0)

#routras1=big_smallrainzero+etopo0*kerroin routras1=big_smallrainzero+etopo0*kerroin

#routras<- #routras1<-rain+deltaetopo*1.475 #routras1<-rain+deltaetopo*kerroin #zeroras1<-rain-etopo*kerroin

#return (routras)

#er1=geterrormean(routras1) #print ("----") #print (er1) return(routras1) }

rain_raster_test_loop_1<-function(stak) { sek1<- seq(0.0,5,by=0.25) para0=9999999 er0=9999999 for(para1 in sek1) { #er1=rain_raster_calcutest1(stak, para1) runk=rain_raster_calcutest1(stak, para1) er1=geterrormean(rasteri)

print (para1) print (er1) if(er1<er0) { para0=para1 er0=er1 } }

print( paste("Minimum is:",str(para0)," ", str(er0)) )

}

tpiw <- function(x, w) { m <- matrix(1/(w^2-1), nc=w, nr=w) m[ceiling(0.5 * length(m))] <- 0 f <- focal(x, m) x - f }

preprocess_rasters<-function() { print("PPR") #smallrainy=smallrainame1

if(dsobject==0) { smallrainy=smallrainame1 rainy=refrainame1 smallrain0<-raster(smallrainy)

}

## lgm! if(dsobject==1) { smallrainy<<-smallrainame1 icerain0<<-raster(smallrainy) rainy<<-refrainame2 icemask0<<-raster(icemaskname1) names(icemask0)<<-"Icemask00"

#icerain0[is.na(icerain0)] <- 0 icerain0[icerain0<0] <<- 0

#icerain1<<-icerain0*3600*24*365 icerain1<<-icerain0

icerain1[icerain1>5000] <<- 5000 icerain1[icerain1<0] <<- 1 icerain1<-log(icerain1)

smallrain0<<-icerain0 smallrain1<<-icerain1

}

1. refrain0<<-raster(rainy)

etopo0<<-raster(iname1)

if(calculate_distance==1) { distans0<<-raster(iname2) }

aspect0<<-raster(iname3) tpi0<<-raster(iname4) lonx0<<-raster(iname5) hillshade0<<-raster(iname6) slope0<<-raster(iname7) windir0<<-raster(iname8) bluretopo0<<-raster(iname9) latx0<<-raster(iname10)

etopox0<<-raster(ixname1)

smallrain0_big0<<-resample(smallrain0,etopo0)

if(dsobject==1) { icemask_big0<<-resample(icemask0,etopo0) names(icemask_big0)<-"Icemask0" icemask_big2<<-(icemask_big0-1)*-1 }

1. 1. secunda

  aspect1<<-cos(aspect0)

names(aspect1)<<-"Aspect_Cos"

etopo1<<-etopo0 etopo1[etopo1 < 1] <- 1 etopo1<<-log(etopo1) names(etopo1)<<-"Elevation"

landmask0<<-etopo0 landmask0[landmask0 > 0] <<- 1 landmask0[landmask0 < 1] <<- NA

   plot(landmask0)

if(calculate_distance==1) { distans1<<-distans0 distans1[distans1 < 1] <<- 1 distans1<<-log(1/distans1) names(distans1)<<-"Distance" }

lonx1<<-lonx0 lonx1[lonx1 < 1] <<- 1 lonx1<<-log(1/lonx1) names(lonx1)<<-"Lonx"

bhillshade_slope=6 bhillshade_aspect=260

#aspect1=cos(aspect0-windir0)

bslope <<- terrain(bluretopo0, opt='slope') baspect <<- terrain(bluretopo0, opt='aspect')

bhill <<- hillShade(bslope, baspect, bhillshade_slope, bhillshade_aspect)

## windrose w/zulu basis bhill2 <<- ( hillShade(bslope, baspect, bhillshade_slope, bhillshade_aspect)+ hillShade(bslope, baspect, bhillshade_slope, (bhillshade_aspect-20))+ hillShade(bslope, baspect, bhillshade_slope, (bhillshade_aspect+20))+ hillShade(bslope, baspect, bhillshade_slope*2, (bhillshade_aspect-40))+ hillShade(bslope, baspect, bhillshade_slope*2, (bhillshade_aspect+40))+ hillShade(bslope, baspect, 85, (bhillshade_aspect+180)) )/6

bhillshade1<<-bhill halli1<<-(bhillshade1-minValue(bhillshade1))/(maxValue(bhillshade1)-minValue(bhillshade1)) ## zulu basis rainfall correction

#plot(halli1)

#halli2<-(halli1*0.2)+0.9 #halli2=halli1*0.01

1. halli2<<-2+log((0.1+halli1)*10)

   halli2<<-halli1*0.5+0.75

#plot(halli2) names(halli2)<<-"HillshadeX" halli22<<-halli2 names(halli2)<<-"HillshadeX2"

#halli3=log(1+halli2*100) halli3<<-halli2*-1-2

names(halli3)<<-"HSY"

#baspect2=cos(baspect-4.55)*0.05+1.0 #plot(baspect2)

baspect3a<<-((cos(baspect-4.6)+cos(baspect-5.3)+cos(baspect-3.5))/3) baspect2<<-(baspect3a*0.5)+1.0

   plot(halli2, col=rev(viridis(100)))
   
   plot(baspect2, col=rev(parula(20)))

bhak1<<-baspect2*halli2

   ehak1<<-etopo0*halli2
   ebk1<<-etopo0*baspect2
   ehk1<<-etopo0*bhill
   esk1<<-etopo0*bslope      
   ebhk1<<-etopo0*bhill*baspect2  ##
   ebshk1<<-etopo0*bhill*baspect2*bslope  ## no hyvä

   plot(bhak1, col=rev(viridis(100)))
   plot(ehak1, col=rev(viridis(100)))
    plot(ebk1, col=rev(viridis(100)))
   plot(ehk1, col=rev(viridis(100)))
    plot(esk1, col=rev(viridis(100)))
    plot(ebhk1, col=rev(parula(100))) ## hyvä?
 #        plot(ebshk1, col=rev(viridis(100))) ## no hyvä

names(baspect2)<<-"AspectX" names(ebhk1)<<-"HillShadeAspectX"

#hill <- hillShade(blurslope, bluraspect, hillshade_slope, hillshade_aspect)

tpix1 <<- tpiw(etopo0, w=5)

names(tpix1)<<-"TPIX5"

tpix2 <<- tpiw(etopo0, w=11)

names(tpix2)<<-"TPIX11"

tpix3 <<- tpiw(etopo0, w=15)

names(tpix3)<<-"TPIX15"

tpix4 <<- tpiw(etopo0, w=21)

names(tpix4)<<-"TPIX21"

plot(tpix4)

tpix5 <<- tpiw(etopo0, w=51)

names(tpix5)<<-"TPIX51"

plot(tpix5, col=inferno(50) )

tri1<<-terrain(etopo0, opt="TRI")

names(tri1)<<-"TPI1"

plot(tri1, col=viridis(100))

    roughness1<<-terrain(etopo0, opt="roughness")

names(roughness1)<<-"ROUGHNESS1"

plot(roughness1)

    flowdir1<<-terrain(etopo0, opt="flowdir")

names(flowdir1)<<-"FLOWDIR1"

plot(flowdir1)

}

downscale_rainfall<-function() { print("DS R") ## downscaling suff

#out2<-downscale_raster(smallrain0, etopo0*distans0*lonx0*aspect0,1) #rastafari1<-stack(etopo0,distans0,tpi0,lonx0)

#rastafari1<-stack(etopo1,distans1,tpi0,lonx1) ## kohtalainen, 17.1

## ok jos #rastafari1<-stack(etopo1,bluretopo0, distans1,lonx1,tpi0, baspect2)

#rastafari1<-stack(etopo1, distans1,lonx1,tpi0,halli2)

#rastafari1<-stack(etopo1, distans1,lonx1,latx0,halli2) #rastafari1<-stack(etopo0,halli2)

#names(halli2)<-"Elevation"

#rastafari1<-stack(halli2, baspect2,distans0) ## melko hyvä

#rastafari1<-stack(etopo0,lonx1,latx0,distans1) # 17

#rastafari1<-stack(etopo0,lonx1,latx0,distans1, halli2) # 17 lm

#rastafari1<-stack(etopo0,lonx1,latx0,distans1, halli2, tpix1) # 17 lm log smallrain

#rastafari1<-stack(etopo0,lonx1,latx0,distans1, halli2, tpix1,baspect2) # TOSI HYVÄ KUVIO20 lm log smallrain

#rastafari1<-stack(etopo0,lonx1,latx0,distans1, halli2, tpix1,baspect2) # TOSI HYVÄ KUVIO20 lm log smallrain

## ETOPO= LONX LATX DISTANS HILLSHADEX TPIX ASPECTX

##blurelev1 <- focal(rout3, w=matrix(1,nrow=7,ncol=7), fun=mean, pad=TRUE, na.rm = TRUE)

#rastafari1<-stack(etopo0,etopo0)

#rastafari1<-stack(etopo0,lonx1,latx0,distans1, icemask_big0)

#rastafari1<-stack(etopo0, icemask_big0) ## melko hyvä jääkaudelle #rastafari1<-stack(etopo0, icemask_big0, distans1) ## melko hyvä jääkaudelle #rastafari1<-stack(etopo0, icemask_big0, distans1, halli2) ## melko hyvä jääkaudelle

if(dsobject==1) {

   print(" Ice age ...")

# rastafari1<-stack(etopo0, icemask_big0, distans1, halli2) ## melko hyvä jääkaudelle ## rastafari1<-stack(etopo1, icemask_big0) #rastafari1<-stack(etopo0, icemask_big0, distans1, halli2, baspect2)

##rastafari1<-stack(etopo0, icemask_big0, distans1, halli2, baspect2, tpix1) nok

1. rastafari1<-stack(etopo0, icemask_big0, distans1, halli2) ## kohtalainen?
   1. rastafari1<-stack(g0,g1+g2+g3, g6) ## near ok

    #g1<-etopo0*icemask_big2
       
    g0<-etopo0*1
    
    g1<-etopo0*tpix1
         
    g2<-etopo0*tpix2
    g3<-etopo0*tpix3
    g4<-etopo0*baspect2
    g5<-etopo0*halli2
    g6<-sin(aspect0)
    
    names(g0)<-"Elevation"     
    names(g1)<-"ElevationX1"
    names(g2)<-"ElevationX2"
    names(g3)<-"ElevationX3"
    names(g4)<-"ElevationX4"
    names(g5)<-"ElevationX5"
    names(g6)<-"ElevationX6"
 
                 

1. 1. rastafari1<-stack(g0,g1+g2+g3, g6) ## near ok

1. rastafari1<-stack(etopo0, icemask_big0, distans1,g1+g2+g3,halli2,baspect2) ## near ok

1. rastafari1<-stack(etopo0, icemask_big0, hak1) ## kohtal

rastafari1<-stack(etopo0, icemask_big0, ebhk1) ## kohtal

print(names(rastafari1))

#plot(rastafari1$Elevation)

#rain_raster_test_loop_1(rastafari1) #out3<-rain_raster_calcutest1(rastafari1, 0.0)

#out3a<-downscale_dissever(log(smallrain0), rastafari1,"glm",1.0)

1. out3a<-downscale_dissever(log(icerain1+1), rastafari1,"lm",1.0)

#out3a<-downscale_dissever(icerain1, rastafari1,"rf",0.1)

#     lupu1<-log(icerain0+1)

1. names(lupu1)<-"Rain"
2. pazka

out3a<<-downscale_dissever(smallrain0, rastafari1,"lm",1.0)

1. out3=exp(out3a)

out3<<-out3a

#out4<-((smallrain0_big0*2)+out3)/3

  }

1. out3<-(out3*landmask0)

writeout(out3,outname1,"Rainfall", "mm/yr", "IceAgeRain")

#writeout(smallrain0,origoname1,"Rain Orig.", "mm/yr", "RainO") writeout(icerain1,origoname1,"Ice age rain", "mm/yr", "RainI")

#writeout(out4,outname1)

1. error3<-geterrormean(out3)
2. print ("Errori3")
3. print (error3)

print(minValue(out3))

   print(maxValue(out3))

   plot(out3, col=rev(viridis(50)))
   	

if(stage2==1) { ## try make accurate #out2<-downscale_raster(out3,etopox0,1)

        names(out3)<- "Elevation"
        

1. out2<<-downscale_raster(out3,etopox0,1)

sg1=tpiw(etopo0,3) sg2<-sin(aspect0) sg3<-hillshade0

names(sg1)<- "SG1" names(sg2)<- "SG2" names(sg3)<- "SG3"

# nok rastafari2<-stack(etopo0, sg2) ##

out2<<-downscale_dissever(out3, rastafari2,"lm",1.0)

writeout(out2,outname2,"Rain Stage2", "mm/yr", "RainST2")

#error3<-geterrormean(out2) #print ("Errori3") #print (error2)

       print(minValue(out2))
       print(maxValue(out2))

1. error2<-geterrormean(out2)
2. print ("Errori2")
3. print (error2)

     }

}

1. 1. loadvar 3d

load_trace_var3d_3<-function(ivars1, lons1, lats1,variaabeli, unitti, age1, fage,numyears, fmonth1)

{

1. ivars1 <- ncvar_get(input1, variaabeli)
2. lats1 <- ncvar_get(input1, "lat")
3. lons1 <- ncvar_get(input1, "lon")

print("Input variable processing ...")

  	varlocation1<<-((fage-age1)*12*-1)+fmonth1
  

if (numyears==1) { vaari1 <-ivars1[,,varlocation1] }

if (numyears>1) { ## vaari10<-ivars1[,,varlocation1]

vaari10=0 for (nn in 0:(numyears-1)) { vaari10 <-vaari10+ivars1[,,varlocation1+nn*12] }

vaari1=vaari10/numyears }

   if (variaabeli=="TS")
   {

vaari1<-vaari1-273.15

    }
    

vaari1<- apply(t(vaari1),2,rev)

rrvar1<-raster (vaari1)

rrvar1@extent<-extent(0, 360, -90, 90)

  print("Write ..")

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

writeRaster(rrvar1, invarfname360_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

   rrvar1_180<-rotate(rrvar1)

   rrvar1_180@extent<-extent(-180, 180, -90, 90)

plot(rrvar1_180, col=rainbow(120))

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

writeRaster(rrvar1_180, invarfname180_1, overwrite=TRUE, format="CDF", varname=variaabeli, varunit=unitti, longname=lonkkunimi, xname="lon", yname="lat")

#stop("Joo")

}

wind360 <- function(x) {

   if(x < 0){ 
       x+360
   }else{
       x
   }

}

trace_wind_location_wanha <- function(dirra2, fage, numyears, fmonth) {

   print("Tracewind ...")
   
   aage2=fage-numyears
   aage1=fage

   aagebase0=22000
   lokat1=(aagebase0-aage1)*12
   nummero1=numyears*12
   layeri1=26
   

1. uname1="trace.01.22000-20001BP.cam2.h0.U.0000101-0200012.nc"
2. vname1="trace.01.22000-20001BP.cam2.h0.V.0000101-0200012.nc"

uname1="D:/datav3/trace.01.22000-20001BP.cam2.h0.U.0000101-0200012.nc" vname1="D:/datav3/trace.01.22000-20001BP.cam2.h0.V.0000101-0200012.nc"

 #	uin1 <- nc_open(paste0("D:/datav3/",uname1)
 	

uin1 <- nc_open(uname1)

	vin1 <- nc_open(vname1) 	
 	
 #	print (uin1)
 	
 #	pazka
 	## layer 26, near ground, time 1
 	

uvar0 <- ncvar_get( uin1, varid='U',start=c(1,1,26,lokat1+fmonth), count=c(-1,-1,1,1)) vvar0 <- ncvar_get( vin1, varid='V',start=c(1,1,26,lokat1+fmonth), count=c(-1,-1,1,1))

uvars0 <- ncvar_get( uin1, varid='U',start=c(1,1,layeri1,lokat1), count=c(-1,-1,1,nummero1)) vvars0 <- ncvar_get( vin1, varid='V',start=c(1,1,layeri1,lokat1), count=c(-1,-1,1,nummero1))

#  print (dim (uvars0))

   vfis0<-uvars0
   vrrs0<-uvars0

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

for (n in 1:12 ) {

  			uu=uvars0[,,n]

vv=vvars0[,,n] # print (dim(uu)) vf00=90.0+180.0*atan2(vv,uu)/3.14159

vv0 <- sqrt(uu*uu+vv*vv) vfis0[,,n]=vf00 vrrs0[,,n]=vv0

vf02=t(vf00)

vf03<-apply(vf02,2,rev)

vfr0=raster(vf03) vfr1<- overlay(vfr0,fun = wind360)

extent(vfr1)<-vext1 names(vfr1)<-"windir"

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

stackfis1 <- stack( stackfis1 , vfr2 )

}

1. uvars2=rowMeans(uvars0, dims = 2)
2. vvars2=rowMeans(vvars0, dims = 2)

1. uvars2=rowMeans(uvars0, dims = 2)
2. vvars2=rowMeans(vvars0, dims = 2)

uvars2=uvar0 vvars2=vvar0

#print(dim (uvars0))

1. print(uvars2)

1. pazka

uvar2=t(uvars2) vvar2=t(vvars2)

1. uvar2=t(uvar0)
2. vvar2=t(vvar0)

uvar1<-apply(uvar2,2,rev) vvar1<-apply(vvar2,2,rev)

u0=raster(uvar1) v0=raster(vvar1)

extent(u0)<-vext1 extent(v0)<-vext1 names(u0)<-"u" names(u0)<-"v"

   sxlats1<-seq(0,360,96/360)
   sylats1<-seq(-90,90,48/180)   
   	
   vfi000=90.0+180.0*atan2(v0,u0)/3.14159

   vfi0 <- overlay(vfi000,fun = wind360)

   vvel0<-sqrt(u0*u0+v0*v0)

extent(vfi0)<-vext1 extent(vvel0)<-vext1

vfi1<-rotate(vfi0) vvel1<-rotate(vvel0)

extent(u0)<-vext1 extent(v0)<-vext1

   u1<-rotate(u0)
   v1<-rotate(v0)
   crs(u1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
   crs(v1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

   windstack1<-stack(u1,v1)

df <- expand.grid(x=seq(-2, 2, .1), y=seq(-2, 2, .1)); k=2; b=4 df$z <- with(df, (y^2)-(x^2) ) library(raster) library(rasterVis) r <- rasterFromXYZ(df)

svg(filename="Std_SVG.svg", width=5, height=4, pointsize=12)

projection(r) <- CRS("+proj=longlat +datum=WGS84") vectorplot(r, par.settings=RdBuTheme)

dev.off()

# plot(vfi0)

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

writeRaster(vfi0, "windir_360.nc", overwrite=TRUE, format="CDF", varname="Windir", varunit="u", longname="Windir", xname="lon", yname="lat")

1. plot(vfi1)

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

1. writeRaster(vvel1, "winvel_180.nc", overwrite=TRUE, format="CDF", varname="Winvel", varunit="u",
2. longname="Winvel", xname="lon", yname="lat")

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

1. writeRaster(vfi1, "windir_180.nc", overwrite=TRUE, format="CDF", varname="Windir", varunit="u",
2. longname="Windir", xname="lon", yname="lat")

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

writeRaster(stackfis1, windpath2, overwrite=TRUE, format="CDF", xname="lon", yname="lat")

   print("Stak ok")

1. 1. test

names(dim(vfis0)) <- c('lon', 'lat', 'time') 

ArrayToNc(list(vfis0 = vfis0),windpath1)

}

load_monthly_wind_and_rain <- function(dirra2, fage, numyears, fmonth) {

   print("Tracewind ...")
   
   aage2=fage-numyears
   aage1=fage

   aagebase0=22000
   lokat1=(aagebase0-aage1)*12
   nummero1=numyears*12
   layeri1=26
   

uname1="D:/datav3/trace.01.22000-20001BP.cam2.h0.U.0000101-0200012.nc" vname1="D:/datav3/trace.01.22000-20001BP.cam2.h0.V.0000101-0200012.nc" rcname1="D:/datav3/trace.01.22000-20001BP.cam2.h0.PRECC.0000101-0200012.nc" rlname1="D:/datav3/trace.01.22000-20001BP.cam2.h0.PRECL.0000101-0200012.nc"

uin1 <- nc_open(uname1)

	vin1 <- nc_open(vname1) 	
 	rcin1 <- nc_open(rcname1)  	
	rlin1 <- nc_open(rlname1) 
	

uvars0 <- ncvar_get( uin1, varid='U',start=c(1,1,layeri1,lokat1), count=c(-1,-1,1,nummero1)) vvars0 <- ncvar_get( vin1, varid='V',start=c(1,1,layeri1,lokat1), count=c(-1,-1,1,nummero1))

rcvars0 <- ncvar_get( rcin1, varid='PRECC',start=c(1,1,lokat1), count=c(-1,-1,nummero1)) rlvars0 <- ncvar_get( rlin1, varid='PRECL',start=c(1,1,lokat1), count=c(-1,-1,nummero1))

  # print (dim(rcvars0))

 #pazka
 
 

#  print (dim (uvars0))

   vfis0<-uvars0
   vrrs0<-uvars0

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

stackrain1<-stack()

for (n in 1:12 ) {

  			uu=uvars0[,,n]

vv=vvars0[,,n]

  			rain00=rcvars0[,,n]+rcvars0[,,n]

# print (dim(uu)) vf00=90.0+180.0*atan2(vv,uu)/3.14159

vv0 <- sqrt(uu*uu+vv*vv) vfis0[,,n]=vf00 vrrs0[,,n]=vv0

vf02=t(vf00) rain01=t(rain00) rain02<-apply(rain01,2,rev)

vf03<-apply(vf02,2,rev)

vfr0=raster(vf03)

rain0=raster(rain02)

vfr1<- overlay(vfr0,fun = wind360)

extent(vfr1)<-vext1 names(vfr1)<-"windir"

extent(rain0)<-vext1 names(rain0)<-"RAIN"

vfr2<-rotate(vfr1) rain2=rotate(rain0)

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

rain2<-rain2*1000*3600*24*30.5

stackfis1 <- stack( stackfis1 , vfr2 ) stackrain1 <- stack( stackrain1 , rain2 )

}

crs(stackfis1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(stackfis1, windpath2, overwrite=TRUE, format="CDF", xname="lon", yname="lat", varname="windir", varunit="deg",

 longname="windir")

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

writeRaster(stackrain1, rainpath2, overwrite=TRUE, format="CDF", xname="lon", yname="lat", varname="Rain", varunit="mm",

 longname="Rain")     
   
   
   print("Stak ok")

1. 1. test

# names(dim(vfis0)) <- c('lon', 'lat', 'time') 

#ArrayToNc(list(vfis0 = vfis0),windpath1)

1. 1. loadwindrain

 }

downscale_one_month_wind_rainfall<-function(wind0, rain0) { print("DS MONTH Rainfall.") graphics.off() if(dsobject==1) {

      windfi1<-resample(wind0, etopo0)
     
     windfi2<-deg2rad(windfi1)
    
    
    plot(windfi1)
#      plot(windfi2)
  #    plot(aspect0)
 #     plot(baspect)
     
     saspect0<<-windfi2-aspect0+3.14159
     sbaspect0<-windfi2-baspect+3.14159
     
     saspect1<-cos(saspect0)
     sbaspect1<-cos(sbaspect0)    
     
     sbaspectopo1<-sbaspect1*etopo0
     
      
    #  plot(saspect0)
    #  plot(sbaspect0)  
    #  plot(saspect1)
    #  plot(sbaspect1)  

     #    plot(sbaspectopo1)  

  #    	bhillshade_slope=6

1. bhillshade_aspect=260

#aspect1=cos(aspect0-windir0)

1. bslope <<- terrain(bluretopo0, opt='slope')
2. baspect <<- terrain(bluretopo0, opt='aspect')

1. bhill <<- hillShade(bslope, baspect, bhillshade_slope, bhillshade_aspect)

   #   	rastafari1<<-stack(etopo0, icemask_big0, ebhk1) ## kohtal

      	rastafari1<<-stack(etopo0, icemask_big0, ebhk1) ## kohtal
   

out3<<-downscale_dissever(rain0, rastafari1,"lm",1.0)

# plot(out3)

# pazka

return(out3)

  }

}

month_rain_wind_loop<-function() { inrain1<-stack(rainpath2) inwindir1<-stack(windpath2)

print(names(inrain1)) print(names(inwindir1))

g0<<-etopo0*1

g1<<-etopo0*tpix1

g2<<-etopo0*tpix2 g3<<-etopo0*tpix3 g4<<-etopo0*baspect2 g5<<-etopo0*halli2 g6<<-sin(aspect0)

names(g0)<<-"Elevation" names(g1)<<-"ElevationX1" names(g2)<<-"ElevationX2" names(g3)<<-"ElevationX3" names(g4)<<-"ElevationX4" names(g5)<<-"ElevationX5" names(g6)<<-"ElevationX6"

# rastafari1<<-stack(etopo0, icemask_big0, ebhk1) ## kohtal

# print(names(rastafari1))

outrainstack1<-stack() outrain1<<-etopo0*0

for (n in 1:12 ) { #raname=paste0("X"+as.string(n) ) #print (raname) rain0=inrain1n wind0=inwindir1n

1. print (dim(rain0))

#print(maxValue(rain0)) raso1<-downscale_one_month_wind_rainfall(wind0, rain0)

outrainstack1<-stack(outrainstack1,raso1)

outrain1<-outrain1+raso1

}

#  outrain1<-outrain1/12
 
   

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

writeRaster(outrainstack1, outrainpath2, overwrite=TRUE, format="CDF", xname="lon", yname="lat",varname="Rain", varunit="mm",

   longname="Rain")     
   

crs(outrain1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(outrain1, outrainpath3, overwrite=TRUE, format="CDF", xname="lon", yname="lat", varname="Rain", varunit="mm",

 longname="Rain")       	

}

rainfarmer<-function() {

hipi=16 # inp raster dim narow=16 ## scaling factor

dipi=narow*hipi

#rainer0<<-raster(smallrainame1) rainer00<<-raster(outname1)

1. etopo000<<-raster(indirname3)

ext1<-c(superlon1,superlon2,superlat1, superlat2) ext2<-c(lon1,lon2,lat1, lat2)

rainer1<-crop(rainer00, ext1)

   names(rainer1)<-"Rain"

etopor0<-crop(etopo_base0, ext1)

  names(etopor0)<-"Elev"
  # print(".")

xsiz1=ysiz1=hipi

xy <- matrix(rnorm(xsiz1*ysiz1),xsiz1,ysiz1)

rast0 <- raster(xy) extent(rast0) <- c(superlon1,superlon2,superlat1, superlat2) crs(rast0) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

etopor1=resample(etopor0,rast0) etopor1[is.na(etopor1[])] <- 0

names(etopor1)<-"Elev" wr1=(etopor1-minValue(etopor1))/(maxValue(etopor1)-minValue(etopor1)) names(wr1)<-"Elev"

#plot(rainer1)

 #  print(".")
   
   

r1<-resample(rainer1,rast0)

 names(r1)<-"Rain"

  plot(r1)
  plot(wr1)
  

r1[is.na(r1[])] <- 0 wr1[is.na(wr1[])] <- 0 #names(r1)

r2<-r1$Rain wr2<-wr1$Elev

print(dim (r2)) print(dim (wr2)) #pazka

1. r=r2[,,1]*3600*24*365

r=r2

#   print(".")
   

r <- matrix(r, nrow = hipi, byrow = TRUE) wr3=wr2[,,1]

#print (r)

#pazka

lon <- seq(superlon1,superlon2 , 360/96) lat <- seq(superlat1,superlat2, 180/48)

   print (dim(r))

dim(r) <- c(hipi, hipi, 1) dim(wr3) <- c(hipi, hipi, 1)

nf <- narow

1. print("K")

rd <- rainfarm(r, 1.7, nf, fsmooth = TRUE)

1. plot(rd)

grid <- lon_lat_fine(lon, lat, nf)

#grid$lon[1:4]

dos1<-rd[,,1]

rout1<-raster(dos1)

plot(rout1)

extent(rout1)<-ext1

rout2<-crop(rout1,ext2)

     names(rout2)<-"rainfarmer"

   print(" .. ")	
   return(rout2)

}

rainfarm_runner<- function() {

print ("Rainfarm ...")

1. etopo000<-raster()

etopo_base0<<-raster(inpath_etopo1)

raka0=rainfarmer()

 print("*")	

for( n in 1:20) { raka1=rainfarmer() raka0=raka0+raka1 }

rafa1=raka0/10

crs(rafa1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(rafa1, filename=rainfarm_name1, overwrite=TRUE, format="CDF", varname="Rainfarm", varunit="mm", longname="Rainfarm", xname="lon", yname="lat")

# print("*")	

}

1. 1. main running ...

print("Teace21ka RAIN downscaler.") print ("Main program running ...")

read_argus()

print("Age") print(fage) print("Numyears") print(numyears) print("Months") print(fmonth)

sealevel_from_age(fage)

if (orography_preprocess==1) { print("Orography ...") preprocess_orography() }

if(capture_variable==1) { print("Processing input variable to snapshot ...") tracelocation10(predata2, variaabeli, fage, numyears, fmonth); }

if(capture_wind_rain_variables==1) { print("Capture monthly rain and wind variable ...") load_monthly_wind_and_rain (datastore1, fage, numyears, fmonth) }

if(capture_elevation==1) { print("Glac elevation ...") load_glac_elev() load_glac_icemask() }

if(capture_elevation==2) {

print("Peltier elevation")
 load_peltier_elev()

}

1. 1. this section nok

1. if(kalk_tables==1)
2. {
   1. process_rasters()
3. calculate_tables()
4. }

1. 1. .. section nok

1. print ("Downscaling ...")

if(capture_rain_orography==1) { acquire_rain_orography_data() }

if(rain_monthly_downscaling==1) { preprocess_rasters() month_rain_wind_loop() }

1. 1. do downscaling

if(normal_ds==1) { print("Normal etopo1 area downscaling.") ## default area DS

1. normal_ts_downscale()

preprocess_rasters() downscale_rainfall()

}

if (run_rainfarm==1) { print("Run rainfarm.") rainfarm_runner() }

if(global_ds==1) { print("World downscaling.")

1. world_ts_downscale()

}

if(accurate_ds==1) { print ("Customized dwonscaling.")

   inris1<<-out3	

rs1<-raster(kustomorofilee1)

out4<-downscale_raster(inris1,rs1,1) writeout(out4,"./indata_out/accurate.nc","Rainfall", "mm/yr", "IceAgeRain")

1. custom_ts_downscale_1(kustomorofilee1,outvarfname1 )

}

if(plot_result==1) { print("Plotting of draft.")

 plottaus1("./indata_out/raintest.nc",plotfname1 )

}

1. try(system("python maplot1.py", intern = TRUE, ignore.stderr = TRUE))

print("Program run done.")