File:Fairbanks 20000 bp climate diagram.svg

1. 2. chelsa trace21ka climate data download
2. and climate diagram w/external python script
3. 1. extract monthly data
4. 1. in this version calibration of alaska inland temperatures
   2. with zagoskin lake chironomid data
   5. 20.10.2021 0000.0002

library(raster) library(sp) library("readxl") library("tidyverse")

download_chelsa_trace21ka_file<-function(linkdir1,variable1,month1, kyrs1) { linkbase1="https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/" filebase1="CHELSA_TraCE21k_" suffix1="_V1.0.tif"

filename1=paste0(filebase1,variable1,"_",as.character(month1),"_-",as.character(kyrs1),suffix1) print(filename1)

inlink1=paste0(linkbase1,linkdir1,filename1) infn1=paste0(filename1)

print(inlink1) print(infn1)

download.file(inlink1,infn1,method="curl", mode="wb") }

download_chelsa_data_files<-function(kyrs1) { indir1="./indata" ## dir.remove(indir1) dir.create(indir1) setwd(indir1)

for (n in 1:12) { download_chelsa_trace21ka_file("pr/","pr",n, kyrs1) }

for (n in 1:12) { download_chelsa_trace21ka_file("tasmax/","tasmax",n, kyrs1) }

for (n in 1:12) { download_chelsa_trace21ka_file("tasmin/","tasmin",n, kyrs1) }

}

load_climate_data<-function(sitename,age, lon1, lat1) { outfilename=paste0(tolower(sitename),".csv")

agecents1=age/100

precips=c() temps1=c() temps2=c() temps=c()

for (n in 1:12) { variable1="pr" precip1<-get_value(variable1, lon1, lat1,agecents1, n) precips=c(precips, precip1) }

for (n in 1:12) { variable1="tasmax" temp1<-get_value(variable1, lon1, lat1,agecents1, n) temps1=c(temps1, temp1) }

for (n in 1:12) { variable1="tasmin" temp2<-get_value(variable1, lon1, lat1,agecents1, n) temps2=c(temps2, temp2) }

temps0=(temps1+temps2)/2 temps1=(temps0/10)-273.15 temps=round(temps1, digits = 1)

months=1:12

print(precips) print(temps)

daffy1<-data.frame(months, temps, precips)

names(daffy1)<-c("Month","T","P") print(daffy1)

write.table(daffy1, file = outfilename, sep = ";",row.names = TRUE, col.names = NA)

return(temps[7])

}

create_climate_plot<-function(sitename,age, lon1, lat1) {

load_climate_data(sitename,age, lon1, lat1)

placename1=sitename yearr1=as.character(age)

1. 1. we use external python script to plot

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

print(sj1)

system(sj1)

}

