Skip to content

Commit

Permalink
fix time axis issue, add new option to specify input file, Ed input o…
Browse files Browse the repository at this point in the history 
…n plot
  • Loading branch information
pmathiot committed Nov 27, 2019
1 parent bfe38cc commit eec9991
Show file tree
Hide file tree
Showing 5 changed files with 121 additions and 43 deletions.
2 changes: 1 addition & 1 deletion SCRIPT/mk_mht.bash
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,4 +39,4 @@ else
fi

# extract only 26.5
ncks -d y,793,793 $FILEOUT nemo_${RUN_NAME}o_${FREQ}_${TAG}_mht_265.nc
ncks -O -d y,793,793 $FILEOUT nemo_${RUN_NAME}o_${FREQ}_${TAG}_mht_265.nc
112 changes: 76 additions & 36 deletions SCRIPT/plot_time_series.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,11 @@
import matplotlib.colors as colors
import sys
import re
import datetime as dt
import argparse
import pandas as pd
import matplotlib.dates as mdates
import matplotlib.ticker as ticker

# def class runid
class run(object):
Expand All @@ -28,45 +30,49 @@ def load_time_series(self, cfile, cvar):
for kf,cf in enumerate(cfile):
try:
ncid = nc.Dataset(cf)
except Exception as e:
print 'issue in trying to load file : '+cf
print e
sys.exit(42)

ncvtime = ncid.variables[ctime]
if 'units' in ncvtime.ncattrs():
cunits = ncvtime.units
else:
cunits = "seconds since 1900-01-01 00:00:00"

# define calendar
if 'calendar' in ncvtime.ncattrs():
ccalendar = ncvtime.calendar
else:
ccalendar = "noleap"
time = nc.num2date(ncid.variables[ctime][:].squeeze(), cunits, ccalendar)
ncvtime = ncid.variables[ctime]
if 'units' in ncvtime.ncattrs():
cunits = ncvtime.units
else:
cunits = "seconds since 1900-01-01 00:00:00"

# convert to proper datetime object
if isinstance(time,list):
ntime = time.shape[0]
else:
ntime = 1

timeidx=[None]*ntime
for itime in range(0, ntime):
# define calendar
if 'calendar' in ncvtime.ncattrs():
ccalendar = ncvtime.calendar
else:
ccalendar = "noleap"
time = nc.num2date(ncid.variables[ctime][:].squeeze(), cunits, ccalendar)

# convert to proper datetime object
if isinstance(time,list):
timeidx[itime] = np.datetime64(time[itime],'D')
ntime = time.shape[0]
else:
timeidx[itime] = np.datetime64(time,'D')
ntime = 1

# build series
cnam=get_varname(cf,cvar)
df[kf] = pd.Series(ncid.variables[cnam][:].squeeze()*self.sf, index = timeidx, name = self.name)
timeidx=[None]*ntime
for itime in range(0, ntime):
if isinstance(time,list):
timeidx[itime] = np.datetime64(time[itime],'D')
else:
timeidx[itime] = np.datetime64(time,'D')

# build series
cnam=get_varname(cf,cvar)
df[kf] = pd.Series(ncid.variables[cnam][:].squeeze()*self.sf, index = timeidx, name = self.name)

except Exception as e:
print 'issue in trying to load file : '+cf
print e
sys.exit(42)


# build dataframe
self.ts = pd.DataFrame(pd.concat(df)).sort_index()
self.mean = self.ts[self.name].mean()
self.std = self.ts[self.name].std()
self.min = self.ts[self.name].min()
self.max = self.ts[self.name].max()

