File:Climate diagram of Paris 40750 bp 1.svg

1. 2. acquire some hadcm3b 60ka climate data
   3. "R" 4.03
   4. v. 0002.02
   5. 17.10.2021
   6. WARNING: script in alpha stage

1. install_libs1=1

1. if(install_libs==1)
2. {
3. install.packages("raster")
4. install.packages("ncdf4")
5. install.packages("abind")
6. install.packages("Cairo")
7. install.packages("svglite")
8. }

library(raster) library(ncdf4) library(abind) library(svglite)

## hadcm3b 60ka files path

1. 1. hadbasepath<<-"D:/varasto_iceagesimu"

hadbasepath<<-"./predata"

hadbaseyear=-1 hadprocesspath<-"./data_processing/"

lones1=0 latis1=0

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

1. print (hadnames1[1])
2. print (hadnames1[2])
3. print(deltayears1)
4. print(deltamonths1)
5. print(deltayears2)
6. print(deltamonths2)

twoo1=0

if(hadbaseyears[1]==hadbaseyears[2]) {

1. 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)

1. 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) {

1. 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(beginyr1, yrs1, varr1) { yrmid1=beginyr1+(yrs1/2)

slaici1<-had_twoslicer(beginyr1,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")

}

load_climate<-function() {

1. beginyear1=36000
   1. beginyear1=40200

beginyear1=40750

years1=100

month1=7

varr1="tas" varr2="pr"

1. 1. paris
2. beginyear1=40750

targetname1="paris" targetlat1=48.856667 targetlon1=2.351111

1. 1. selerika 64.66666,147.833333
   2. selerika 64° 40' N, 147° 45' E
2. targetname1="selerika"
3. targetlat1=64.666667
4. targetlon1=147.833333

1. targetname1="zyryanka"
2. targetlat1=65.75
3. targetlon1=150.9

1. targetname1="seymchan"
2. targetlat1=62.930833
3. targetlon1=152.385

1. targetname1="sungir"
2. targetlat1=56.175833
3. targetlon1=40.509167

1. 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)

placename1=targetname1 yearr1=as.character(beginyear1)

sj1=paste0("python hadiag1.py ",placename1," ",yearr1)

print(sj1)

system(sj1)

}

raster_experiment_1<-function() { beginyear1=40650 years1=100 month1=7 varr1="tas" varr2="pr"

load_had_rasters_var(beginyear1, years1, varr1) load_had_rasters_var(beginyear1, years1, varr2) }

load_python_draw_climate<-function(beginyear1, targetname1, targetlon1, targetlat1) {

years1=33 ## num of yrs to average

1. month1=7

varr1="tas" varr2="pr"

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)

placename1=targetname1 yearr1=as.character(beginyear1)

sj1=paste0("python hadiag1.py ",placename1," ",yearr1)

print(sj1)

system(sj1)

}

1. 2. Main proggis

print("HadCM3B 60ka simulation climate data.")

beginyear1=40750 targetname1="paris" targetlat1=48.856667 targetlon1=2.351111

load_python_draw_climate(beginyear1, targetname1, targetlon1, targetlat1)

print("Program run done.")

1. 1. raster_experiment_1()

1. load_climate()

1. 3. drawing climate diagram in python 3
   4. from input csv file
   5. version 2.1101
   6. 17.10.2021

import matplotlib.pyplot as plt import numpy as np import pandas as pd from scipy import interpolate import sys

print ('Argument List:', str(sys.argv))

pohjanimi=sys.argv[1] ika=sys.argv[2] isonimi=pohjanimi.capitalize()

print(pohjanimi, isonimi, ika)

1. quit(-1)

1. pohjanimi="paris"
2. ika="40750"

captioni=isonimi+", "+ika+" BP" maxrainfall=120 mintemperature=-40 maxtemperature=20

datafilename=pohjanimi+".csv" savename=pohjanimi+"_"+ika+"_climate_diagram.svg"

figsizex=12 figsizey=8

x0 = [] y0 = [] y20= []

x = [] y = [] y2= []

dfin0=pd.read_csv(datafilename, sep=";") lst1 = ['Month','T','P']

dfin1 = dfin0[dfin0.columns.intersection(lst1)]

x0=dfin1['Month'] y0=dfin1['T'] y20=dfin1['P']

x.append(0) y.append(y0[11]) y2.append(y0[11])

for n in range(0, 12): x.append(x0[n]) y.append(y0[n]) y2.append(y20[n])

x.append(13) y.append(y0[0]) y2.append(y0[0])

print(x)

1. print(y)
2. print (type(x))
3. print (type(y))

1. quit(0)

yearprecip=0 yeartemp=0

for n in range(1, 13): yearprecip=yearprecip+y2[n] yeartemp=yeartemp+y[n] print (n,y[n],y2[n])

size1=22 size2=26 size3=30

yeartemp=round((yeartemp/12.0),1) mintemp=min(y) maxtemp=max(y) yearprecip=round(yearprecip,0) maxprecip=max(y2) minprecip=min(y2)

print(yearprecip) print(minprecip) print(maxprecip)

print(yeartemp) print(mintemp) print(maxtemp)

ymax1=int((maxprecip+60)/20)*20 ymax2=int((maxtemp+15)/5)*5 ymin2=int((mintemp-10)/5)*5

x_sm = np.array(x) y_sm = np.array(y) x_smooth = np.linspace(x_sm.min(), x_sm.max(), 200) funk1 = interpolate.interp1d(x_sm, y_sm, kind="quadratic") y_smooth = funk1(x_smooth)

fig, ax1 = plt.subplots()

1. plt.rcParams["figure.figsize"] = (12,16)

ax1.axis((1,12,0,ymax1))

ax1.bar(x, y2, color='#0000ff', label="Precip. mm", width=0.9, align="center")

ax1.set_ylabel('Precipitation mm', color='#00007f', fontsize=size2)

for tl in ax1.get_yticklabels():

tl.set_color('b')
tl.set_fontsize(size1)

ax2 = ax1.twinx() ax2.set_ylabel('Temperature °C', color='#7f0000', fontsize=size2)

ax2.axis((1,12,ymin2, ymax2))

1. ax2.plot(x,y, label='Temperature °C',color="#ff0000", linewidth=7)

ax2.plot(x_smooth,y_smooth, label='Temperature °C',color="red", linewidth=10)

for t2 in ax2.get_yticklabels():

t2.set_color('r')
t2.set_fontsize(size1)

ax1.set_xlabel('Month', color="darkgreen", fontsize=size2)

for tix in ax1.get_xticklabels():

tix.set_color("Black")
tix.set_fontsize(size1)

ax1.set_title(captioni, fontsize=size3)

ax2.text(1, ymax2-4, " P annual "+str(int(yearprecip))+ " mm", color="#00007f", fontsize=size1) ax2.text(1, ymax2-8, " T year "+str(yeartemp) + " °C", color="#7f0000",fontsize=size1) ax2.text(1, ymax2-12, " T max "+str(maxtemp)+ " °C", color="#7f0000", fontsize=size1) ax2.text(1, ymax2-16, " T min "+str(mintemp) + " °C", color="#7f0000",fontsize=size1)

fig = plt.gcf() fig.set_size_inches(figsizex, figsizey, forward=True)

plt.plot()

plt.savefig(savename, format="svg", dpi = 100)

plt.show()