read_chironomid_datafile<-function(zagoname1,aake1) { cin1=read.csv(zagoname1, sep=";") #print(cin1) ager1=as.vector(cin1['Age']) temper1=as.vector(cin1['July_T']) lusmu1=which(ager1>aake1) a2=lusmu1[1] a1=a2-1 ake1=cin1[a1,1] t1=cin1[a1,2] ake2=cin1[a2,1] t2=cin1[a2,2] agestep=ake2-ake1 tstep=t2-t1 koo=tstep/agestep dake=aake1-ake1 tempera=t1+dake*koo return(tempera) }

read_climate_datafile<-function(filename1) { cin1=read.csv(filename1, sep=";") return(cin1) }

write_climate_datafile<-function(daffy1, filename1) { #cin1=read.csv(filename1, sep=";") if (file.exists(filename1)) { file.remove(filename1) } write.table(daffy1, filename1, sep = ";",row.names = TRUE, col.names = NA) return(daffy1) }

filter_temperature_data<-function (daffy1, tofset1) { print(" TO FILT") print(daffy1)

df1 = data.frame(Month = numeric(), T = numeric(),P=numeric() )

   names(df1)<-c("Month","T","P")

print (daffy1$July_T)

print(tofset1)

for (n in 1:12) {

kuku1<-daffy1[n,2]

## JN WARNING

tofu1<-daffy1[n,3]+tofset1 #tofu1<-(daffy1[n,8]*1.0)+tofset1

tofu2<-round(tofu1,2) pepe1<-daffy1[n,4]

df1[nrow(df1) + 1,] = c(kuku1,tofu2,pepe1)

}

print("FILTETED")

print(df1)

return(df1)

}

filter_climate_datafile<-function(filename1, tofset1) { print("Filtering climate data by TJuly offset ...") print(tofset1) daffy1<-read_climate_datafile(filename1)

daffy2<-filter_temperature_data(daffy1, tofset1) print("Filtered data:") print(daffy2)

write_climate_datafile(daffy2, filename1) print(filename1)

return(daffy2)

}

load_zagoskin_lake_midge_data_from_net_xls<-function() { #targeturl1="https://www.ncei.noaa.gov/pub/data/paleo/insecta/chironomidae/northamerica/usa/alaska/zagoskin-burial2009.txt" targeturl1="https://www.ncei.noaa.gov/pub/data/paleo/insecta/chironomidae/northamerica/usa/alaska/zagoskin-burial2009.xls" targetfile1="zagoskin_input.xls"

download.file(targeturl1, targetfile1, method="curl")

excel1=read_excel(targetfile1)

path1 <- targetfile1

excel1<-path1 %>% excel_sheets() %>% set_names() %>% map(read_excel, path = path1)

#excel1

#str(excel1)

zagoskin01<-excel1[3]

zagoskin02<-zagoskin01$`Zagoskin Lake`

zages01<-zagoskin02[2] ztemps01<-zagoskin02[3]

#zages01 #ztemps01

zago_data_1<-data.frame(zages01,ztemps01)

names(zago_data_1)<-c("Age", "July_T")

write.table(zago_data_1,"zagoskinlake.csv",sep=";")

#print(zago_data_1)

return(zago_data_1) }

get_value<-function(variable1, lon1, lat1,agecents1, month1) { filebase1="CHELSA_TraCE21k_" suffix1="_V1.0.tif"

filename1=paste0(filebase1,variable1,"_",as.character(month1),"_-",as.character(agecents1),suffix1)

inf1=paste0("./indata/", filename1) rin1<-raster(inf1)

#print(filename1) #print(inf1)

xy1 <- SpatialPoints(coords = cbind(lon1,lat1))

#print(rin1) #print(xy1)

## must use raster::, because tidyr has also extract()! value1<-raster::extract(rin1, xy1)

 return(value1)

}

1. 1. main proge
2. 1. set 1, if you want download data, must do tht in first time
3. USAGE is some complicated
   1. by first, to download date, set download_?=1 , plot_data=0
4. drawer does not func
5. download_data=1
6. download_alaska_midge_data=1
7. plot_data=0
   1. by second, to made plot, set plot_data=0, drawer is ok
8. download_data=0
9. download_alaska_midge_data=0
10. plot_data=1

download_data=0 ## obtain chelsa trace 21 kadata ad 2021.10 download_alaska_midge_data=0

1. 1. secont, set plot 1

plot_data=1

filter_by_zagoskin=1 ## filter plotted alaska temperatures w/zagoskin lake chironomid data

1. age=20000
   1. age=14000

age=20000

sitename="Fairbanks" lon1=-147.723056 lat1=64.843611

1. sitename="Selerika"
2. lat1=64.66666
3. lon1=147.833333

1. sitename="Kitluk"
2. lat1=66.50939782345924
3. lon1=-164.42294254249194

1. sitename="Paris"
2. lat1=48.856613
3. lon1=2.352222

1. sitename="Yukagir"
2. lat1=71.75
3. lon1=139.883333

1. 1. sitename="Kiev"
   2. lat1=50.45
   3. lon1=50.45

1. sitename="Les_Eyzies"
2. lat1=44.936667
3. lon1=1.018056

1. sitename="Rouffignac"
2. lat1=45.008889
3. lon1=0.9875

1. 1. must download chelsa trace21ka files, if you have not

if(download_alaska_midge_data==1) { print("Download Alaska midge paleo data ...") load_zagoskin_lake_midge_data_from_net_xls() }

if(download_data==1) { print("Download Chelsa Trace21ka data ...")

agecents1=age/100 download_chelsa_data_files(agecents1) }

if(plot_data==1) {

1. Sys.sleep(20)

if(filter_by_zagoskin==1) {

1. 1. filter alaska data by zakoskun laje Age, T_July chironomid
   2. to get better july temp

print("Filtering climate data by offset ...") print("Plotting, correction with Zagoskin lake data ....")

zagotemp1=read_chironomid_datafile("./zagoskinlake.csv",age)

print(zagotemp1)

zagosimu1=load_climate_data("zagoskin",age, -162.10752764569278, 63.44888166191213)

print(zagosimu1)

filename1=paste0(tolower(sitename),".csv") ## loaf fairbaks etc load_climate_data(sitename,age, lon1, lat1)

tofset1<-zagotemp1-zagosimu1 #tofset1<-round(tofset1,1)

print (tofset1)

print(sitename)

filter_climate_datafile(filename1, tofset1)

placename1=sitename yearr1=as.character(age)

## we use external python script to plot

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

print(sj1)

system(sj1)

}

if(filter_by_zagoskin==0) { print("Plotting, no correction w/Zagoskin lake data ....") create_climate_plot(sitename,age, lon1, lat1) }

} # ... plot data

1. 3. drawing climate diagram in python 3
   4. from input csv file
   5. version 2.1102
   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+"_bp_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()