def __str__(self):
return 'runid = {}, name = {}, line = {}, color = {}'.format(self.runid, self.name, self.line, self.color)
Expand Down Expand Up @@ -108,9 +114,10 @@ def find_key(char, fid):
def load_argument():
# deals with argument
parser = argparse.ArgumentParser()
parser.add_argument("-runid", metavar='runid list' , help="used to look information in runid.db" , type=str, nargs='+' , required=True)
parser.add_argument("-runid", metavar='runid list' , help="used to look information in runid.db" , type=str, nargs='+' , required=True )
parser.add_argument("-f" , metavar='file list' , help="file list to plot (default is runid_var.nc)" , type=str, nargs='+' , required=False)
parser.add_argument("-var" , metavar='var list' , help="variable to look for in the netcdf file ./runid_var.nc", type=str, nargs='+' , required=True)
parser.add_argument("-var" , metavar='var list' , help="variable to look for in the netcdf file ./runid_var.nc", type=str, nargs='+' , required=True )
parser.add_argument("-varf" , metavar='var list' , help="variable to look for in the netcdf file ./runid_var.nc", type=str, nargs='+' , required=False)
parser.add_argument("-title", metavar='title' , help="subplot title (associated with var)" , type=str, nargs='+' , required=False)
parser.add_argument("-dir" , metavar='directory of input file' , help="directory of input file" , type=str, nargs=1 , required=False, default=['./'])
parser.add_argument("-sf" , metavar='scale factor', help="scale factor" , type=float, nargs=1 , required=False, default=[1])
Expand Down Expand Up @@ -225,6 +232,11 @@ def main():
plt.figure(figsize=np.array([210, 210]) / 25.4)

# need to deal with multivar
mintime=dt.date.max
maxtime=dt.date.min
ymin=-sys.float_info.max
ymax=sys.float_info.max

for ivar, cvar in enumerate(args.var):
ax[ivar] = plt.subplot(nvar, 1, ivar+1)
# load obs
Expand All @@ -246,6 +258,16 @@ def main():
if len(cfile)==0:
print 'no file found with this pattern '+args.dir[0]+'/'+runid+'/'+fglob
sys.exit(42)
elif args.varf:
# in case only one file pattern given
if len(args.varf) == 1 :
fglob = args.varf[0]
else :
fglob = args.varf[ivar]
cfile = glob.glob(args.dir[0]+'/'+runid+'/'+fglob)
if len(cfile)==0:
print 'no file found with this pattern '+args.dir[0]+'/'+runid+'/'+fglob
sys.exit(42)
else:
cfile = glob.glob(args.dir[0]+'/'+runid+'_'+cvar+'.nc')
if len(cfile)==0:
Expand All @@ -254,14 +276,32 @@ def main():

run_lst[irun].load_time_series(cfile, cvar)
ts_lst[irun] = run_lst[irun].ts
lg = ts_lst[irun].plot(ax=ax[ivar], legend=False, style=run_lst[irun].line,color=run_lst[irun].color,label=run_lst[irun].name, linewidth=2, rot=0)
lg = ts_lst[irun].plot(ax=ax[ivar], legend=False, style=run_lst[irun].line,color=run_lst[irun].color,label=run_lst[irun].name, x_compat=True, linewidth=2, rot=0)
#
# limit of time axis
mintime=min([mintime,ts_lst[irun].index[0].to_pydatetime().date()])
maxtime=max([maxtime,ts_lst[irun].index[-1].to_pydatetime().date()])

# set title
if (args.title):
ax[ivar].set_title(args.title[ivar],fontsize=20)

# rm xaxis label
ax[ivar].xaxis.set_major_locator(mdates.YearLocator())

# set x axis
nlabel=5
ndays=(maxtime-mintime).days
nyear=ndays/365
if nyear < 10:
nyt=1
elif 10<=nyear<50:
nyt=5
elif 50<=nyear<100:
nyt=10
else:
nyt=100
nmt=ts_lst[irun].index[-1].to_pydatetime().date().month
ndt=ts_lst[irun].index[-1].to_pydatetime().date().day

ax[ivar].xaxis.set_major_locator(mdates.YearLocator(nyt,month=nmt,day=1))
ax[ivar].tick_params(axis='both', labelsize=16)
if (ivar != nvar-1):
ax[ivar].set_xticklabels([])
Expand All @@ -270,9 +310,9 @@ def main():

for lt in ax[ivar].get_xticklabels():
lt.set_ha('center')
# ax[ivar].set_xticklabels(ha='center')

rmin[ivar],rmax[ivar]=get_ybnd(run_lst,obs_min[ivar],obs_max[ivar])
print rmin[ivar], rmax[ivar]
ax[ivar].set_ylim([rmin[ivar],rmax[ivar]])
ax[ivar].grid()

Expand Down
29 changes: 29 additions & 0 deletions param.bash
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,4 +12,33 @@ EXEPATH=${HOME}/VALSO/
# SCRIPT location
SCRPATH=${HOME}/VALSO/SCRIPT/

#RUNVALSO=0
#RUNVALGLO=0
#RUNALL=1
#
#if [[ $RUNVALSO == 1 ]]; then
# runACC=1 #acc ts
# runMLD=1 #mld ts
# runBSF=1 #gyre ts
# runBOT=1 #bottom TS ts
#elif [[ $RUNVALGLO == 1 ]]; then
# runAMOC=1
# runAMHT=1
# runSIE=1
# runSST=1
#elif [[ $RUNALL == 1 ]]; then
# runACC=1 #acc ts
# runMLD=1 #mld ts
# runBSF=1 #gyre ts
# runBOT=1 #bottom TS ts
# runAMOC=1
# runAMHT=1
# runSIE=1
# runSST=1
#else
# echo 'need to define what you want in param.bash; exit 42'
# exit 42
#fi


module load gcc/8.1.0 mpi/mpich/3.2.1/gnu/8.1.0 hdf5/1.8.20/gnu/8.1.0 netcdf/4.6.1/gnu/8.1.0
15 changes: 11 additions & 4 deletions run_all.bash
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,13 +28,15 @@ for RUNID in `echo $RUNIDS`; do
TAG=${YEAR}1201-$((YEAR+1))1201
TAG09=${YEAR}0901-${YEAR}1001
TAG02=${YEAR}0201-${YEAR}0301
TAG02=${YEAR}0301-${YEAR}0401

# get data
mooVyid=$(sbatch --job-name=moo_${YEAR}_V --output=${JOBOUT_PATH}/moo_${YEAR}_V ${SCRPATH}/get_data.bash $CONFIG $RUNID 1y $TAG grid-V | awk '{print $4}') # for mk_trp and mk_psi
mooUyid=$(sbatch --job-name=moo_${YEAR}_U --output=${JOBOUT_PATH}/moo_${YEAR}_U ${SCRPATH}/get_data.bash $CONFIG $RUNID 1y $TAG grid-U | awk '{print $4}') # for mk_trp and mk_psi
mooTyid=$(sbatch --job-name=moo_${YEAR}_T --output=${JOBOUT_PATH}/moo_${YEAR}_T ${SCRPATH}/get_data.bash $CONFIG $RUNID 1y $TAG grid-T | awk '{print $4}') # for mk_bot.bash
mooT09mid=$(sbatch --job-name=moo_${YEAR}_T --output=${JOBOUT_PATH}/moo_${YEAR}_T09 ${SCRPATH}/get_data.bash $CONFIG $RUNID 1m $TAG09 grid-T | awk '{print $4}') # for mk_mxl.bash
mooT02mid=$(sbatch --job-name=moo_${YEAR}_T --output=${JOBOUT_PATH}/moo_${YEAR}_T02 ${SCRPATH}/get_data.bash $CONFIG $RUNID 1m $TAG02 grid-T | awk '{print $4}') # for mk_mxl.bash
mooT03mid=$(sbatch --job-name=moo_${YEAR}_T --output=${JOBOUT_PATH}/moo_${YEAR}_T03 ${SCRPATH}/get_data.bash $CONFIG $RUNID 1m $TAG03 grid-T | awk '{print $4}') # for mk_mxl.bash

# run cdftools
# scheduler option
Expand All @@ -54,15 +56,15 @@ for RUNID in `echo $RUNIDS`; do

# mld
# VALSO
sbatchrunopt="--dependency=afterany:$mooTmid --job-name=SO_mxl_${TAG}_${RUNID} --output=${JOBOUT_PATH}/mxl_${TAG}.out"
sbatchrunopt="--dependency=afterany:$mooT09mid --job-name=SO_mxl_${TAG}_${RUNID} --output=${JOBOUT_PATH}/mxl_${TAG}.out"
sbatch ${sbatchschopt} ${sbatchrunopt} ${SCRPATH}/mk_mxl.bash $CONFIG $RUNID $TAG09 1m > /dev/null 2>&1 &
njob=$((njob+1))

# botT/S
# VALSO
#sbatchrunopt="--dependency=afterany:$mooTyid --job-name=SO_bot_${TAG}_${RUNID} --output=${JOBOUT_PATH}/bot_${TAG}.out"
#sbatch ${sbatchschopt} ${sbatchrunopt} ${SCRPATH}/mk_bot.bash $CONFIG $RUNID $TAG 1y > /dev/null 2>&1 &
#njob=$((njob+1))
sbatchrunopt="--dependency=afterany:$mooTyid --job-name=SO_bot_${TAG}_${RUNID} --output=${JOBOUT_PATH}/bot_${TAG}.out"
sbatch ${sbatchschopt} ${sbatchrunopt} ${SCRPATH}/mk_bot.bash $CONFIG $RUNID $TAG 1y > /dev/null 2>&1 &
njob=$((njob+1))

