source: lmdz_wrf/trunk/tools/obs-sim_Comparison.bash @ 1982

Last change on this file since 1982 was 1975, checked in by lfita, 7 years ago

Adding:

  • Shell script to compare observations with simulations using PyNCplot
  • Property svn:executable set to *
File size: 23.9 KB
Line 
1#!/bin/bash
2# Shell script to compare observations with simulations using PyNCplot
3#
4#   observations: different files from create_OBSnetcdf.py and/or UWyoming_snd_nc.py
5#   simulations: multiple outputs from different models and runs
6#
7
8# files from scratch
9fscratch=false
10
11# figures from scratch
12gscratch=true
13
14# Working directory
15wdir='/home/lluis/estudios/ChemGBsAs/tests/199501'
16
17## U. Wyoming file (respect wdir)
18snddir='obs/snd'
19sndorigfn='UWyoming_snd_87576.nc'
20snddimt=6
21sndstn='Ezeiza'
22
23# Kind of graphical file output
24kfig='png'
25
26# Color for SkewT_evol
27col=255
28col2=200
29fraccol=16
30fraccol2=16
31
32# Simulation dir
33simdir='/home/lluis/PY'
34
35# Variables from the simulations
36simvars='XLONG,XLAT,T,P,PB,PH,PHB,QVAPOR,U,V,MU,MUB,SINALPHA,COSALPHA,T2,U10,V10,'
37simvars=${simvars}'PSFC,Q2'
38
39# Diagnotics for the soundigns
40wrfsnddiags='WRFt:WRFtda:WRFp:WRFua:WRFva:ws'
41
42# Figures output directory
43ofigdir=${wdir}/figs
44
45# Sounding comparison vars
46sndvars='ta:tda:ws'
47
48# Observations surface files
49sfcobsdir='obs/sfc_CAM/'
50obsHf='OBSnetcdf_'
51# Surface variables [varn1]@[minv1]@[maxv1]:[varn2]@[minv2]@[maxv2]:...
52#sfcvars='tas:wds:wss:psl:clt:tdas'
53sfcvars='tas@FROMobs@FROMobs:wds@0.@6.28:wss@0.@@FROMobs:tdas@FROMobs@FROMobs'
54#wrfsfcdiags='wds:wss:WRFpsl_ecmwf:clt:WRFtdas'
55wrfsfcdiags='T2:wds:wss:WRFtdas'
56
57# map values
58mapv='cyl,f'
59
60# map coverage for cont-disc
61mapcover='sponge,0.1,0.1'
62
63# Format of time-axis for time-series
64fmtTts='exct,12,h|%d$^{%H}$'
65
66####### ###### ##### #### ### ## #
67
68function num_hex(){
69# Fucntion to transform a number to a hexagonal basis
70#   num: number to transform
71  fname='num_hex'
72  num=$1
73
74  if test ${num} -gt 255; then
75    echo ${errormsg}
76    echo "  "${fname}": value "${num}" too large !!"
77    exit
78  fi
79  if test ${num} -gt 16; then
80    xviA=`expr ${num} / 16`
81    numA=`expr ${xviA} '*' 16`
82    xviB=`expr ${num} - ${numA}`
83  else
84    xviA=0
85    xviB=${num}
86  fi
87
88  if test ${xviA} -gt 9; then
89    case ${xviA} in
90      10)
91        hexA='A'
92      ;;
93      11)
94        hexA='B'
95      ;;
96      12)
97        hexA='C'
98      ;;
99      13)
100        hexA='D'
101      ;;
102      14)
103        hexA='E'
104      ;;
105      15)
106        hexA='F'
107      ;;
108    esac
109  else
110    hexA=${xviA}
111  fi
112
113  if test ${xviB} -gt 9; then
114    case ${xviB} in
115      10)
116        hexB='A'
117      ;;
118      11)
119        hexB='B'
120      ;;
121      12)
122        hexB='C'
123      ;;
124      13)
125        hexB='D'
126      ;;
127      14)
128        hexB='E'
129      ;;
130      15)
131        hexB='F'
132      ;;
133    esac
134  else
135    hexB=${xviB}
136  fi
137
138  echo ${hexA}${hexB}
139}
140
141#######    #######
142## MAIN
143    #######
144ofilefigs=${wdir}/allins_figures.inf
145echo "*** INstruction from all figures" > ${ofilefigs}
146
147mkdir -p ${ofigdir}
148
149# Soundings plots
150cd ${wdir}/${snddir}
151
152# Time characteristics
153tunits=`python ${pyHOME}/nc_var.py -o ivattrs -v time -f ${sndorigfn} | grep -v allvattrs | grep units | awk '{print $3}'`
154Tu=`echo ${tunits} | tr '!' ' ' | awk '{print $1}'`
155Dref0=`echo ${tunits} | tr '!' ' ' | awk '{print $3}'`
156Dref=${Dref0:0:4}${Dref0:5:2}${Dref0:8:2}
157
158# Sounding
159
160# Getting wind components
161ofile='sounding_uava.nc'
162if ${fscratch}; then rm ${ofile}; fi
163if test ! -f ${ofile}; then
164  python ${pyHOME}/diagnostics.py -f ${sndorigfn} -d 'time@time,pres@pres' -v 'uavaFROMwswd|ws@wd'
165  if test $? -ne 0; then
166    echo ${errmsg}
167    echo "  python failed!!"
168    echo python ${pyHOME}/diagnostics.py -f ${sndorigfn} -d 'time@time,pres@pres' -v 'uavaFROMwswd|ws@wd'
169    exit
170  fi 
171  mv diagnostics.nc ${ofile}
172fi
173
174it=0
175tatda_evol_values=''
176tatda_evol_files=''
177sndtimes=''
178snddates=''
179sndFdates=''
180while test ${it} -le ${snddimt}; do
181  timev=`python $pyHOME/nc_var.py -o varout -S time:${it} -f ${sndorigfn} -v time | grep NC | awk '{printf ("%d", $2)}'`
182  timeS=`date +%Y"/"%m"/"%d"!"%H -d"${Dref} ${timev} ${Tu}"`
183  timeS2=`date +%d"$^{"%H"}$" -d"${Dref} ${timev} ${Tu}"`
184  timeS3=`date +%Y"-"%b -d"${Dref} ${timev} ${Tu}"`
185  timeS4=`date +%Y -d"${Dref} ${timev} ${Tu}"`
186  echo ${it}" "${timev}" "${Dref}" '"${timev}"' "${Tu}" "$(echo ${timeS} | tr '!' ' ')
187  timefS=`date +%Y%m%d%H%M%S -d"${Dref} ${timev} ${Tu}"`
188  # Wind rose
189  values='pres|-1;time|'${it}':linepoint;multicol;pres;auto;auto;rainbow;auto:'
190  values=${values}${sndstn}'!sounding!WindRose!on!'${timeS}'!local!time:png:cardinals:False:WindRose:True'
191  ofign=${ofigdir}'/WindRose_'${timefS}'.png'
192  if ${gscratch}; then rm ${ofign} >& /dev/null; fi
193  if test ! -f ${ofign}; then
194    echo "  Wind Rose on: "${timeS}
195    python ${pyHOME}/drawing.py -o draw_WindRose -S ${values} -v ua,va -f ${ofile} >& /dev/null
196    if test $? -ne 0; then
197      echo ${errmsg}
198      echo "  python failed!!"
199      echo python ${pyHOME}/drawing.py -o draw_WindRose -S ${values} -v ua,va -f ${ofile}
200      exit
201    fi
202    mv WindRose.png ${ofign}
203    echo "* "${ofign} >> ${ofilefigs}
204    echo python ${pyHOME}/drawing.py -o draw_WindRose -S ${values} -v ua,va -f ${ofile} >> ${ofilefigs}
205  fi
206
207  # Sounding
208  echo ${sndorigfn}
209  values='time|'${it}',pres|-1:auto:auto:'${sndstn}'!sounding!on!'${timeS}'!local!time:png:True'
210  ofign=${ofigdir}'/SkewT-logP_'${timefS}'.png'
211  if ${gscratch}; then rm ${ofign} >& /dev/null; fi
212  if test ! -f ${ofign}; then
213    echo "  Sounding on: "$(echo ${timeS} | tr '!' ' ')
214    python ${pyHOME}/drawing.py -o draw_SkewT -S ${values} -v ta,tda,pres -f ${sndorigfn} >& /dev/null
215    if test $? -ne 0; then
216      echo ${errmsg}
217      echo "  python failed!!"
218      echo python ${pyHOME}/drawing.py -o draw_SkewT -S ${values} -v ta,tda,pres -f ${sndorigfn}
219      exit
220    fi
221    mv SkewT.png ${ofign}
222    echo "* "${ofign} >> ${ofilefigs}
223    echo python ${pyHOME}/drawing.py -o draw_SkewT -S ${values} -v ta,tda,pres -f ${sndorigfn} >> ${ofilefigs}
224  fi
225
226  colH=`num_hex ${col}`
227  colH2=`num_hex ${col2}`
228  if test ${it} -eq 0; then
229    sndtimes=${timev}
230    snddates=${timeS}
231    sndFdates=${timefS}
232    tatda_evol_files=${sndorigfn}':pres|-1,time|'${it}':ta,pres;'
233    tatda_evol_files=${tatda_evol_files}${sndorigfn}':pres|-1,time|'${it}':tda,pres'
234    tatda_evol_labs='ta_'${timeS2}',tda'
235    tatda_evol_cols='#'${colH}${colH2}${colH2}',#'${colH2}${colH2}${colH}
236  else
237    sndtimes=${sndtimes}','${timev}
238    snddates=${snddates}','${timeS}
239    sndFdates=${sndFdates}','${timefS}
240    tatda_evol_files=${tatda_evol_files}';'${sndorigfn}':pres|-1,time|'${it}':ta,pres'
241    tatda_evol_files=${tatda_evol_files}';'${sndorigfn}':pres|-1,time|'${it}':tda,pres'
242    if test $(expr ${it} % 2) -eq 0; then 
243      tatda_evol_labs=${tatda_evol_labs}','${timeS2}',None'
244    else 
245      tatda_evol_labs=${tatda_evol_labs}',None,None'
246    fi
247    tatda_evol_cols=${tatda_evol_cols}',#'${colH}${colH2}${colH2}',#'${colH2}${colH2}${colH}
248  fi
249  col=`expr ${col} - ${fraccol}`
250  col2=`expr ${col2} - ${fraccol2}`
251
252  #exit
253
254  it=`expr ${it} + 1`
255done
256
257# multi_sounding
258tatda_evol_values='auto:auto:multilines!'${tatda_evol_labs}'!'${tatda_evol_cols}'!'
259tatda_evol_values=${tatda_evol_values}'-!,!2:0,auto:!sounding!evolution!on!'
260tatda_evol_values=${tatda_evol_values}${timeS3}'!local!time:png:True'
261ofign=${ofigdir}'/SkewT-logP_evol.png'
262if ${gscratch}; then rm ${ofign} >& /dev/null; fi
263if test ! -f ${ofign}; then
264  echo "  Sounding on: "$(echo ${timeS} | tr '!' ' ')
265  python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${tatda_evol_values} -f ${tatda_evol_files} >& /dev/null
266  if test $? -ne 0; then
267    echo ${errmsg}
268    echo "  python failed!!"
269    python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${tatda_evol_values} -f ${tatda_evol_files}
270    exit
271  fi
272  mv multi_SkewT.png ${ofign}
273  echo "* "${ofign} >> ${ofilefigs}
274  echo python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${tatda_evol_values} -f ${tatda_evol_files} >> ${ofilefigs}
275fi
276
277# comparison
278cd ${wdir}
279sndlon=`ncdump -h ${wdir}/${snddir}/${sndorigfn} | grep Station_longitude | awk '{print $3}'`
280sndlat=`ncdump -h ${wdir}/${snddir}/${sndorigfn} | grep Station_latitude | awk '{print $3}'`
281sndts=`echo ${sndtimes} | tr ',' ' '`
282
283# Joining files
284ofile0='simout_vars.nc'
285if test ! -f ${ofile0}; then
286  values=${simdir}',Time,WRFtime'
287  HMT='wrfout,'${timeS4},'00'
288  python ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S ${values} -f ${HMT} \
289   -v ${simvars}
290  if test $? -ne 0; then
291    echo ${errmsg}
292    echo "  python failed!!"
293    echo python ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S ${values}    \
294    -f ${HMT} -v ${simvars}
295    exit
296  fi
297  mv netcdf_fold_concatenated_HMT.nc ${wdir}/${ofile0}
298fi
299
300gridsnd=`python ${pyHOME}/nc_var.py -o get_point -f ${wdir}/${ofile0}                \
301  -S 'XLONG:XLAT:Time|0' -v ${sndlon}','${sndlat}`
302xsnd=`echo ${gridsnd} | tr ',' ' ' | awk '{print $1}'`
303ysnd=`echo ${gridsnd} | tr ',' ' ' | awk '{print $2}'`
304echo "Equivalent simulated sounding grid point: "${gridsnd}
305
306# Computing sounding diagnostics
307diagns=`echo ${wrfsnddiags} | tr ':' ' '`
308diagd='time@time,bottom_top@P,south_north@XLAT,west_east@XLONG'
309diagvals=''
310idiag=1
311for diagn in ${diagns}; do
312  diagv=`cat $pyHOME/diagnostics.inf | grep ${diagn}',' | tr ',' ' ' |               \
313    awk '{print $2"|"$3}'`
314  if test ${idiag} -eq 1; then
315    diagvals=${diagv}
316  else
317    diagvals=${diagvals}','${diagv}
318  fi
319  idiag=`expr ${idiag} + 1`
320done
321
322ofile='simout_snddiags.nc'
323if test ! -f ${ofile}; then
324  python ${pyHOME}/diagnostics.py -f ${ofile0} -d ${diagd} -v ${diagvals}
325  if test $? -ne 0; then
326    echo ${errmsg}
327    echo "  python failed!!"
328    echo python ${pyHOME}/diagnostics.py -f ${ofile0} -d ${diagd} -v ${diagvals}
329    exit
330  fi
331  mv diagnostics.nc simout_snddiags.nc
332  #To deg
333  python ${pyHOME}/nc_var.py -o valmod -S subc,273.15 -f ${ofile} -v ta
334  #To Pa
335  python ${pyHOME}/nc_var.py -o valmod -S divc,100. -f ${ofile} -v pres
336fi
337
338# Getting values at the sounding point
339ofile2='simout_vars_sndpt.nc'
340if test ! -f ${ofile2}; then
341  values='time,0,-1,1@bottom_top,0,-1,1@south_north,'${ysnd}','${ysnd}',1@'
342  values=${values}'west_east,'${xsnd}','${xsnd}',1'
343  python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${wdir}/${ofile}         \
344    -S ${values} -v all
345  if test $? -ne 0; then
346    echo ${errmsg}
347    echo "  python failed!!"
348    echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${wdir}/${ofile}  \
349    -S ${values} -v all
350    exit
351  fi
352  ofilen=`ls -rt1 | tail -n 1`
353  mv ${ofilen} ${ofile2}
354fi
355
356# Time information from simulation
357simdimt=`python ${pyHOME}/nc_var.py -o itime -S CFtime, -v time -f simout_snddiags.nc\
358  | grep dimt | awk '{print $2}'`
359simtunits=`python ${pyHOME}/nc_var.py -o ivattrs -v time -f simout_snddiags.nc |     \
360  grep -v allvattrs | grep units | awk '{print $3}'`
361simTu=`echo ${tunits} | tr '!' ' ' | awk '{print $1}'`
362simDref0=`echo ${tunits} | tr '!' ' ' | awk '{print $3}'`
363simDref=${simDref0:0:4}${simDref0:5:2}${simDref0:8:2}
364
365# Plotting
366complabs='$ta^{obs}$,$ta^{sim}$,$tda^{obs}$,$tda^{sim}$'
367compcols='#FF6464,#FF6464,#6464FF,#6464FF'
368it=0
369for its in ${sndts}; do
370  timestep=`python ${pyHOME}/nc_var.py -o get_time -f ${wdir}/${ofile2}              \
371    -S ${its}';'${tunits} -v time`
372  if test $? -ne 0; then
373    echo ${errmsg}
374    echo "  python failed!!"
375    echo python ${pyHOME}/nc_var.py -o get_time -f ${wdir}/${ofile2}                 \
376      -S ${its}';'${tunits} -v time
377  fi
378
379  dateS=`echo ${snddates} | tr ',' '\n' | head -n ${it} | tail -n 1`
380  datefS=`echo ${sndFdates} | tr ',' '\n' | head -n ${it} | tail -n 1`
381  echo ${its}": it="${timestep}" "$(echo ${dateS} | tr '!' ' ')" "${datefS}
382
383  values='auto:auto:multilines!'${complabs}'!'${compcols}'!'
384  values=${values}'-,-.,-,-.!,!2:0,auto:!sounding!evolution!on!'
385  values=${values}${dateS}'!local!time:png:True'
386
387  compfiles=${wdir}/${snddir}/${sndorigfn}':pres|-1,time|'${it}':ta,pres;'
388  compfiles=${compfiles}${ofile2}':bottom_top|-1,time|'${timestep}':ta,pres;'
389  compfiles=${compfiles}${wdir}/${snddir}/${sndorigfn}':pres|-1,time|'${it}':tda,pres;'
390  compfiles=${compfiles}${ofile2}':bottom_top|-1,time|'${timestep}':tda,pres'
391
392  ofign=${ofigdir}'/SkewT-logP_obs-sim_'${datefS}'.png'
393  if ${gscratch}; then rm ${ofign} >& /dev/null; fi
394  if test ! -f ${ofign}; then
395    echo "  Sounding on: "$(echo ${dateS} | tr '!' ' ')
396    python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${values} -f ${compfiles} >& /dev/null
397    if test $? -ne 0; then
398      echo ${errmsg}
399      echo "  python failed!!"
400      echo python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${values} -f ${compfiles}
401      exit
402    fi
403    mv multi_SkewT.png ${ofign}
404    echo "* "${ofign} >> ${ofilefigs}
405    echo python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${values} -f ${compfiles} >> ${ofilefigs}
406
407  fi
408  it=`expr ${it} + 1`
409done
410
411# Plotting snd comparison
412
413# Soundings
414sndvs=`echo ${sndvars} | tr ':' ' '`
415for sndv in ${sndvs}; do
416  ofign=${ofigdir}'/SkewT-logP_obs-sim_evol_'${sndv}'.png'
417
418  cfiles='simout_vars_sndpt.nc;'${sndv}';time;pres;time|-1,pres|-1@'
419  cfiles=${cfiles}'obs/snd/UWyoming_snd_87576.nc;'${sndv}';time;pres;time|-1,pres|-1'
420  cvalues=${sndv}';y;auto|exct,12,h|%d$^{%H}$|time!($[DD]^{[HH]}$);Vfix,auto,50.,auto'
421  cvalues=${cvalues}';auto;Srange,Srange;auto;obs!&!sim!Ezeiza!airport!sounding;png;'
422  cvalues=${cvalues}'flip@y;None;yes'
423  if ${gscratch}; then rm ${ofign} >& /dev/null; fi
424  if test ! -f ${ofign}; then
425    python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc_time -f ${cfiles} -S ${cvalues}
426    if test $? -ne 0; then
427      echo ${errmsg}
428      echo "  python failed!!"
429      echo python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc_time -f ${cfiles} -S ${cvalues}
430      exit
431    fi
432    mv 2Dshad_obs-sim_comparison_time.png ${ofign}
433    echo "* "${ofign} >> ${ofilefigs}
434    echo python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc_time -f ${cfiles} -S ${cvalues} >> ${ofilefigs}
435  fi
436done
437
438# Plotting map comparison
439cd $wdir
440
441# Joining all observations
442obssfcfile='all_single-stations.nc'
443if test ! -f ${wdir}/${sfcobsdir}/${obssfcfile}; then
444  python ${pyHOME}/nc_var.py -o join_singlestation_obsfiles -S                       \
445    ${sfcobsdir}':OBSnetcdf' -v all
446  if test $? -ne 0; then
447    echo ${errmsg}
448    echo "  python failed!!"
449    echo  python ${pyHOME}/nc_var.py -o join_singlestation_obsfiles \
450      -S ${sfcobsdir}':OBSnetcdf' -v all
451  fi
452  mv joined_singlestations.nc ${wdir}/${sfcobsdir}/${obssfcfile}
453fi
454
455# Computing surface diagnostics
456diagns=`echo ${wrfsfcdiags} | tr ':' ' '`
457diagd='time@time,south_north@XLAT,west_east@XLONG'
458diagvals=''
459idiag=1
460iadv=1
461for diagn in ${diagns}; do
462  diagv=`cat $pyHOME/diagnostics.inf | grep WRF | grep ${diagn}',' | tr ',' ' ' |    \
463    awk '{print $2"|"$3}'`
464  Ldiagv=`expr length $(echo ${diagv} | tr ' ' '@')0`
465  if test ${Ldiagv} -lt 2; then
466    echo ${errmsg}
467    echo "  diagnostic '"${diagn}"' (for WRF) does not exist !!"
468    if test ${iadv} -eq 1; then advs=${diagn}; iadv=`expr ${iadv} + 1`
469    else advs=${adv}','${diagn}; fi
470  else
471    if test ${idiag} -eq 1; then
472      diagvals=${diagv}
473      idiag=`expr ${idiag} + 1`
474    else
475      diagvals=${diagvals}','${diagv}
476    fi
477  fi
478done
479
480osfcdiagfile='simout_sfcdiags.nc'
481if test ! -f ${osfcdiagfile}; then
482  python ${pyHOME}/diagnostics.py -f ${ofile0} -d ${diagd} -v ${diagvals}
483  if test $? -ne 0; then
484    echo ${errmsg}
485    echo "  python failed!!"
486    echo python ${pyHOME}/diagnostics.py -f ${ofile0} -d ${diagd} -v ${diagvals}
487    exit
488  fi
489  mv diagnostics.nc ${osfcdiagfile}
490  #To rad
491  python ${pyHOME}/nc_var.py -o valmod -S mulc,0.0174532925199 -f ${osfcdiagfile} -v wds
492  python ${pyHOME}/nc_var.py -o varaddattr -S 'units|rad' -f ${osfcdiagfile} -v wds
493fi
494
495#Adding non-diagnostic variables
496if test ${iadv} -ne 1; then
497  varsadd=`echo ${advs} | tr ',' ' '`
498  for vadd in ${varsadd}; do
499    python ${pyHOME}/nc_var.py -o fvaradd -S ${ofile0},${vadd} -f ${osfcdiagfile}
500    if test $? -ne 0; then
501      echo ${errmsg}
502      echo "  python failed!!"
503      echo python ${pyHOME}/nc_var.py -o fvaradd -S ${osfcdiagfile},${vadd} -f ${ofile0}
504      exit
505      rm ${osfcdiagfile}
506    fi
507    # CF varname
508    CFvarn=`python $pyHOME/generic.py -o variables_values -S ${vadd} | tr ':' ' ' |  \
509      awk '{print $1}'`
510    if test $? -ne 0; then
511      echo ${errmsg}
512      echo "  python failed!!"
513      echo python $pyHOME/generic.py -o variables_values -S ${vadd}
514      exit
515      rm ${osfcdiagfile}
516    fi
517    python $pyHOME/nc_var.py -o chvarname -S ${CFvarn} -f ${osfcdiagfile} -v ${vadd}
518    if test $? -ne 0; then
519      echo ${errmsg}
520      echo "  python failed!!"
521      echo python ${pyHOME}/nc_var.py -o chvarname -S ${CFvarn} -f ${osfcdiagfile} -v ${vadd}
522      exit
523      rm ${osfcdiagfile}
524    fi
525  done
526fi
527
528sfcs=`echo ${sfcvars} | tr ':' ' '`
529
530# getting min,max values
531iv=1
532for sfcpv in ${sfcs}; do
533  sfcv=`echo ${sfcpv} | tr '@' ' ' | awk '{print $1}'`
534  nsfc=`echo ${sfcpv} | tr '@' ' ' | awk '{print $2}'`
535  xsfc=`echo ${sfcpv} | tr '@' ' ' | awk '{print $3}'`
536  if test ${nsfc} = 'FROMobs'; then
537    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
538      -f ${wdir}/${sfcobsdir}/${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
539    nsfc=`echo ${fstats} | awk '{print $3}'`
540  elif test ${nsfc} = 'FROMsims'; then
541    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
542      -f ${osfcdiagfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
543    nsfc=`echo ${fstats} | awk '{print $3}'`
544  elif test ${nsfc} = 'FROMobssims'; then
545    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
546      -f ${wdir}/${sfcobsdir}/${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
547    nsfc1=`echo ${fstats} | awk '{print $3}'`
548    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
549      -f ${osfcdiagfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
550    nsfc2=`echo ${fstats} | awk '{print $3}'`
551    echo ${nsfc1}" "${nsfc2} | awk '{if ($1 < $2) {print $1;} else {print $2;}}'
552  fi
553  if test ${xsfc} = 'FROMobs'; then
554    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
555      -f ${wdir}/${sfcobsdir}/${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
556    xsfc=`echo ${fstats} | awk '{print $4}'`
557  elif test ${xsfc} = 'FROMsims'; then
558    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
559      -f ${osfcdiagfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
560    xsfc=`echo ${fstats} | awk '{print $4}'`
561  elif test ${nsfc} = 'FROMobssims'; then
562    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
563      -f ${wdir}/${sfcobsdir}/${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
564    nsfc1=`echo ${fstats} | awk '{print $3}'`
565    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
566      -f ${osfcdiagfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
567    nsfc2=`echo ${fstats} | awk '{print $3}'`
568    echo ${nsfc1}" "${nsfc2} | awk '{if ($1 > $2) {print $1;} else {print $2;}}'
569  fi
570  if test $iv -eq 1; then
571    newsfcvars=${sfcv}'@'${nsfc}'@'${xsfc}
572  else
573    newsfcvars=${newsfcvars}':'${sfcv}'@'${nsfc}'@'${xsfc}
574  fi
575  iv=`expr ${iv} + 1`
576done
577sfcs=`echo ${newsfcvars} | tr ':' ' '`
578echo "new sfc variables: "${sfcs}
579
580# Time-series
581
582# stations
583# FROM: https://stackoverflow.com/questions/2961635/using-awk-to-print-all-columns-from-the-nth-to-the-last
584stsids=`python ${pyHOME}/nc_var.py -o varout -v id -f \
585  ${wdir}/${sfcobsdir}/${obssfcfile} -S 'Nstations:-1|time:0' | \
586  awk '{$1=""; printf("%d:",substr($0,2))}'`
587stsnames=`python ${pyHOME}/nc_var.py -o varout -v stsname -f \
588  ${wdir}/${sfcobsdir}/${obssfcfile} -S Nstations:-1 | \
589  awk '{$1=""; print substr($0,2)}' | tr '\n' ':' | tr ' ' '!'`
590stslons=`python ${pyHOME}/nc_var.py -o varout -v stslon -f \
591  ${wdir}/${sfcobsdir}/${obssfcfile} -S Nstations:-1 | \
592  awk '{$1=""; print substr($0,2)}' | tr '\n' ':'`
593stslats=`python ${pyHOME}/nc_var.py -o varout -v stslat -f \
594  ${wdir}/${sfcobsdir}/${obssfcfile} -S Nstations:-1 | \
595  awk '{$1=""; print substr($0,2)}' | tr '\n' ':'`
596
597Nsts=`echo ${stsnames} | tr ':' ' ' | wc -w | awk '{print $1}'`
598ist=1
599while test ${ist} -le ${Nsts}; do
600  sti=`echo ${stsids} | tr ':' '\n' | head -n ${ist} | tail -n 1`
601  stn=`echo ${stsnames} | tr ':' '\n' | head -n ${ist} | tail -n 1`
602  stl=`echo ${stslons} | tr ':' '\n' | head -n ${ist} | tail -n 1`
603  stL=`echo ${stslats} | tr ':' '\n' | head -n ${ist} | tail -n 1`
604  gridsfc=`python ${pyHOME}/nc_var.py -o get_point -f ${osfcdiagfile}              \
605    -S 'XLONG:XLAT:Time|0' -v ${stl}','${stL}`
606  xgrid=`echo ${gridsfc} | tr ',' ' ' | awk '{print $1}'`
607  ygrid=`echo ${gridsfc} | tr ',' ' ' | awk '{print $2}'`
608  sfcstats=${sti}'@'${stn}'@'${stl}'@'${stL}'@'${gridsfc}
609  ist_1=`expr ${ist} - 1`
610
611  for sfcpv in ${sfcs}; do
612    sfcv=`echo ${sfcpv} | tr '@' ' ' | awk '{print $1}'`
613    nsfc=`echo ${sfcpv} | tr '@' ' ' | awk '{print $2}'`
614    xsfc=`echo ${sfcpv} | tr '@' ' ' | awk '{print $3}'`
615
616    ofign=${ofigdir}'/obs-sim_evol_'${sfcv}_${sti}'.'${kfig}
617    files=${wdir}/${sfcobsdir}/${obssfcfile}'%time|-1;Nstations|'${ist_1}','
618    files=${files}${osfcdiagfile}'%time|-1;south_north|'${ygrid}';west_east|'${xgrid}
619    values='ststimes,time;y;time!([DD]$^{[HH]}$);auto;obs,sim;'${sfcv}';'${sfcv}
620    values=${values}'|evolution|at|'$(echo ${stn} | tr '!' '|')';'${nsfc}','${xsfc}';'
621    values=${values}'time|'${fmtTts}';0|12;'
622    values=${values}${kfig}';-;k,b;.;2.;2.;all;-1;True'
623    if ${gscratch}; then rm ${ofign} >& /dev/null; fi
624    if test ! -f ${ofign}; then
625      python ${pyHOME}/drawing.py -o draw_lines_time -f ${files} -S ${values} -v ${sfcv}
626      if test $? -ne 0; then
627        echo ${errmsg}
628        echo "  python failed!!"
629        python ${pyHOME}/drawing.py -o draw_lines_time -f ${files} -S ${values} -v ${sfcv}
630        exit
631      fi
632      mv lines_time.${kfig} ${ofign}
633      echo "* "${ofign} >> ${ofilefigs}
634      echo python ${pyHOME}/drawing.py -o draw_lines_time -f ${files} -S ${values} -v ${sfcv} >> ${ofilefigs}
635    fi
636    #exit
637  done
638  ist=`expr $ist + 1`
639  #exit
640done
641exit
642
643# Cont-disc comparison
644it=0
645while test ${it} -le ${simdimt}; do
646  simtstep=`python ${pyHOME}/nc_var.py -o varout -f ${wdir}/${osfcdiagfile} -S       \
647    'time:'${it} -v time | grep NC | awk '{printf ("%d",$2)}'`
648  if test $? -ne 0; then
649    echo ${errmsg}
650    echo "  python failed!!"
651    echo python ${pyHOME}/nc_var.py -o varout -f ${wdir}/${osfcdiagfile} -S          \
652      'time:'${it} -v time
653  fi
654
655  timeS=`date +%Y"/"%m"/"%d"!"%H -d"${simDref} ${simtstep} ${simTu}"`
656  timeS2=`date +%d"$^{"%H"}$" -d"${simDref} ${simtstep} ${simTu}"`
657  timeS3=`date +%Y"-"%b -d"${simDref} ${simtstep} ${simTu}"`
658  timeS4=`date +%Y -d"${simDref} ${simtstep} ${simTu}"`
659  echo ${it}" "${simtstep}" "${simDref}" "${simTu}" "$(echo ${timeS} | tr '!' ' ')
660  timefS=`date +%Y%m%d%H%M%S -d"${simDref} ${simtstep} ${simTu}"`
661
662  # Looking for obs time
663  obsf=${wdir}/${sfcobsdir}/${obsfile}
664  oit=`python ${pyHOME}/nc_var.py -o get_time -f ${wdir}/${sfcobsdir}/${obssfcfile}  \
665    -S ${simtstep}';'${tunits} -v ststimes | grep -v differ`
666
667  for sfcnxv in ${sfcs}; do
668    sfcv=`echo ${sfcnxv} | tr '@' ' ' | awk '{print $1}'`
669    nsfc=`echo ${sfcnxv} | tr '@' ' ' | awk '{print $2}'`
670    xsfc=`echo ${sfcnxv} | tr '@' ' ' | awk '{print $3}'`
671
672    ofign=${ofigdir}'/map_obs-sim_'${timefS}'_'${sfcv}'.png'
673
674    CFvarvals=`python $pyHOME/generic.py -o variables_values -S ${sfcv}`
675    cbar=`echo ${CFvarvals} | tr ':' ' ' | awk '{print $7}'`
676
677    cfiles=${wdir}/${osfcdiagfile}';'${sfcv}';XLONG;XLAT;Time|'${it}',time|'${it}','
678    cfiles=${cfiles}'west_east|-1,south_north|-1@'${wdir}/${sfcobsdir}/${obssfcfile}
679    cfiles=${cfiles}';'${sfcv}';stslon;stslat;time|'${oit}',lon|-1,lat|-1'
680    cvalues=${sfcv}':west_east,south_north:auto:'${cbar}',auto,auto:'${nsfc}','
681    cvalues=${cvalues}${xsfc}':''auto:obs!&!sim!'${sfcv}'!on!'${timeS}'!UTC:png:None:'
682    cvalues=${cvalues}${mapv}':'${mapcover}':yes'
683
684    if ${gscratch}; then rm ${ofign} >& /dev/null; fi
685    if test ! -f ${ofign}; then
686      python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc -f ${cfiles} -S ${cvalues}
687      if test $? -ne 0; then
688        echo ${errmsg}
689        echo "  python failed!!"
690        echo python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc -f ${cfiles} -S ${cvalues}
691      fi
692      mv 2Dshad_obs-sim_comparison.png ${ofign} 
693      echo "* "${ofign} >> ${ofilefigs}
694      echo python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc -f ${cfiles} -S ${cvalues} >> ${ofilefigs}
695
696    fi
697    #exit
698  done
699  #exit
700  it=`expr ${it} + 1`
701done
Note: See TracBrowser for help on using the repository browser.