File:40750bp europe precipitation annual 1.png
Original file (1,472 × 976 pixels, file size: 1.01 MB, MIME type: image/png)
Captions
Summary
[edit]Description40750bp europe precipitation annual 1.png |
English: Annual precipitation in Europe, 40750 years ago. |
Date | |
Source | Own work |
Author | Merikanto |
The source of data to produce this image is
Armstrong et al. 2019: A simulated Northern Hemisphere land based climate dataset for the past 60,000 years.
Article https://www.nature.com/articles/s41597-019-0277-1#citeas
Data on CEDA archive
http://data.ceda.ac.uk/badc/deposited2018/HadCM3B_60Kyr_Climate/data/precip
"R" code to downscale this image.
Uses topography files Etopo1, Tarasov GLAC1D dataset
Panoply visualization.
Code:
- Hadcm3b 60kyr netcdf downscaler
- Peltier Glac1D Etopo1
-
- Language "R" 4.0.2
- v 1015.05
- 28.9.2020
-
- Externel requirements:
- You need in most cases Etopo1, Peltier and Tarasov Glac 1d dataset to
- run this script
-
- WARNING1: Maybe inaccuracy due to inaccurate input orography data
- WARNING2: plotting is under development stage, must visualize w/Panoply or other tools
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("pals")
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("Rfast")
}
library(raster)
library(rgdal)
library(ncdf4)
library(lattice)
library(maptools)
library(rgeos)
library(spatialEco)
library(dissever)
library(RColorBrewer)
library(viridis)
library(pals)
library(data.table)
library(stringr)
library(rlist)
library(pipeR)
library(Rfast)
library(sp)
library(reshape2)
library(dplyr)
- 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=1
orography_preprocess=1
- get dordogne srtm
get_srtm_data=0
- gtopo30
get_gtopo30_data=0
- output area type
- 0: normal local output
- 1: global output, according to Glac!D
- NOT OK 2: custom orography, local NOK
- acquire trace21ka temperature data, default 1
capture_temperature=0
capture_rainfall=1
- na fill method 0,1, 2 or 4
- method 4 resample
- method 2 is very slow
hadcm3_60ka_seafill_method=1
- if method 4, sample cols and rows
- samplecols1=50
- samplerows1=25
samplecols1=75
samplerows1=35
- acquire hadcm3 60 kyr data
- Warning: you must edit procedure to acquire right dataser
capture_hadcm3_60ka_data=1
- acquire trace 26k data
capture_trace26k_data=0
- acquire glac1d elevation data,
- default 2
- 1 Tarasov glac1d (26 ka)
- 2 Peltier 1ce 6G-D (122 ka) NOTE: must use, when age over 26ka!
capture_elevation=2
- generate rain shadows, assume westerly wind by default
make_rain_shadows=1
calculate_distance=0
- either global_ds or normal_ds must set 1
- normal etopo1 area downscaling normal_ds=3
-
- 1: temperature downscaling: normal_ds=1
- 3: rainfall downscaling: normal_ds=3
normal_ds=3
- 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=0
- 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=1
kustomorofilee1="./predata/dordogne.nc"
- kustomorofilee1="./predata/europe30.nc"
- year to render
fage=40740
- fage=13000
- fage=44950
- fage=38000
- month of year to render
- month 0 sum of all monnths: tex annual rainfall
- month 13 average of all months
fmonth=0
- age in tarasov ice sheet
fage2=fage
- fage=56000
- number of years to average, since "fage"
numyears=20
- sea level, below current: eq 120 means height of -120 to current
- not needed set by default, calculated from age w/ thisscript
sealevel=-70
out_text3="X"
out_text4= "Europe, "
if(capture_temperature==1)
{
out_text_1="TS"
out_text_2="°C"
if(fmonth==7) {
out_text_3=paste0("Temperature of July, ",out_text4)
}
if (fmonth==1) {
out_text_3=paste0("Temperature of January, ",out_text4)
}
}
if(capture_rainfall==1)
{
out_text_1="Precip"
out_text_2="mm"
if(fmonth==7) {
out_text_3=paste0("Precipitation of July, ",out_text4)
}
if (fmonth==1) {
out_text_3=paste0("Precipitation of July, ",out_text4)
}
if(fmonth==0) {
out_text_3=paste0("Annual precipitation , ",out_text4)
}
}
out_text_3=paste0(out_text_3,toString(fage))
out_text_3=paste0(out_text_3, " years ago")
- out_text_1="Precipitation"
- out_text_2="mm"
- out_text_3="Annual precipitation, Europe, "
- out_text_3=paste0(out_text_3,toString(fage))
- out_text_3=paste0(out_text_3, " years ago")
- area to render, lon1, lon2 lat1 lat2 rectangle
- europe
lon1=-15.0
lon2=40.0
lat1=30.0
lat2=70.0
- beringia
- lon1=140
- lon2=240
- lat1=50.0
- lat2=80.0
- east asia
- lon1=80
- lon2=200
- lat1=40.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
if(fage>25999)
{
- NOTE: must use, when age over 26ka!
if(capture_elevation==1)
{
capture_elevation=2
}
}
if(normal_ds==1)
{
if(capture_rainfall==1)
{
normal_ds=3
make_rain_shadows=1
}
}
filename1="KKK"
tslocation1=-999
elevlocation1=-999
- calculate location in tarasov ice sheet data
elevlocation1<<-round( (fage2-26000)/100)*-1
variaabeli="TS"
unitti="Celcius"
lonkkunimi="Temperature"
lonvar1="lon"
latvar1="lat"
- hadcm3b 60ka files path
- hadbasepath<-"D:/varasto_iceagesimu"
hadmonths=1
hadoperaatio=1 ## 1: tex tjuly or precipjuly, 3 annual rain
hadfilename="test.nc"
hadbaseyear=0
- hadvariaabeli="pr"
hadvariaabeli="tas"
unitti="Celcius"
hadlonkkunimi="Near-Surface Air Temperature"
hadlonvar1="longitude"
hadlatvar1="latitude"
hadmonth=fmonth
hadyear=fage
hadyears=numyears
haditem=15006 ## !!! will override in many cases!
## hadcm3b 60ka files path
hadbasepath<<-"D:/varasto_iceagesimu"
hadbaseyear=-1
hadprocesspath<-"./data_processing/"
lones1=0
latis1=0
inputdataname1=""
- non control vars
capture_hadcm3_temperature=0
capture_hadcm3_rainfall=0
- input data dirs
- etopo1 dir
predata="./predata/"
- Trace21ka TS dir
predata2="./predata2/"
dir.create("./data_processing")
dir.create("./plotz")
- outpath11<-"./data_processing/area.nc"
- outpath12<-"./data_processing/area_neworog.nc"
- glac1D HDC file
inpath_tarasov1<-"./predata/TOPicemsk.GLACD26kN9894GE90227A6005GGrBgic.nc"
- peltier ice thickness 122 kyr
inpath_peltier1<-"./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<-"./predata/ETOPO1_Ice_c_gdal.grd"
inpath_etopo2<- "./predata/etopo1_360.nc"
data_processing="./data_processing/"
invarfname1<-"./data_processing/tempin.nc"
invarfname2<-"./data_processing/precipin.nc"
inorofname2<-"./data_processing/area.nc"
inorofname3<-"./data_processing/area_neworog.nc"
outorofname1<-"./data_processing/area_glacial.nc"
outlandseamaskfname1<-"./data_processing/landsea_glacial.nc"
outlandseamaskfname2<-"./data_processing/landsea_glacial.gif"
outvarfname1<-"./data_result/result.nc"
outvarfname1b<-"./data_result/result_land.nc"
outvarfname2<-"./data_result/result2.nc"
outvarfname_accurate1<-"./data_result/accurate.nc"
inorofname360_1<-"./data_processing/oroin.nc"
invarfname360_1<-"./data_processing/tempin.nc"
inicefname360_1<-"./data_processing/maskin.nc"
inorofname180_1<-"./data_processing/oroin_180.nc"
invarfname180_1<-"./data_processing/tempin_180.nc"
inicefname180_1<-"./data_processing/maskin_180.nc"
invarfname180_2<-"./data_processing/precipin_180.nc"
smallrainame1<-"./data_processing/smallrain.nc"
- outname1<-"./data_processing/chelsa_current_annual_rain.nc"
- outname2<-"./data_processing/chelsa_lgm_ccsm4_annual_rain.nc"
outname3<-"./data_processing/etopo1.nc"
outname4<-"./data_processing/lons.nc"
outname5<-"./data_processing/lats.nc"
outname6<-"./data_processing/distance.nc"
outname7<-"./data_processing/slope.nc"
outname8<-"./data_processing/aspect.nc"
outname9<-"./data_processing/hillshade.nc"
outname10<-"./data_processing/tpi.nc"
outname11<-"./data_processing/westness.nc"
outname12<-"./data_processing/blurelev.nc"
outname13<-"./data_processing/windir.nc"
outname14<-"./data_processing/etopo_big.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"
- 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;
}
- subroutines
sealevel_from_age<-function(inage)
{
- agetable1<-c(60000,24000,22000,20000,18000,16000, 15000,14000,13000,12000,11000, 10000,8000,7000,6000 ,0)
- leveltable1<-c(-65,-125,-130,-125,-120,-115, -110,-80,-75,-65,-60, -45,-15,-10,-5 ,0 )
agetable1<-c(110000, 65000,60000,55000,50000,45000, 40000,37500, 35000, 30000, 24000,22000,20000,18000,16000, 15000,14000,13000,12000,11000, 10000,8000,7000,6000 ,0)
leveltable1<-c(-40, -85,-75, -60,-70,-80,-75, -65, -80, -75, -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_temperature<<-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
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,
zname="Zname", longname=longnamex, xname="lon", yname="lat")
}
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)
}
generate_hadfilename<-function(hadbaspath1, yrr1, varr1)
{
#hadfilename="D:/varasto_iceagesimu/bias_regrid_tas_55_57.5kyr.nc"
hadfilenamex1=hadbaspath1
hadfilenamex1<-paste0(hadfilenamex1,"/bias_regrid_")
hadfilenamex1<-paste0(hadfilenamex1,varr1)
hadfilenamex1<-paste0(hadfilenamex1,"_")
yrr1=yrr1-(2500/2)+1
a=yrr1/100
b=a
c=round(b/25,0)
d=(c*25)/10
e=d+2.5
hadbaseyear<<-d*1000
#print (b)
#print (c)
#print (d)
#print (e)
hadfilenamex1<-paste0(hadfilenamex1,toString(d))
hadfilenamex1<-paste0(hadfilenamex1,"_")
hadfilenamex1<-paste0(hadfilenamex1,toString(e))
hadfilenamex1<-paste0(hadfilenamex1,"kyr.nc")
return(hadfilenamex1)
}
- 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")
}
tracelocation3 <- function(dirra1, variaabeli, fage, numyears, fmonth)
{
aage2=fage-numyears
aage1=fage
lista0<-list.files(path=dirra1,pattern="trace*", ignore.case=TRUE)
pitu<-length(lista0)
- print(fage1,)
#print("Kup")
#print (pitu)
# lista0
fagek1=as.integer(aage1)
fagek2=as.integer(aage2)
if(fagek2>fagek1)
{
huba=fagek1
fagek1=fagek2
fagek2=huba
}
mara1=-999
mara2=-999
for (n in 1:(pitu))
{
row1<-lista0[n]
jono1<-row1
#print (jono1)
s1<-gsub('\\.', '_', jono1[1])
s2<-gsub('\\-', '_', s1)
s3<-gsub('BP', , s2)
q1 <-str_split(s3, "_")
sr1<-q11
age10<-as.integer(sr1[3])
age20<-as.integer(sr1[4])
age1<-age10[1]
age2<-age20[1]
if (fagek2==400)
{
fagek2<-0
}
if(fagek1<age1)
{
if(fagek1>age2)
{
mara1<-n
}
}
if(fagek1==age1)
{
mara1<-n
# break
}
if(fagek2<age1)
{
if(fagek2>age2)
{
mara2<-n
break
}
}
if(fagek2==age1)
{
mara2<-n
}
- print (paste0(str(age1)," ",str(age2)))
}
print("Loop done")
print(mara1)
print(mara2)
- Haista
ivars0=0
ivars1=0
mara=0
natta=0
agari1=0
- jn warning debug
for (mara in mara1:mara2)
{
row1<-lista0[mara]
print(row1)
jono1<-row1
s1<-gsub('\\.', '_', jono1[1])
s2<-gsub('\\-', '_', s1)
s3<-gsub('BP', , s2)
q1 <-str_split(s3, "_")
sr1<-q11
age1<-as.numeric(sr1[3])
age2<-as.numeric(sr1[4])
# print (mara)
filename1<-paste0(predata2, row1)
print(filename1)
putin1 <- nc_open(filename1)
print(class(putin1))
# stop("k")
# ivars0<-0
ivars0 <- ncvar_get(putin1, variaabeli)
lones1<- ncvar_get(putin1, lonvar1)
latis1<- ncvar_get(putin1, latvar1)
if(natta>0)
{
ivars1<-abind(ivars1,ivars0,along=3)
#ivars1=c(ivars1,ivars0)
}
else
{
agari1<-as.numeric(age1)
ivars1<-ivars0
}
natta=natta+1
nc_close(putin1)
}
elevlocation1<<-round( (fage-26000)/100)*-1
class(ivars1)
dim (ivars1)
# varlocation1<<-((fage-age1)*12*-1)+fmonth
- orig
- load_trace_var3d_3(ivars1, lons1, lats1,variaabeli, unitti,agari1, fage, numyears, fmonth )
load_trace_var3d_3(ivars1, lons1, lats1,variaabeli, unitti,agari1, fage ,numyears, fmonth)
}
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")
print(number2)
elev0<-elev[,,number2]
str(elev0)
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"
1
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)
if(capture_elevation==2)
{
print("Peltier test code")
vali_raster1<-resample( fine_raster1, coarse_raster1)
coarse_raster1<-coarse_raster1+vali_raster1
}
- coarseoro<- resample(p1, coarse_raster)
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")
outras2<-outras1
outras2[outras2[]<=0.0] <- 0
outras2[outras2[]>0.0] <- 1.0
plot(outras2, col=rev(plasma(256)))
writeRaster(outras2, outlandseamaskfname1, overwrite=TRUE, format="CDF", varname="Land", varunit="unit",
longname="Land", xname="lon", yname="lat")
- writeRaster(outras2, outlandseamaskfname2, overwrite=TRUE, format="CDF", varname="Land", varunit="unit",
- longname="Land", 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)
- out_text_1="TS"
- out_text_2="oC"
- out_text_3="Temperature"
writeout(outras2,outvarfname1,out_text_1, out_text_2, out_text_3)
mask1<-raster(outlandseamaskfname1)
outras3<-outras2
outras3<-outras3*mask1
writeout(outras3,outvarfname2,out_text_1, out_text_2, out_text_3)
#writeRaster(outras3,outvarfname2 , overwrite=TRUE, format="GIF", 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
print("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.")
print(orofilee1)
print(tempfilee1)
rasnum1<<- 1
coarse_raster=raster(tempfilee1)
fine_raster=raster(orofilee1)
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<-"./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("Temperature of July, ",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(-15,-10,-5,0,5,10,15,20,25,30,35)
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(-15,35,5)
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(rainbow(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()
}
fill.na <- function(x) {
width1<-67
center = 0.5 + (width1*width1/2)
if( is.na(x)[center] ) {
return( round(mean(x, na.rm=TRUE),0) )
} else {
return( round(x[center],0) )
}
}
get_europe_gtopo30<-function()
{
# france
#ext1 <- extent(-4,8, 39, 51)
#dordogne
#ext1 <- extent(-2,4,42 , 48)
#europe
ext1 <- extent(-15,40,30,70)
outname1<-"./predata/europe30.nc"
r1<-raster("./predata/gt30w020n90.tif")#
r2<-raster("./predata/gt30e020n90.tif") #
r4<-raster("./predata/gt30w020n40.tif") #
r3<-raster("./predata/gt30e020n40.tif")
#plot(r1)
mosaik1 <- merge(r1,r2,r3,r4)
kropped1<-crop(mosaik1, ext1)
#plot(kropped1)
crs(kropped1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
writeRaster(kropped1, filename=outname1, format="CDF", overwrite=TRUE)
}
load_hadcm3_slice <- function(hadfilename, hadvariaabeli, haditem)
{
print("Trying to load hadcm3 filee ")
## varlocation1<<-((fage-age1)*12*-1)+fmonth1
varlocation1=haditem
print(hadfilename)
hadcount=1
putin1 <- nc_open(hadfilename)
print(class(putin1))
lones1<- ncvar_get(putin1, hadlonvar1)
latis1<- ncvar_get(putin1, hadlatvar1)
t <- ncvar_get(putin1, "time")
lenlones1<-length(lones1)
lenlatis1<-length(latis1)
vaari0<- ncvar_get(putin1,hadvariaabeli, start=c(1,1,haditem), count=c(lenlones1,lenlatis1,1) )
nc_close(putin1)
print(dim(vaari0))
padding1 = matrix(0,lenlones1,180)
vaari1<-cbind(padding1,vaari0)
#vaari1<-vaari0
print(dim(vaari1))
vaari1<- apply(t(vaari1),2,rev)
rrvar00<-raster (vaari1)
rrvar00@extent<-extent(0, 360, -90, 90)
width1 = 67
rrvar1<- focal(rrvar00, w = matrix(1,width1,width1), fun = fill.na,
pad = TRUE, na.rm = FALSE)
summary(getValues(rrvar1))
plot(rrvar1, col=rainbow(120))
print("invarfname:")
print(invarfname360_1)
print("Write hadcm3 file ..")
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")
# print("Test break")
# stop()
}
load_hadcm3_year_month <- function(hadfilename, hadvariaabeli, hadyear, hadmonth, hadbeginyear)
{
hadyeardelta=hadyear-hadbeginyear
haddelta=hadyeardelta*12+hadmonth
haditem=30001-haddelta
load_hadcm3_slice(hadfilename, hadvariaabeli, haditem)
}
load_average_hadcm3_years_months <- function(hadfilename, hadvariaabeli,hadbeginyear, hadyear, hadmonth, hadyears, hadmonths, operaatio, metoden1)
{
## operaatio==0: sum
## operaatio==1: averege of data matrices
print ("Loading many HADCM3 slices.")
print(hadfilename)
putin1 <- nc_open(hadfilename)
lones1<- ncvar_get(putin1, hadlonvar1)
latis1<- ncvar_get(putin1, hadlatvar1)
t <- ncvar_get(putin1, "time")
lenlones1<-length(lones1)
lenlatis1<-length(latis1)
plastic1<-matrix(0,lenlones1,lenlatis1)
kohtia=0
hadyearend=hadyear+hadyears
hadmonthend=hadmonth+hadmonths
#print ("Hadyear")
#print (hadyear)
#print (hadyearend)
#print (hadmonth)
#print (hadmonthend)
hadyeardelta1=hadyear-hadbeginyear
haddelta1=hadyeardelta1*12
haditem1=30001-haddelta1
hadpususize1=hadyears*12
#print ("Haddeltas ")
#print( hadyeardelta1)
#print(haddelta1)
#print(haditem1)
print("Get had data ...")
hadpusu1<-ncvar_get(putin1,hadvariaabeli, start=c(1,1,haditem1), count=c(lenlones1,lenlatis1,hadpususize1) )
nc_close(putin1)
#print (hadpusu1[7])
hadex1=0
ammax=hadyears-1
kammax=hadmonth+hadmonths-1
#print("Extract haddata:")
hadpusupitu1=length(hadpusu1)
#print (hadpusupitu1)
#print (hadpususize1)
#print (ammax)
#print (kammax)
#print (hadmonth)
#print ("luup")
vuozia1=0
for (m in (1:ammax) )
{
#print("R")
for (n in (hadmonth:kammax) )
{
#print(".")
#print(n)
hadex1=m*12+n
# print (hadex1)
vaari0<-hadpusu1[,,hadex1]
plastic1<-plastic1+vaari0
kohtia=kohtia+1
}
vuozia1=vuozia1+1
}
print("Kohtia:")
print(kohtia)
print("Vuozia1:")
print (vuozia1)
print("Operaatio:")
print(operaatio)
if(operaatio==0)
{
## operaatio 0, sum
## one year precip
print("sum")
vaari3<-plastic1
}
if(operaatio==1)
{
print("one month, many years")
## one month, many years
vaari3<-(plastic1/kohtia)
# JN WARNING koe
}
if(operaatio==2)
{
print("annualaverage, one year")
## annual
#print("Operaatio 3: temperature annual")
vaari3<-plastic1/12
}
if(operaatio==3)
{
print("annual sum, many years")
#print("Operaatio 3: rainfall annual")
vaari3<-plastic1/vuozia1
}
# plot(vaari3)
vaari4<-vaari3
vaari5<-vaari3
vaari6<-vaari3
vaari7<-vaari3
dima1<-dim(vaari3)
w1=dima1[1]
h1=dima1[2]
print (w1)
print (h1)
if(metoden1==2)
{
for (ix in (1:w1) )
{
beg1=1
viksu1=-1
for (iy in (1:h1) )
{
piksu1=vaari3[ix,iy]
if(is.na(piksu1))
{
if(beg1==0)
{
vaari4[ix,iy]=viksu1
vaari3[ix,iy]=1
}
else
{
vaari4[ix,iy]=piksu1
}
}
else
{
beg1=0
viksu1=piksu1
beg1=0
}
}
beg1=1
viksu1=-1
for (iy in seq(h1, 1, by=-1 ) )
{
piksu1=vaari3[ix,iy]
if(is.na(piksu1))
{
vaari4[ix,iy]=viksu1
}
else ## ei na
{
beg1=0
viksu1=piksu1
beg1=0
}
}
}
- sekond
for (iy in (1:h1) )
{
beg1=1
viksu1=-1
for (ix in (1:w1) )
{
piksu1=vaari5[ix,iy]
if(is.na(piksu1))
{
if(beg1==0)
{
vaari6[ix,iy]=viksu1
vaari5[ix,iy]=1
}
else
{
vaari6[ix,iy]=piksu1
}
}
else
{
beg1=0
viksu1=piksu1
beg1=0
}
}
beg1=1
viksu1=-1
for (ix in seq(w1, 1, by=-1 ) )
{
piksu1=vaari5[ix,iy]
if(is.na(piksu1))
{
vaari6[ix,iy]=viksu1
}
else ## ei na
{
beg1=0
viksu1=piksu1
beg1=0
}
}
}
vaari7<-(vaari4+vaari6)/2
}
- plot(vaari3)
#print(dim(vaari3))
padding1 = matrix(0,lenlones1,180)
vaari1<-cbind(padding1,vaari7)
#vaari1<-vaari0
#print(dim(vaari1))
vaari1<- apply(t(vaari1),2,rev)
rrvar00<-raster (vaari1)
rrvar00@extent<-extent(0, 360, -90, 90)
- fill nas with nearest non-na
if(metoden1==1)
{
## orig 67
width1 = 67
##width1=157
rrvar1<- focal(rrvar00, w = matrix(1,width1,width1), fun = fill.na,
pad = TRUE, na.rm = FALSE)
}
if(metoden1==0)
{
rrvar0<-rrvar00
rrvar1<-rrvar0
}
- summary(getValues(rrvar1))
plot(rrvar1, col=rainbow(8))
print("invarfname:")
print(invarfname360_1)
print("Write hadcm3 file ..")
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")
print("HadCM3 data file loaded .")
}
- 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
}
- generate rain shadows, assume westerly wind by default
generate_rain_shadows<-function()
{
## params
hillshade_slope=2
hillshade_aspect=260
## direction of global, assumed wind
tirektion=(260*6.28)/360
orography1<-raster(outorofname1)
## copy raster
windir<-orography1
windir<-windir*0
windir<-windir+tirektion
latras <- lonras <- orography1
xy <- coordinates(orography1)
lonras[] <- xy[, 1]
latras[] <- xy[, 2]
#print("write windir ...")
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")
#print(" write lonras")
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")
slope <- terrain(orography1, opt='slope')
aspect <- terrain(orography1, opt='aspect')
hill <- hillShade(slope, aspect, hillshade_slope, hillshade_aspect)
#plot(hill, col=grey(0:100/100), legend=FALSE, main='France')
tpi5 <- tpiw(orography1, 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(orography1, 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")
#plot(rasa1)
if(calculate_distance==1)
{
cacheras1<-orography1
cacheras1[cacheras1 > 1] <- NA
cacheras1[cacheras1 < 0] <- 1
print( "Calculating distance. Maybe very slooow. \n Wait ...")
distrasa1 <- distance(cacheras1)
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")
}
}
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 <- 3 # Minimum number of iterations
max_iter <- 7 # Maximum number of iterations
p_train <- samplerate # Subsampling of the initial data
print("Dissever() begin ...")
oma_juttu <- dissever(coarse = coarse_raster, fine = raster_stack, method = dismethod, p = p_train, min_iter = min_iter,max_iter = max_iter, verbose=1)
print("Dissever() end.")
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)
}
downscale_rainfall<-function()
{
print("Downscaling rainfall ...")
#outname3<-"./data_processing/etopo1.nc"
#outname4<-"./data_processing/lons.nc"
#outname5<-"./data_processing/lats.nc"
#outname6<-"./data_processing/distance.nc"
#outname7<-"./data_processing/slope.nc"
#outname8<-"./data_processing/aspect.nc"
#outname9<-"./data_processing/hillshade.nc"
#outname10<-"./data_processing/tpi.nc"
#outnam44e11<-"./data_processing/westness.nc"
#outname12<-"./data_processing/blurelev.nc"
#outname13<-"./data_processing/windir.nc"
#outname14<-"./data_processing/etopo_big.nc"
etopo0<<-raster(outorofname1)
aspect0<<-raster(outname8)
tpi0<<-raster(outname10)
lonx0<<-raster(outname4)
hillshade0<<-raster(outname9)
slope0<<-raster(outname7)
windir0<<-raster(outname13)
bluretopo0<<-raster(outname12)
latx0<<-raster(outname5)
#print("Krop smallrain")
ext1<-extent(etopo0)
#print(ext1)
#smallrain00<-raster(invarfname180_1)
smallrain00<-raster(invarfname180_2)
smallrain0<-crop(smallrain00, ext1)
plot(smallrain0, col=rev(viridis(100)))
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
#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)
dsobject=1
if(dsobject==1)
{
print("DS 1")
# preprocess_rasters=1
plot(etopo0)
g0<-etopo1
g1<-etopo1
#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
# rastafari1<-stack(etopo0, ebhk1) ## kohtal
rastafari1<-stack(g0,g1) ## kohtal
print(names(rastafari1))
#plot(rastafari1$Elevation)
#out3a<<-downscale_dissever(smallrain0, rastafari1,"lm",1.0)
out3<<-downscaler(smallrain0,g0,0)
out3[out3 < 1] <- 1
plot(out3, col=rev(viridis(100)))
print(minValue(out3))
print(maxValue(out3))
#writeout(out3,outvarfname1,"Rainfall", "mm/yr", "IceAgeRain")
rainout1<-"./data_result/outrain1.nc"
#writeout(out3,rainout1,"Rainfall", "mm/yr", "IceAgeRain")
# out_text_1="Rainfall"
# out_text_2="mm/yr"
# out_text_3="IceAgeRain"
writeout(out3,rainout1,out_text_1, out_text_2, out_text_3)
}
}
get_dordogne_srtm<-function()
{
## France/Dorgogne area SRTM netcdf file creation
## "R" script
srtm1 <- getData('SRTM', lon=-5, lat=40)
srtm2 <- getData('SRTM', lon=-5, lat=45)
srtm3 <- getData('SRTM', lon=-5, lat=50)
srtm4 <- getData('SRTM', lon=-5, lat=55)
srtm6 <- getData('SRTM', lon=0, lat=40)
srtm7 <- getData('SRTM', lon=0, lat=45)
srtm8 <- getData('SRTM', lon=0, lat=50)
srtm9 <- getData('SRTM', lon=0, lat=55)
srtm10 <- getData('SRTM', lon=5, lat=40)
srtm11 <- getData('SRTM', lon=5, lat=45)
srtm12 <- getData('SRTM', lon=5, lat=50)
srtm13 <- getData('SRTM', lon=5, lat=55)
srtm14 <- getData('SRTM', lon=10, lat=40)
srtm15 <- getData('SRTM', lon=10, lat=45)
srtm16 <- getData('SRTM', lon=10, lat=50)
srtm17<- getData('SRTM', lon=10, lat=55)
mosaik1 <- merge( srtm1, srtm2,srtm3, srtm4, srtm6,srtm7,srtm8, srtm9,srtm10, srtm11,srtm12,srtm13,
srtm14, srtm15,srtm16,srtm17)
#dordogne
ext1 <- extent(-2,4,42 , 48)
kropped1<-crop(mosaik1, ext1)
#plot(kropped1)
crs(kropped1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
writeRaster(kropped1, filename='./predata/dordogne.nc', format="CDF", overwrite=TRUE)
system("del srtm*.*")
}
hadcm3_loadslice <- function(temp_name, var_name, hadyear)
{
putin1 <- nc_open(temp_name)
lones1<<- ncvar_get(putin1, "longitude")
latis1<<- ncvar_get(putin1, "latitude")
t <- ncvar_get(putin1, "time")
lenlones1<-length(lones1)
lenlatis1<-length(latis1)
deltayears1=hadyear-hadbaseyear
deltamonths1=deltayears1*12
item1=30000-deltamonths1-12+1
months1=12
temp_pusu1<-ncvar_get(putin1,var_name, start=c(1,1,item1), count=c(lenlones1,lenlatis1,months1) )
nc_close(putin1)
taimi1=t[1]
return(temp_pusu1)
}
generate_hadfilename<-function(hadbaspath1, yrr1, varr1)
{
hadfilenamex1=hadbaspath1
hadfilenamex1<-paste0(hadfilenamex1,"/bias_regrid_")
hadfilenamex1<-paste0(hadfilenamex1,varr1)
hadfilenamex1<-paste0(hadfilenamex1,"_")
a=as.integer(yrr1/2500)
b=a*2500
c=b/1000
d=c+2.5
hadbaseyear<<-b
hadfilenamex1<-paste0(hadfilenamex1,toString(c))
hadfilenamex1<-paste0(hadfilenamex1,"_")
hadfilenamex1<-paste0(hadfilenamex1,toString(d))
hadfilenamex1<-paste0(hadfilenamex1,"kyr.nc")
return(hadfilenamex1)
}
load_had_slices<-function(beginyr1, yrs1, varr1)
{
endyr1=beginyr1+yrs1-1
print("Loading haccm3 slices, wait ...")
markki1=0
yyyy1=0
for (yrr1 in (beginyr1:endyr1))
{
hadfilename=generate_hadfilename(hadbasepath, yrr1, varr1)
print(yrr1)
print (hadfilename)
slice00=hadcm3_loadslice(hadfilename, varr1, yrr1)
if(markki1==0)
{
baseslice1<-slice00
}
else
{
# add slices
baseslice1<-baseslice1+slice00
}
markki1=1
yyyy1=yyyy1+1
}
#print(head(baseslice1))
baseslice1=baseslice1/yyyy1
return (baseslice1)
}
draw_climate_diagram<-function(lampot, sadem)
{
#mydata <- read.csv("kiova2.txt", header=FALSE, sep=";")
labeli='Paris, 40750 BP'
nimi="paris_40750bp"
datanimi=paste(nimi,".txt");
kuvanimi=paste(nimi,".svg");
prmax=100
prmin=0
tmax = 20.0
tmin=-25.0
tstep=5
widthi=10
heighti=16
asteikko<-c(" "," ","3"," "," ","6"," ", " ","9"," "," ","12" )
svg(kuvanimi, width=widthi, height=heighti)
deltapr=prmax-prmin
deltatee<-(tmax-tmin)
ratio<-deltapr/deltatee
y2offset= -1*ratio*tmin
total_sadem=sum(sadem)
avg_lampotila=sum(lampot)/12
avg_lampotila=(round(avg_lampotila)*10)/10
max_lampotila=max(lampot)
min_lampotila=min(lampot)
par(mar=c(6,6,6,6),cex.axis=2,cex.lab=2.5)
b<-barplot(sadem, names.arg=asteikko,
col="blue", border="blue",ylim=c(prmin, prmax),
cex.axis=2.5, cex.names=2.5 )
lines(b, (lampot*ratio)+y2offset, col="Red",lwd=8)
right.axis.ticks<- seq(from =tmin , to=tmax , by=tstep)
axis(4,at=(right.axis.ticks*ratio)+y2offset,labels=paste0(right.axis.ticks),las=2,
cex.axis=2.5)
mtext(side = 2, line = 3, 'Precipitation', cex=2.5, col="darkblue")
mtext(side = 4, line = 3, 'Temperature', cex=2.5, col="darkred")
mtext(side = 1, line = 3, 'Month', cex=2.5, col="darkgreen")
text(7,(prmax-2),cex=3.5, labeli);
text(1,(prmax-8),cex=2.4, pos=4, paste("Tavg=",avg_lampotila, " C" ));
text(1,(prmax-12),cex=2.4, pos=4, paste("Tamax=",max_lampotila, " C" ));
text(1,(prmax-16),cex=2.4, pos=4,paste("Tamin=",min_lampotila, " C" ));
text(1,(prmax-20),cex=2.4, pos=4,paste("Pra=",total_sadem, " mm" ));
}
get_had_climate_data<-function(beginyears, years, targetname1, lat1, lon1)
{
print("Loading data, wait ...")
varr1="tas"
varr2="pr"
tempsit1<-load_had_trapezoid(beginyears, years, varr1, lon1, lat1)
precsit1<-load_had_trapezoid(beginyears, years, varr2, lon1, lat1)
#print (tempsit1)
#print (precsit1)
months1<-1:12
tempsit1<-round(tempsit1, digits = 1)
precsit1<-round(precsit1, digits = 1)
tavg1<-sum(tempsit1)/12.0
pannual1<-sum(precsit1)
df1<-data.frame(months1, tempsit1,precsit1)
names(df1)<-c("Month", "T", "P")
coutname1=paste0(targetname1, ".csv")
#write.csv2(df1,coutname1)
#write.table(df1,file=coutname1,sep=";")
write.table(df1,file=coutname1,sep=";",row.names=FALSE)
print("Monthly data:")
print(df1)
print ("Climate averages:")
print (tavg1)
print (pannual1)
draw_climate_diagram(tempsit1, precsit1)
}
had_twoslicer<-function(beyr1,yrs1,varr1)
{
enyr1=beyr1+yrs1
print(beyr1)
print(enyr1)
hadbasepath1<<-hadbasepath
hadnames1<-vector(mode="character", length=2)
hadbaseyears<-rep(0,2)
#print (hadbaseyears)
hadnames1[1]<-generate_hadfilename(hadbasepath1, beyr1, varr1)
hadbaseyears[1]<-hadbaseyear
hadnames1[2]<-generate_hadfilename(hadbasepath1, enyr1, varr1)
hadbaseyears[2]<-hadbaseyear
deltayears1=(beyr1-hadbaseyears[2])*-1
deltamonths1=deltayears1*12
item1=30000-deltamonths1-12+1
months1=12
deltayears2=enyr1-hadbaseyears[2]
deltamonths2=deltayears2*12
item2=30000-deltamonths2-12+1
months2=12
- print (hadnames1[1])
- print (hadnames1[2])
- print(deltayears1)
- print(deltamonths1)
- print(deltayears2)
- print(deltamonths2)
twoo1=0
if(hadbaseyears[1]==hadbaseyears[2])
{
- print("Twoo 1")
twoo1=1
}
if(deltayears1>-1)
{
twoo1=1
}
if(twoo1==0)
{
putin1 <- nc_open(hadnames1[1])
lones1<<- ncvar_get(putin1, "longitude")
latis1<<- ncvar_get(putin1, "latitude")
t <- ncvar_get(putin1, "time")
lenlones1<-length(lones1)
lenlatis1<-length(latis1)
temp_pusu1<-ncvar_get(putin1,varr1, start=c(1,1,item1), count=c(lenlones1,lenlatis1,deltamonths1) )
nc_close(putin1)
- print("put in 2")
putin2 <- nc_open(hadnames1[2])
lones1<<- ncvar_get(putin2, "longitude")
latis1<<- ncvar_get(putin2, "latitude")
t2 <- ncvar_get(putin2, "time")
lenlones1<-length(lones1)
lenlatis1<-length(latis1)
temp_pusu2<-ncvar_get(putin2,varr1, start=c(1,1,item2), count=c(lenlones1,lenlatis1,deltamonths2) )
nc_close(putin2)
# print("put in 2")
dima1=dim(temp_pusu2)
pusu3=abind(temp_pusu1,temp_pusu2,along=3)
} #two file buffers
else
{
# print("Twoo 1 ...")
deltayears2=beyr1-hadbaseyears[1]
deltamonths2=deltayears2*12
item2=30000-deltamonths2-12+1
months2=(enyr1-beyr1)*12
#print(deltayears2)
#print(deltamonths2)
#print(item2)
#print("put in 2")
putin2 <- nc_open(hadnames1[2])
lones1<<- ncvar_get(putin2, "longitude")
latis1<<- ncvar_get(putin2, "latitude")
t2 <- ncvar_get(putin2, "time")
lenlones1<-length(lones1)
lenlatis1<-length(latis1)
pusu3<-ncvar_get(putin2,varr1, start=c(1,1,item2), count=c(lenlones1,lenlatis1,months2) )
nc_close(putin2)
}
dima3=dim(pusu3)
#print(dima1)
#print(dima3)
as1<- array(rep(0, 720*180*12), dim=c(720, 180, 12))
ylimit1=dima3[3]
#print (dim(as1))
hhh1=0
maxima1<-(ylimit1/12)-1
print (maxima1)
for( m in 1:maxima1)
{
for( n in 1:12)
{
has1<-pusu3[,,m*12+n]
as1[,,n]<-as1[,,n]+pusu3[,,m*12+n]
}
hhh1=hhh1+1
}
as1<-as1/hhh1
return(as1)
}
load_had_trapezoid<-function(beginyr1, yrs1, varr1, lon1, lat1)
{
- slice00=load_had_slices(beginyr1, yrs1, varr1)
slice00<-had_twoslicer(beginyr1,yrs1,varr1)
dima1=dim(slice00)
#print (dima1)
max1=dima1[1]
may1=dima1[2]
londex2=which(lones1 >= lon1 )[1]
latdex2=which(latis1 >= lat1 )[1]
londex1=londex2-1
latdex1=latdex2-1
if(londex1<1) londex1=max1
if(latdex1<1) latdex1=may1
abslon1=lones1[londex1]
abslat1=latis1[latdex1]
abslon2=lones1[londex2]
abslat2=latis1[latdex2]
#print("lons")
#print(abslon1)
#print(abslon2)
#print(abslat1)
#print(abslat2)
#print (max1)
#print (may1)
#print (lones1[0])
vektor1<-1:12
vektor1<-vektor1*0
n=7
for (n in 1:12)
{
## attempt to process trapezoid
value1=slice00[londex1,latdex1, n]
value2=slice00[londex1,latdex2, n]
value3=slice00[londex2,latdex1, n]
value4=slice00[londex2,latdex2, n]
rulon1=abslon2-abslon1
rulat1=abslat2-abslat1
daata1<-c(value1,value2,value3,value4)
matrix <- matrix(daata1, nrow=2, ncol=2)
r <- raster(matrix)
## lon lat
extent(r) <- c(abslon1, abslon2, abslat1,abslat2)
## reso 100x100
s <- raster(nrow=100, ncol=100)
extent(s)<-extent(r)
s <- resample(r, s, method='bilinear')
xy <- cbind(lon1,lat1)
resultt1<-extract(r, xy)
vektor1[n]=resultt1
}
return(vektor1)
}
load_had_raster<-function(beginyr1, yrs1, varr1, month1)
{
#slaici1=load_had_slices(beginyear1, yrs1, varr1)
slaici1<-had_twoslicer(beginyr1,yrs1,varr1)
dima1=dim(slaici1)
print (dima1)
if(month1==0)
{
markki=0
yyyy1=0
## select all months
for (n in 1:12)
{
vaari0=slaici1[,,month1]
if(markki1==0)
{
baseslice1<-slice00
}
else
{
# add slices
baseslice1<-baseslice1+slice00
}
merkki=1
yyyy1=yyyy1+1
}
vaari0=baseslice1/yyyy1
}
else
{
vaari0=slaici1[,,month1]
}
print (dim(vaari0))
padding1 = matrix(0,720,180)
vaari1<-cbind(padding1,vaari0)
vaari1<- apply(t(vaari1),2,rev)
rrvar1<-raster (vaari1)
rrvar1@extent<-extent(0, 360, -90, 90)
crs(rrvar1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
rvarfilename1=paste0(hadprocesspath, "global_360_",varr1,"_",month1,"-nc")
longname1=paste0(varr1," ",toString(beginyr1) )
writeRaster(rrvar1, rvarfilename1, overwrite=TRUE, format="CDF", varname=varr1, varunit="unit",
longname=longname1, xname="lon", yname="lat")
}
load_had_rasters_var<-function(byr1, yrs1, varr1)
{
yrmid1=byr1+(yrs1/2)
slaici1<-had_twoslicer(byr1,yrs1,varr1)
dima1=dim(slaici1)
print (dima1)
for (n in 1:12)
{
print (n)
vaari0=slaici1[,,n]
padding1 = matrix(0,720,180)
vaari1<-cbind(padding1,vaari0)
vaari1<- apply(t(vaari1),2,rev)
rrvar1<-raster (vaari1)
rrvar1@extent<-extent(0, 360, -90, 90)
crs(rrvar1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
if(n==1)
{
rs1=stack(rrvar1)
}
else
{
rs1=stack(rs1, rrvar1)
}
}
plot(rs1)
rvarfilename1=paste0(hadprocesspath, "global_360_",varr1,"_",yrmid1)
longname1=paste0(varr1," ",toString(yrmid1) )
writeRaster(rs1, rvarfilename1, overwrite=TRUE, format="CDF", varname=varr1, varunit="unit",
longname=longname1, xname="lon", yname="lat")
inputdataname1<<-rvarfilename1
}
load_climate<-function()
{
- beginyear1=40700
beginyear1=33950
years1=100
month1=7
varr1="tas"
varr2="pr"
- paris
- targetname1="paris"
- targetlat1=48.856667
- targetlon1=2.351111
targetname1="sungir"
targetlat1=56.175833
targetlon1=40.509167
- targetlon1=0.0
print("-----------------------------")
print("Age:")
hage1<-beginyear1+(years1/2)
print(hage1)
print("Target:")
print(targetname1)
print (targetlon1)
print (targetlat1)
print ("")
get_had_climate_data(beginyear1, years1,targetname1,targetlat1, targetlon1)
}
hadcm3_raster_fillna_method_2<-function(vaari3)
{
printf(" FillNA method 2. Sooow, wait patiently ...")
vaari4<-vaari3
vaari5<-vaari3
vaari6<-vaari3
dima1<-dim(vaari3)
w1=dima1[1]
h1=dima1[2]
#print (w1)
#print (h1)
for (ix in (1:w1) )
{
beg1=1
viksu1=-1
for (iy in (1:h1) )
{
piksu1=vaari3[ix,iy]
if(is.na(piksu1))
{
if(beg1==0)
{
vaari4[ix,iy]=viksu1
vaari3[ix,iy]=1
}
else
{
vaari4[ix,iy]=piksu1
}
}
else
{
beg1=0
viksu1=piksu1
beg1=0
}
}
beg1=1
viksu1=-1
for (iy in seq(h1, 1, by=-1 ) )
{
piksu1=vaari3[ix,iy]
if(is.na(piksu1))
{
vaari4[ix,iy]=viksu1
}
else ## ei na
{
beg1=0
viksu1=piksu1
beg1=0
}
}
}
- sekond
for (iy in (1:h1) )
{
beg1=1
viksu1=-1
for (ix in (1:w1) )
{
piksu1=vaari5[ix,iy]
if(is.na(piksu1))
{
if(beg1==0)
{
vaari6[ix,iy]=viksu1
vaari5[ix,iy]=1
}
else
{
vaari6[ix,iy]=piksu1
}
}
else
{
beg1=0
viksu1=piksu1
beg1=0
}
}
beg1=1
viksu1=-1
for (ix in seq(w1, 1, by=-1 ) )
{
piksu1=vaari5[ix,iy]
if(is.na(piksu1))
{
vaari6[ix,iy]=viksu1
}
else ## ei na
{
beg1=0
viksu1=piksu1
beg1=0
}
}
}
vaari7<-(vaari4+vaari6)/2
return(vaari7)
}
fill.na <- function(x) {
width1<-67
center = 0.5 + (width1*width1/2)
if( is.na(x)[center] ) {
return( round(mean(x, na.rm=TRUE),0) )
} else {
return( round(x[center],0) )
}
}
fill.na7 <- function(x) {
width1<-7
center = 0.5 + (width1*width1/2)
if( is.na(x)[center] ) {
return( round(mean(x, na.rm=TRUE),0) )
} else {
return( round(x[center],0) )
}
}
get_hadgcm3_global_raster<-function(byr1, yrs1, var1, month1, varunit1, ovarnam1, ofisuffix1, method1)
{
- mont 0 sum of all year
- month 13 average of all year
- method 4: resample ratser in input, attempting
- toreduce original downscaling out
print ("Get global raster ...")
load_had_rasters_var(byr1, yrs1, var1)
inputdataname10<-paste0(inputdataname1,".nc")
stak1<-stack(inputdataname10)
if(month1==0)
{
r00 <- calc(stak1, sum)
}
if(month1==13)
{
r00 <- calc(stak1, sum)/12
}
else
{
if(month1>0)
{
r00<-stak1@layersmonth1
}
}
oname1=paste0(hadprocesspath,"/",ofisuffix1,"0_360.nc")
writeout(r00, oname1, var1, "C", var1)
rr00_180<-rotate(r00)
## rr00_180@extent<-extent(-180, 180, -90, 90)
plot(rr00_180, col=rainbow(64))
oname2=paste0(hadprocesspath,"/",ofisuffix1,"0_180.nc")
writeout(rr00_180, oname2, var1, "C", var1)
inras1<-raster(oname1)
inras2<-raster(oname2)
#varunit1="C"
if(method1==0)
{
sampledname1=paste0(hadprocesspath,"/",ofisuffix1,".nc")
writeout(inras1,sampledname1, var1, varunit1, var1)
sampledname2=paste0(hadprocesspath,"/",ofisuffix1,"_180.nc")
writeout(inras2,sampledname2, var1, varunit1, var1)
}
if(method1==1)
{
## orig 67
width1 = 67
##width1=157
outras1<- focal(inras1, w = matrix(1,width1,width1), fun = fill.na, pad = TRUE, na.rm = FALSE)
outras2<- focal(inras2, w = matrix(1,width1,width1), fun = fill.na, pad = TRUE, na.rm = FALSE)
sampledname1=paste0(hadprocesspath,"/",ofisuffix1,".nc")
writeout(outras1,sampledname1, var1,varunit1, var1)
sampledname2=paste0(hadprocesspath,"/",ofisuffix1,"_180.nc")
writeout(outras2,sampledname2, var1,varunit1, var1)
}
if(method1==2)
{
outras1<-hadcm3_raster_fillna_method_2(inras1)
outras2<-hadcm3_raster_fillna_method_2(inras2)
sampledname1=paste0(hadprocesspath,"/",ofisuffix1,".nc")
writeout(outras1,sampledname1, var1,varunit1, var1)
sampledname2=paste0(hadprocesspath,"/",ofisuffix1,"_180.nc")
writeout(outras2,sampledname2, var1,varunit1, var1)
}
if(method1==4)
{
sampler1 <- raster(ncol=samplecols1, nrow=samplerows1)
outras01<-resample(inras1, sampler1)
outras02<-resample(inras2, sampler1)
width1=7
outras1<- focal(outras01, w = matrix(1,width1,width1), fun = fill.na7, pad = TRUE, na.rm = FALSE)
outras2<- focal(outras02, w = matrix(1,width1,width1), fun = fill.na7, pad = TRUE, na.rm = FALSE)
plot(outras1)
sampledname1=paste0(hadprocesspath,"/",ofisuffix1,".nc")
writeout(outras1,sampledname1, var1,varunit1, var1)
sampledname2=paste0(hadprocesspath,"/",ofisuffix1,"_180.nc")
writeout(outras2,sampledname2, var1,varunit1, var1)
}
}
test_code_raster_experiment_1<-function()
{
beginyear1=45000-5
years1=10
month1=7
varr1="tas"
varr2="pr"
method1=4
ovarnam1="TS"
varunit1="C"
ofisuffix1="tempin"
#get_hadgcm3_global_raster(beginyear1, years1, varr1, month1, method1)
get_hadgcm3_global_raster(beginyear1, years1, varr1, month1,varunit1, ovarnam1,ofisuffix1,method1)
}
- Main proggis
-
- main running ...
print("TS 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)
{
"Orography ..."
preprocess_orography()
}
if(get_srtm_data==1)
{
get_dordogne_srtm()
}
if(get_gtopo30_data==1)
{
get_europe_gtopo30()
}
if(capture_hadcm3_60ka_data==1)
{
if(capture_temperature==1)
{
capture_hadcm3_temperature=1
}
if(capture_rainfall==1)
{
capture_hadcm3_rainfall=1
}
if(capture_hadcm3_temperature==1)
{
print("Processing hadcm3 temperature data to snapshot ...")
hadvariaabeli="tas"
#hadfilename<<-generate_hadfilename(hadbasepath, hadyear, hadvariaabeli)
#print (hadyear)
#load_average_hadcm3_years_months(hadfilename, hadvariaabeli,hadbaseyear, hadyear,hadmonth, hadyears, hadmonths, hadoperaatio,hadcm3_60ka_seafill_method)
#byr1=hadbaseyear-(hadyears/2)
byr1=fage
years1=numyears
month1=fmonth
varr1=hadvariaabeli
method1=hadcm3_60ka_seafill_method
ovarnam1="TS"
varunit1="C"
ofisuffix1="tempin"
#print ("BYR")
#print (byr1)
get_hadgcm3_global_raster(byr1, years1, varr1, month1,varunit1, ovarnam1,ofisuffix1,method1)
}
if(capture_hadcm3_rainfall==1)
{
print("Processing hadcm3 rainfall data to snapshot ...")
hadvariaabeli="pr"
#hadfilename<<-generate_hadfilename(hadbasepath, hadyear, hadvariaabeli)
#load_average_hadcm3_years_months(hadfilename, hadvariaabeli,hadbaseyear, hadyear,1, hadyears, 12, 3, hadcm3_60ka_seafill_method)
byr1=fage
years1=numyears
month1=fmonth
varr1=hadvariaabeli
method1=hadcm3_60ka_seafill_method
ovarnam1="PR"
varunit1="C"
ofisuffix1="precipin"
#print ("BYR")
#print (byr1)
get_hadgcm3_global_raster(byr1, years1, varr1, month1,varunit1, ovarnam1,ofisuffix1,method1)
}
}
if(capture_trace26k_data==1)
{
if(capture_temperature==1)
{
print("Processing input temperature variable to snapshot ...")
tracelocation3(predata2, variaabeli, fage, numyears, fmonth);
}
if(capture_rainfall==1)
{
print("Processing input rainfall variable to snapshot ...")
tracelocation3(predata2, variaabeli, 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()
}
if(make_rain_shadows==1)
{
if(capture_elevation==2)
{
# force method to 0 because of peltier topography
downscale_orography(0)
}
else
{
downscale_orography(method_ds_oro)
}
## generate rain shadows, assume westerly wind by default
print("Generating rain shadows.")
generate_rain_shadows()
}
print ("Downscaling ...")
- do ts downscaling
if(normal_ds==1)
{
print("Normal etopo1 area downscaling.")
## default area DS
normal_ts_downscale()
}
if(normal_ds==3)
{
print("Rainfall downscaling.")
#out_text_1="Rainfall"
#out_text_2="mm/yr"
#out_text_3="IceAgeRain"
## default area DS
downscale_rainfall()
}
if(global_ds==1)
{
print("World downscaling.")
world_ts_downscale()
}
if(accurate_ds==1)
{
print ("Customized dwonscaling.")
custom_ts_downscale_1(kustomorofilee1,outvarfname1 )
}
if(plot_result==1)
{
print("Plotting of draft.")
plottaus1(outvarfname1,plotfname1 )
}
- try(system("python maplot1.py", intern = TRUE, ignore.stderr = TRUE))
print (fage)
print (fmonth)
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.
- ↑ TY - JOUR AU - Armstrong, Edward AU - Hopcroft, Peter O. AU - Valdes, Paul J. PY - 2019 DA - 2019/11/07 TI - A simulated Northern Hemisphere terrestrial climate dataset for the past 60,000 years JO - Scientific Data SP - 265 VL - 6 IS - 1 AB - We present a continuous land-based climate reconstruction dataset extending back 60 kyr from 0 BP (1950) at 0.5° resolution on a monthly timestep for 0°N to 90°N. It has been generated from 42 discrete snapshot simulations using the HadCM3B-M2.1 coupled general circulation model. We incorporate Dansgaard-Oeschger (DO) and Heinrich events to represent millennial scale variability, based on a temperature reconstruction from Greenland ice-cores, with a spatial fingerprint based on a freshwater hosing simulation with HadCM3B-M2.1. Interannual variability is also added and derived from the initial snapshot simulations. Model output has been downscaled to 0.5° resolution (using simple bilinear interpolation) and bias corrected. Here we present surface air temperature, precipitation, incoming shortwave energy, minimum monthly temperature, snow depth, wind chill and number of rainy days per month. This is one of the first open access climate datasets of this kind and can be used to study the impact of millennial to orbital-scale climate change on terrestrial greenhouse gas cycling, northern extra-tropical vegetation, and megaflora and megafauna population dynamics. SN - 2052-4463 UR - https://doi.org/10.1038/s41597-019-0277-1 DO - 10.1038/s41597-019-0277-1 ID - Armstrong2019 ER -
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 05:31, 9 September 2020 | 1,472 × 976 (1.01 MB) | Merikanto (talk | contribs) | Update of script | |
12:45, 7 September 2020 | 2,000 × 1,545 (2.18 MB) | Merikanto (talk | contribs) | Uploaded own work with UploadWizard |
You cannot overwrite this file.
File usage on Commons
There are no pages that use this file.