# moc
# VALGLO
Expand Down Expand Up @@ -94,6 +96,11 @@ for RUNID in `echo $RUNIDS`; do
sbatch ${sbatchschopt} ${sbatchrunopt} ${SCRPATH}/mk_sie.bash $CONFIG $RUNID $TAG02 1m > /dev/null 2>&1 &
njob=$((njob+1))

# VALGLO/VALSO
sbatchrunopt="--dependency=afterany:$mooT03mid --job-name=GLO_sie_${TAG}_${RUNID} --output=${JOBOUT_PATH}/sie_${TAG03}.out"
sbatch ${sbatchschopt} ${sbatchrunopt} ${SCRPATH}/mk_sie.bash $CONFIG $RUNID $TAG03 1m > /dev/null 2>&1 &
njob=$((njob+1))

# sie
# VALGLO/VALSO
sbatchrunopt="--dependency=afterany:$mooT09mid --job-name=GLO_sie_${TAG}_${RUNID} --output=${JOBOUT_PATH}/sie_${TAG09}.out"
Expand Down
6 changes: 4 additions & 2 deletions run_plot_VALGLO.bash
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
#!/bin/bash

set -x

if [ $# -eq 0 ] ; then echo 'need a [KEYWORD] (will be inserted inside the figure title and output name) and a list of id [RUNIDS RUNID ...] (definition of line style need to be done in RUNID.db)'; exit; fi

module load scitools/production-os41-1
Expand Down Expand Up @@ -48,12 +50,12 @@ convert ${KEY}_fig06.png -crop 1240x1040+0+0 tmp06.png

# ARC SIE
echo 'plot 02 SIE time series'
python2.7 SCRIPT/plot_time_series.py -noshow -runid $RUNIDS -f '*sie*0201*.nc' -var NExnsidc SExnsidc -sf 0.001 -title "SIE Arctic m02 [1e6 km2] : ${KEY}" "SIE Antarctic m02 [1e6 km2] : ${KEY}" -dir ${DATPATH} -o ${KEY}_fig07 -obs OBS/ARC_sie02_obs.txt OBS/ANT_sie02_obs.txt
python2.7 SCRIPT/plot_time_series.py -noshow -runid $RUNIDS -varf '*sie*0301*.nc' '*sie*0901*.nc' -var NExnsidc NExnsidc -sf 0.001 -title "SIE Arctic m03 [1e6 km2] : ${KEY}" "SIE Arctic m09 [1e6 km2] : ${KEY}" -dir ${DATPATH} -o ${KEY}_fig07 -obs OBS/ARC_sie03_obs.txt OBS/ARC_sie09_obs.txt
if [[ $? -ne 0 ]]; then exit 42; fi
convert ${KEY}_fig07.png -crop 1240x1040+0+0 tmp07.png

echo 'plot 09 SIE time series'
python2.7 SCRIPT/plot_time_series.py -noshow -runid $RUNIDS -f '*sie*0901*.nc' -var NExnsidc SExnsidc -sf 0.001 -title "SIE Arctic m09 [1e6 km2] : ${KEY}" "SIE Antarctic m09 [1e6 km2] : ${KEY}" -dir ${DATPATH} -o ${KEY}_fig08 -obs OBS/ARC_sie09_obs.txt OBS/ANT_sie09_obs.txt
python2.7 SCRIPT/plot_time_series.py -noshow -runid $RUNIDS -varf '*sie*0201*.nc' '*sie*0901*.nc' -var SExnsidc SExnsidc -sf 0.001 -title "SIE Antarctic m02 [1e6 km2] : ${KEY}" "SIE Antarctic m09 [1e6 km2] : ${KEY}" -dir ${DATPATH} -o ${KEY}_fig08 -obs OBS/ANT_sie02_obs.txt OBS/ANT_sie09_obs.txt
if [[ $? -ne 0 ]]; then exit 42; fi
convert ${KEY}_fig08.png -crop 1240x1040+0+0 tmp08.png

Expand Down

0 comments on commit eec9991

Please sign in to comment.