Commit 35a29e3b authored by saeed's avatar saeed
Browse files

m

parent a8076a34
......@@ -7,30 +7,77 @@ static int isleapyear(int year){
static int getlastday(int year, int month){
int result = 0;
static int lastday_noleap[] = {31,28,31,30,31,30,31,31,30,31,30,31};
static int lastday_leap [] = {31,29,31,30,31,30,31,31,30,31,30,31};
if (isleapyear(year) ) return lastday_leap [month-1];
if (!isleapyear(year)) return lastday_noleap[month-1];
return 0;
if (isleapyear(year) ) result = lastday_leap [month-1];
if (!isleapyear(year)) result = lastday_noleap[month-1];
return result;
}
static void get_yearmonday(char *date, int *y, int *m, int *d){
char *st = date;
if (date){
char tmp[5] = {'\0'};
for (int i=0; i < 4;++i){
tmp[i] = *st;
st++;
}
tmp[4] = '\0';
*y = atoi(tmp);
char tmp1[3] = {'\0'};
for (int i=0; i < 2;++i){
tmp1[i] = *st;
st++;
}
tmp1[2] = '\0';
*m = atoi(tmp1);
char tmp2[3] = {'\0'};
int main(void){
for (int i=0; i < 2;++i){
tmp2[i] = *st;
st++;
}
tmp2[2] = '\0';
*d = atoi(tmp2);
}
}
FILE *file = NULL;
file = fopen("start_end.dat","w+");
int y1 = 2015;
int m1 = 2;
int d1 = 16;
int y2 = 2016;
int m2 = 10;
int d2 = 13;
int yt = y1, mt = m1;
int lastday = 0;
int nfirst = 0;
if (y1 < y2){
while (yt < y2){
lastday= getlastday(yt, mt);
int main(int argc, char *argv[]){
if (argc == 3){
char *startdate = argv[1];
char *enddate = argv[2];
int y1 = 0;
int m1 = 0;
int d1 = 0;
int y2 = 0;
int m2 = 0;
int d2 = 0;
get_yearmonday(startdate,&y1,&m1,&d1);
get_yearmonday(enddate,&y2,&m2,&d2);
FILE *file = NULL;
file = fopen("start_end.dat","w+");
int yt = y1, mt = m1;
int lastday = 0;
int nfirst = 0;
if (y1 < y2){
while (yt < y2){
lastday= getlastday(yt, mt);
if (mt == m1 && nfirst == 0){
fprintf(file,"%04d%02d%02d\t%04d%02d%02d\n",yt, mt,d1,yt, mt,lastday);
nfirst = 1;
......@@ -57,7 +104,7 @@ int main(void){
mt = m2;
lastday = d2;
fprintf(file,"%04d%02d%02d\t%04d%02d%02d\n",yt, mt,1,yt, mt,lastday);
}
}
if (y1 == y2){
if (m1 < m2){
......@@ -78,7 +125,7 @@ int main(void){
}
fclose(file);
}
return 0;
}
import datetime
import sys
import numpy as np
import math
import matplotlib.pyplot as plt
from matplotlib.pyplot import *
#-------------------------------------------------------------------------------
def strDate(day):
cyear, cmonth, cday="%04d"%int(day.year),"%02d"%int(day.month),"%02d"%int(day.day)
return "".join([cyear, cmonth, cday])
def getDateRange(fromt, tot):
initd = datetime.datetime(int(fromt[0:4]),int(fromt[4:6]),int(fromt[6:8]),int(fromt[8:10]))
totd = datetime.datetime(int(tot[0:4]),int(tot[4:6]),int(tot[6:8]),int(tot[8:10]))
step = datetime.timedelta(hours=1)
dateticks = []
while initd <= totd:
dateticks.append(strDate(initd))
initd = initd + step
return dateticks
def filldata(date_obs,fromt, tot):
initd = datetime.datetime(int(fromt[0:4]),int(fromt[4:6]),int(fromt[6:8]))
totd = datetime.datetime(int(tot[0:4]),int(tot[4:6]),int(tot[6:8]))
step = datetime.timedelta(hours=24)
dmcor = {}
while initd <= totd:
if strDate(initd) not in date_obs:
dmcor[strDate(initd)] = {}
dmcor[strDate(initd)] = -9999.0
initd = initd + step
return dmcor
def getera5obs(cera5,cbeg,cend):
fera5 = open(cera5,"r")
vobs_era5= []
vtime_era5 = []
for vl in fera5:
vobs_era5.append(float(vl.split()[1]))
vtime_era5.append(vl.split()[0])
fera5.close()
### fill era5
dera5=dict(zip(vtime_era5,vobs_era5))
dmiss = filldata(vtime_era5,cbeg,cend)
if dmiss.keys():
dera5.update(dmiss)
era_fill = []
ktime=sorted(dera5.keys())
for nt in range(len(ktime)):
era_fill.append(dera5[ktime[nt]])
eram = np.ma.masked_where(np.array(era_fill)==-9999.0,np.array(era_fill))
return eram, ktime
def getprecobs(cprec,cbeg,cend):
fprec = open(cprec,"r")
vobs_prec= []
vtime_prec = []
for vl in fprec:
vobs_prec.append(float(vl.split()[1]))
vtime_prec.append(vl.split()[0])
fprec.close()
dprec=dict(zip(vtime_prec,vobs_prec))
dmiss = filldata(vtime_prec,cbeg,cend)
if dmiss.keys():
dprec.update(dmiss)
prec_fill = []
ktime=sorted(dprec.keys())
for nt in range(len(ktime)):
prec_fill.append(dprec[ktime[nt]])
#prec_fill
precm = np.ma.masked_where(np.array(prec_fill)==-9999.0,np.array(prec_fill))
return precm,ktime
def getreportype(cera5):
st1 = cera5.replace("_thin.txt","")
st2 = st1.replace("odb_era5_","")
report_type = st2.replace("_cycle12","")
return report_type
def getcycle(cera5):
st1 = cera5.replace("_thin.txt","")
il = st1.find("_cycle")
return st1[il+len("_cycle"):]
def getreportype_era5(cera5):
il = cera5.find("odb_era5_")
st1=""
result=""
if il >= 0: st1=cera5[il+len("odb_era5_"):]
i2=st1.find("_")
if i2 >=0:
result=st1[0:i2]
return result
def getreportype_prec(cprec):
il = cprec.find("odb_precise_")
st1=""
result=""
if il >= 0: st1=cprec[il+len("odb_precise_"):]
i2=st1.find("_")
if i2 >=0:
result=st1[0:i2]
return result
def check_same_reporttype_cycle(cera5,cprec):
result = 0
rep_era5 = getreportype_era5(cera5)
rep_prec = getreportype_prec(cprec)
cycle_era5 = getcycle(cera5)
cycle_prec = getcycle(cprec)
print (rep_era5, rep_prec, cycle_era5, cycle_prec)
if rep_era5 != rep_prec or cycle_era5 != cycle_prec:
result=1
return result
def prechassamereporttype(cprec,report_type):
ind = cprec.find(report_type)
if ind == -1 :
result = 1
else :
result = 0
return result
def report_type_info(report_type):
fcsv = open("reportype.csv","r")
drtype = {}
for vl in fcsv:
vinfo = []
row = vl.split(";")
vinfo.append(row[1])
vinfo.append(row[2])
drtype[row[0]] = vinfo
fcsv.close()
return drtype[report_type][0]
def plot_comparsion(ktime,eram,precm,report_type,cycle,sreport_type,cbeg,cend):
print (ktime)
print (eram)
print (precm)
fig=figure(num=1,figsize=(12.0,6.0),dpi=300,facecolor='w',edgecolor='k')
ax = axes([0.08, 0.08, 0.8, 0.8], facecolor= '1.0')
ax.set_xlabel('Time',fontsize='15',weight='bold')
ax.set_ylabel ('number of observations',fontsize='15',weight='bold')
ax.grid()
jtar = [i for i in range(len(ktime))]
ax.xaxis.set_ticks(jtar)
ax.plot(precm,'-',color='r')
ax.plot(eram,'-',color='b')
ax.xaxis.set_ticklabels(ktime,rotation="vertical")
ax.legend(["PRECISE","ERA5"], numpoints=1, prop=dict(size='small'),loc="best")
ax.set_title("report type " + report_type + "(" +sreport_type+") "+"cycle"+cycle,fontsize="15")
figout = "numobs_era5_prec_reporttype"+ str(report_type) + "_cycle"+cycle+"_"+cbeg+"_"+cend+".png"
print ("creating file ", figout)
savefig(figout, bbox_inches='tight',dpi=300,facecolor='w',edgecolor='w',orientation='portrait')
close(1)
if __name__ == '__main__':
if len(sys.argv) != 5:
print ("Usage "+sys.argv[0]+" era5_data precise_data startdate(YYYYMMDD) endate(YYYYMMDD)")
sys.exit()
cera5 = sys.argv[1]
cprec = sys.argv[2]
cbeg = sys.argv[3]
cend = sys.argv[4]
err = check_same_reporttype_cycle(cera5,cprec)
if err == 1:
print (" era5 aand precise does not have the same report type or cycle so we exit")
sys.exit()
eram, ktime = getera5obs(cera5,cbeg,cend)
precm, ktime = getprecobs(cprec,cbeg,cend)
report_type = getreportype_era5(cera5)
cycle = getcycle(cera5)
print (report_type, cycle)
string_report_type = report_type_info(str(report_type))
plot_comparsion(ktime,eram,precm,report_type,cycle,string_report_type,cbeg,cend)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment