File:Europe rainfall trace21ka 21000BP annual.svg
Original file (SVG file, nominally 1,650 × 1,275 pixels, file size: 1.46 MB)
Captions
Summary
[edit]DescriptionEurope rainfall trace21ka 21000BP annual.svg |
English: Annual rainfall in Europe, 21000 yearsa ago. |
|||
Date | ||||
Source | Own work | |||
Author | Merikanto | |||
SVG development InfoField |
|
This image is based on trace21ka post-processed data, Tarasov Glac1d dataset and etopo1,
This image is based on Trace21ka dataset on Earth System Grid.
TraCE-21ka: Simulation of Transient Climate Evolution over the last 21,000 years
https://www.cgd.ucar.edu/ccr/TraCE/ https://www.earthsystemgrid.org/project/trace.html https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm.trace.html
Otto-Bliesner, B. & Liu, Zhensheng & He, Feng & Brady, Esther & Clark, P. & Carlson, Andrew & Tomas, Robert & Erickson, D. & Jacob, Robert. (2008). Transient Simulation of Climate Evolution over the Last 21,000 years (TraCE21,000): First Results. AGU Fall Meeting.
Downscaling : zulu basis rain shadows, assumes mostly wesrerly winds,
Data is processed with R script(s) and visualized with NASA Panoply.
-
- Trace21ka Ccsm3 Glac1D Etopo1 Rainfall "variable" downscaler
- zulu basis rain shadows for Europe
-
- Language "R" 3.6.1 2019
- tested under Windows 10 Rscript
- 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)
- library(OceanView)
library(sp)
library(reshape2)
library(dplyr)
library(REdaS)
library(easyNCDF)
library(numbers)
- library(gdata)
library(abind)
library(Cairo)
- NOTE SET THESE
- set orografy creating off
- argreading=0 ## WARNING NOT IMPLEMENTED default 0
- kalk_distance=0 # for speed 0, if u not neet to calc distance
- kalk_tables=0 # suppress erroris from tables, id kalk_distance=0
- must be 1, if you create output foor certain year, first
- Control variables
arg_reading =0
orography_preprocess=0
- output area type
- 0: normal local output
- 1: global output, according to Glac!D
- NOT OK 2: custom orography, local NOK
- acquire trace21ka variable data, default 1
capture_variable=0
- capture **MONTHLY** wind and rain
capture_wind_rain_variables=1
- acquire glac1d elevation data,
- default 1 ( Tarasov glac1d)
- 2 Peltier 1ce 6G-D 122k
capture_elevation=0
- 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
- monthly downscaling
rain_monthly_downscaling=1
- normal one year etopo1 single area downscaling
normal_ds=0
- rainfarm test utility
run_rainfarm=9
stage2=0 ## 0: no "stage 2", 1: "stage 2"
- either global_ds or normal_ds must set 1
- global area downscaling to glac1d
global_ds=0
- ! afterburner accurate "kustomorofilee1.nc" downscaling, default 0
- need accurate orog file (GTOPO30 slice, SRTM etc.)
accurate_ds=0
- kustomorofilee1="./predata/europe30.nc"
- plotting is under alpha stage!
plot_result=1
- downscaling method from orography
- 0 delta, 1 spatialeco, 2 dissever, 3 temperature delta lapse rate 6,5 c/ 1 km
- method 2 slow
method_ds_oro=1
- downscaling method for temperature
method_ds_var=3
kustomorofilee1="../atrace1/predata/dordogne.nc"
- year to render
- fage=14559
fage=21000
- number of years to average, since "fage"
numyears=1
- month of year to render
fmonth=1
- sea level, below current: eq 120 means height of -120 to current
- not needed set by default, calculated from age w/ thisscript
sealevel=-110
- area to render, lon1, lon2 lat1 lat2 rectangle
- europe
lon1=-15.0
lon2=40.0
lat1=30.0
lat2=70.0
- NOTE : use longitude 200 instead of -160
- beringia
- lon1=160
- lon2=240
- lat1=50.0
- lat2=80.0
- global
- klopaali1=1
- lon1=-180.0
- lon2=180.0
- lat1=-90.0
- lat2=90.0
- dordogne 1
- lon1=-2.0
- lon2=4.0
- lat1=42.0
- lat2=48.0
- dordogne 2
- lon1=0.0
- lon2=2.0
- lat1=44.0
- lat2=45.5
- dordogne 3
- lon1=-1.0
- lon2=2.0
- lat1=44.0
- lat2=45.5
filename1="KKK"
tslocation1=-999
elevlocation1=-999
variaabeli="RAIN"
unitti="mm"
lonkkunimi="Rain"
lonvar1="lon"
latvar1="lat"
- input data dirs
- etopo1 dir
predata="./predata/"
- Trace21ka TS dir
predata2="./predata2/"
datastore1="d:/datav3/"
dir.create("./data_processing")
dir.create("./plotz")
dir.create("./indata_wind/")
dir.create(predata)
- outpath11<-"./data_processing/area.nc"
- outpath12<-"./data_processing/area_neworog.nc"
- glac1D HDC file
inpath_tarasov1<-"../atrace1/predata/TOPicemsk.GLACD26kN9894GE90227A6005GGrBgic.nc"
inpath_peltier1<-"../atrace1/predata/IceT.I6F_C.131QB_VM5a_1deg.nc"
- outpath22<-"./data_processing/oroin.nc"
- ! grid-registreted etopo1
- inpath11<-"./predata/ETOPO1_Ice_g_gdal.grd"
- cell reg maybe inaccurate, but func w/dissever
- inpath_etopo1<-"../trace1/predata/ETOPO1_Ice_c_gdal.grd"
inpath_etopo1<-"D:/datavarasto/etopo1.nc"
- inpath_etopo1<-"./predata/ETOPO1_Ice_g_geotiff.tif"
inpath_etopo2<- "./predata/etopo1_360.nc"
data_processing="./data_processing/"
windpath1="./indata_wind/windirs1.nc"
- wind, rain vars 12 montsh
windpath2="./indata_wind/windirs2.nc"
rainpath2="./indata_rain/rain12.nc"
- 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)
- 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"
}
- unlink("indata_rain")
- unlink("indata_map")
- unlink("indata_out")
- 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"
- france
- reso1min=2.5
- coarsereso1deg=1
- france
- lon1=-5
- lat1=42
- lon2=8
- lat2=50
- 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
- rainfarm doms lat, lon
superlon1=-20
superlon2=70
superlat1=30
superlat2=90
- superlon1=-20
- superlon2=40
- superlat1=30
- superlat2=90
- 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"
- 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)
- 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")
}
- CODE BEGINS
exte1<- extent(lon1,lon2,lat1,lat2)
- global doordinates: pacific or europe centered
- rasnum 1: -180 - 180
- 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 )
- 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
- 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
- plot(rout1)
heightdelta=sealevel*-1
rout2=rout1+heightdelta
rout2[rout2 < 0] <- 0
plot(rout2, col=viridis(100))
- inorofname2<-"./data_processing/area.nc"
- 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")
- Pazka
# rout2
#plot(rout2)
}
- 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)
- 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"
- 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
- print (lones1)
# print (latis1)
# print (ivars1[,,1])
- class(ivars1)
- print(dim (ivars1))
#print (ivars1[0,,])
- ra0=ivars1[,,4]
- print (ra0)
- ra1<-as.raster(ra0)
- 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)
- 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")
}
- load_var_nc
load_glac_icemask<-function()
{
number2<-elevlocation1
- inpath22
input2 <- nc_open(inpath_tarasov1)
- 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)
- 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")
}
- 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 )
- extent(small_sabluna)<-c(-1.875, 358.125, -89.01354, 89.01354 )
- 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')
- 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 )
- 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)
- res(rt_small)<-c(3.75, 3.75)
rlon_big <- rlat_big<- re_big
xy <- coordinates(re_big)
rlon_big[] <- xy[, 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')
- plot(rlat_small)
- re_big<-rotate(re_big)
- extent (re_big)
- plot(re_big)
- plot (rb_out)
- extent(rb_out)
- extent(rt_in)
- 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")
- plot (rt_in)
plot (rt_small)
rt_in
rt_small
- re_big<-rotate(re_big)
- re_small<-rotate(re_small)
- rt_big<-rotate(rt_big)
- rt_small<-rotate(rt_small)
- rlon_big<-rotate(rlon_big)
- rlon_small<-rotate(rlon_small)
- rlat_big<-rotate(rlat_big)
- 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)
}
- 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
}
- delta regression 1,1
if(method==0)
{
orotemp<-orodelta
}
- 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
- plot(outtemp, col=rev(rainbow(256)) )
outtempr<-rotate(outtemp)
plot(outtempr)
return(outtemp)
}
downscale_orography<-function(method)
{
print("Orography ..")
- 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
- if(rasnumi1==1)
## {
- "rotate oro"
- coarse_raster1<-rotate(coarse_raster1)
- }
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"
- names(outras1)<-"Oro"
- sabluna1<-raster(inorofname3)
"sabluna1"
# plot(sabluna1)
# sabluna1[sabluna1<0]<- 0
- sabluna1[sabluna1>0]<- 1
- outras1<-outras1*sabluna1
# outras1<-outras1*kover_koeff1
# plot(outras1)
- elevenaame1<-"./data_processing/area_neworog.nc"
- elevenaame2<-"./data_processing/area_glacial.nc"
- 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)
{
- downscale temperature
print("DS of variable...")
rasnumi1<<-rasnum1
varoko1<<-invarfname360_1;
if(rasnumi1==1)
{
- coarse_raster2<-rotate(coarse_raster2)
varoko1<<-invarfname180_1;
}
coarse_raster2<-raster(varoko1)
fine_raster2<-raster(outorofname1)
print("TS rsaters")
#coarse_raster2
#fine_raster2
- if(klopaali1==1)
- {
- fine_raster2=rotate(fine_raster2)
- }
plot(fine_raster2)
plot(coarse_raster2)
#fine_raster2<-outras1
- fine_raster2<-raster(inorofname2)
- if(rasnumi1==1)
## {
- coarse_raster2<-rotate(coarse_raster2)
- }
plot(fine_raster2, col=viridis(100))
plot(coarse_raster2, col=viridis(100) )
outras2<-downscaler(coarse_raster2, fine_raster2,method)
- names(outras1)<-"TS"
plot(outras2, col=rev(rainbow(256)))
- 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)
- spatialeco
- downscale_orography(1)
- downscale_temperature(1)
- delta
- downscale_orography(0)
- downscale_temperature(3)
}
world_ts_downscale<-function()
{
rasnum1<<- 2
klopaali1<<-1
inorofname1<<-"./data_processing/oroin.nc"
invarfname1<<-"./data_processing/tempin.nc"
- inorofname1<-"./data/rst/elevation.nc"
- 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)
{
- 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)
- 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)
- 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)
- 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
- 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 <- 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
}
- delta regression 1,1
if(method==0)
{
orotemp<-orodelta
}
- 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
- plot(outtemp, col=rev(rainbow(256)) )
- 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
- plot(p1)
- 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
- plot(outtemp, col=rev(rainbow(256)) )
- outtempr<-rotate(outtemp)
#plot(outtempr)
return(outtemp)
}
geterrormean<-function(rasteri)
{
error1<-(rasteri-refrain0)
- 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
}
- 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
}
- 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
- 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
- rastafari1<-stack(etopo0, icemask_big0, distans1, halli2) ## kohtalainen?
- 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"
- rastafari1<-stack(g0,g1+g2+g3, g6) ## near ok
- rastafari1<-stack(etopo0, icemask_big0, distans1,g1+g2+g3,halli2,baspect2) ## near ok
- 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)
- out3a<-downscale_dissever(log(icerain1+1), rastafari1,"lm",1.0)
#out3a<-downscale_dissever(icerain1, rastafari1,"rf",0.1)
# lupu1<-log(icerain0+1)
- names(lupu1)<-"Rain"
- pazka
out3a<<-downscale_dissever(smallrain0, rastafari1,"lm",1.0)
- out3=exp(out3a)
out3<<-out3a
#out4<-((smallrain0_big0*2)+out3)/3
}
- 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)
- error3<-geterrormean(out3)
- print ("Errori3")
- 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"
- 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))
- error2<-geterrormean(out2)
- print ("Errori2")
- print (error2)
}
}
- loadvar 3d
load_trace_var3d_3<-function(ivars1, lons1, lats1,variaabeli, unitti, age1, fage,numyears, fmonth1)
{
- ivars1 <- ncvar_get(input1, variaabeli)
- lats1 <- ncvar_get(input1, "lat")
- 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
- uname1="trace.01.22000-20001BP.cam2.h0.U.0000101-0200012.nc"
- 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 )
}
- uvars2=rowMeans(uvars0, dims = 2)
- vvars2=rowMeans(vvars0, dims = 2)
- uvars2=rowMeans(uvars0, dims = 2)
- vvars2=rowMeans(vvars0, dims = 2)
uvars2=uvar0
vvars2=vvar0
#print(dim (uvars0))
- print(uvars2)
- pazka
uvar2=t(uvars2)
vvar2=t(vvars2)
- uvar2=t(uvar0)
- 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")
- plot(vfi1)
crs(vvel1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
- writeRaster(vvel1, "winvel_180.nc", overwrite=TRUE, format="CDF", varname="Winvel", varunit="u",
- longname="Winvel", xname="lon", yname="lat")
crs(vfi1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
- writeRaster(vfi1, "windir_180.nc", overwrite=TRUE, format="CDF", varname="Windir", varunit="u",
- 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")
- 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")
- test
# names(dim(vfis0)) <- c('lon', 'lat', 'time')
#ArrayToNc(list(vfis0 = vfis0),windpath1)
- 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
- 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)
# 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
- 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)
- 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
- 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
- print("K")
rd <- rainfarm(r, 1.7, nf, fsmooth = TRUE)
- 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 ...")
- 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("*")
}
- 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()
}
- this section nok
- if(kalk_tables==1)
- {
- process_rasters()
- calculate_tables()
- }
- .. section nok
- print ("Downscaling ...")
if(capture_rain_orography==1)
{
acquire_rain_orography_data()
}
if(rain_monthly_downscaling==1)
{
preprocess_rasters()
month_rain_wind_loop()
}
- do downscaling
if(normal_ds==1)
{
print("Normal etopo1 area downscaling.")
## default area DS
- normal_ts_downscale()
preprocess_rasters()
downscale_rainfall()
}
if (run_rainfarm==1)
{
print("Run rainfarm.")
rainfarm_runner()
}
if(global_ds==1)
{
print("World downscaling.")
- 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")
- custom_ts_downscale_1(kustomorofilee1,outvarfname1 )
}
if(plot_result==1)
{
print("Plotting of draft.")
plottaus1("./indata_out/raintest.nc",plotfname1 )
}
- try(system("python maplot1.py", intern = TRUE, ignore.stderr = TRUE))
print("Program run done.")
Licensing
[edit]- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
- share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 14:25, 20 October 2019 | 1,650 × 1,275 (1.46 MB) | Merikanto (talk | contribs) | Sum from monthly dataset | |
10:56, 18 October 2019 | 1,650 × 1,275 (1.4 MB) | Merikanto (talk | contribs) | User created page with UploadWizard |
You cannot overwrite this file.
File usage on Commons
There are no pages that use this file.
Metadata
This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Width | 1650px |
---|---|
Height | 1275px |