File:Fairbanks 40750 climate diagram.svg
Original file (SVG file, nominally 1,080 × 720 pixels, file size: 47 KB)
Captions
Summary
[edit]DescriptionFairbanks 40750 climate diagram.svg |
English: Climate diagram of Fairbanks, 40750 years ago |
Date | |
Source | Own work |
Author | Merikanto |
Camera location | 64° 50′ 00″ N, 147° 43′ 00″ W | View this and other nearby images on: OpenStreetMap | 64.833333; -147.716667 |
---|
- acquire some hadcm3b 60ka climate data
- "R" 4.03
- v. 0002.02
- 17.10.2021
- WARNING: script in alpha stage
-
- install_libs1=1
- if(install_libs==1)
- {
- install.packages("raster")
- install.packages("ncdf4")
- install.packages("abind")
- install.packages("Cairo")
- install.packages("svglite")
- }
library(raster)
library(ncdf4)
library(abind)
library(svglite)
## hadcm3b 60ka files path
- 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
- 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(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()
{
- beginyear1=36000
- beginyear1=40200
beginyear1=40750
years1=100
month1=7
varr1="tas"
varr2="pr"
- paris
- beginyear1=40750
targetname1="paris"
targetlat1=48.856667
targetlon1=2.351111
- selerika 64.66666,147.833333
- selerika 64° 40' N, 147° 45' E
- targetname1="selerika"
- targetlat1=64.666667
- targetlon1=147.833333
- targetname1="zyryanka"
- targetlat1=65.75
- targetlon1=150.9
- targetname1="seymchan"
- targetlat1=62.930833
- targetlon1=152.385
- 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)
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
- 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)
}
- Main proggis
print("HadCM3B 60ka simulation climate data.")
beginyear1=40750
targetname1="paris"
targetlat1=48.856667
targetlon1=2.351111
beginyear1=40750
targetname1="selerika"
targetlat1=64.66666
targetlon1=147.833333
beginyear1=40750
targetname1="fairbanks"
targetlat1=64.833333
- targetlon1=-147.716667
targetlon1=212.283333
load_python_draw_climate(beginyear1, targetname1, targetlon1, targetlat1)
print("Program run done.")
- raster_experiment_1()
- load_climate()
- drawing climate diagram in python 3
- from input csv file
- version 2.1101
- 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)
- quit(-1)
- pohjanimi="paris"
- 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)
- print(y)
- print (type(x))
- print (type(y))
- 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()
- 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))
- 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()
Licensing
[edit]- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
- share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 09:23, 17 October 2021 | 1,080 × 720 (47 KB) | Merikanto (talk | contribs) | Uploaded own work with UploadWizard |
You cannot overwrite this file.
File usage on Commons
The following page uses this file:
Metadata
This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Width | 864pt |
---|---|
Height | 576pt |