DFS calculation updated to include SGNSS. Plot script updated to plot…

… also relative DFS (#46) Co-authored-by: Martin Ridal <smb@ac6-102.bullx>
ACCORD-NWP · May 17, 2024 · 7fb8169 · 7fb8169
1 parent a89cd3b
commit 7fb8169
Show file tree

Hide file tree

Showing 3 changed files with 93 additions and 21 deletions.
diff --git a/README.md b/README.md
@@ -1,5 +1,6 @@
 # AccordDaTools
 Collection of ACCORD DA tools. So far we have jbdiagnose and associated plotting. 
+Also available is DFS calculation and associated plotting.
 
 Desirables:
 - festat
@@ -59,3 +60,27 @@ Example:
 plotjbdiag -l 50 -t vercor -p QQ -r ./ -e IRELAND75L65 
 ```
 
+### DFS
+The dfscomp tool reads the (ASCII) data from an unperturbed CCMA and a perturbed CCMA.
+
+#### Calculate DFS
+Help/usage:
+```
+scripts/dfscomp.sh -h
+```
+
+Example:
+```
+scripts/dfscomp.sh -u $DTG_CCMA_unpert.dat -p $DTG_CCMA.dat -o dfs.dat 
+```
+
+#### Plot DFS
+Help/usage:
+```
+python3 scripts/plotdfs.py -h
+```
+
+Example:
+```
+python3 scripts/plotdfs.py -i dfs.dat
+```
diff --git a/scripts/plotdfs.py b/scripts/plotdfs.py
@@ -10,7 +10,7 @@
 import argparse
 
 #############################################################################################
-def plot_dfsabs_pie(infile) :
+def plot_dfsabs_pie(infile,dfsobsstr) :
 
   obindex,obnum,obdfsabs= np.loadtxt(infile, unpack=True)
 
@@ -22,13 +22,7 @@ def plot_dfsabs_pie(infile) :
 #    11 : AIREP_U #    12 : SATOB_U  #    13 : DRIBU_Z  #    14 : DRIBU_U   #    15 : PILOT_Z
 #    16 : PILOT_U #    17 : AMSUA    #    18 : MHS      #    19 : MWHS2     #    20 : ATMS
 #    21 : IASI    #    22 : CRIS     #    23 : SEVIRI   #    24 : SCATT_U   #    25 : RADAR_Z
-#    26 : RADAR_U #    27 : TEMP_CLS
-  dfsobsstr= ('SYNOP-Z','SYNOP-T2','SYNOP-R2','SYNOP-U10','GNSS-ZTD',
-              'TEMP-U','TEMP-T','TEMP-Z','TEMP-Q','AIREP-T',
-              'AIREP-U','SATOB-U','BUOY-Z','BUOY-U','PILOT-Z',
-              'PILOT-U','AMSU-A','MHS','MWHS','ATMS',
-              'IASI','CRIS','SEVIRI','SCATT-U','RADAR-Z',
-              'RADAR-U','TEMP_CLS')
+#    26 : RADAR_U #    27 : TEMP_CLS #    28 : SGNSS
 
   mask=obnum.nonzero()
 
@@ -45,7 +39,7 @@ def plot_dfsabs_pie(infile) :
   plt.savefig('dfsabs_pie.png',format='png')
 
 #############################################################################################
-def plot_dfsabs_bar(infile) :
+def plot_dfsabs_bar(infile,dfsobsstr) :
 
   obindex,obnum,obdfsabs= np.loadtxt(infile, unpack=True)
 
@@ -57,14 +51,7 @@ def plot_dfsabs_bar(infile) :
 #EW               11: 'AIREP-U', 12: 'SATOB-U', 13: 'BUOY-Z', 14: 'BUOY-U',   15: 'PILOT-Z',
 #EW               16: 'PILOT-U', 17: 'AMSU-A',  18: 'MHS',    19: 'MWHS',     20: 'ATMS',
 #EW               21: 'IASI',    22: 'CRIS',    23: 'SEVIRI', 24: 'SCATT-U',  25: 'RADAR-Z',
-#EW               26: 'RADAR-U', 27: 'TEMP_CLS'}
-
-  dfsobsstr= ('SYNOP-Z','SYNOP-T2','SYNOP-R2','SYNOP-U10','GNSS-ZTD',
-              'TEMP-U','TEMP-T','TEMP-Z','TEMP-Q','AIREP-T',
-              'AIREP-U','SATOB-U','BUOY-Z','BUOY-U','PILOT-Z',
-              'PILOT-U','AMSU-A','MHS','MWHS','ATMS',
-              'IASI','CRIS','SEVIRI','SCATT-U','RADAR-Z',
-              'RADAR-U','TEMP_CLS')
+#EW               26: 'RADAR-U', 27: 'TEMP_CLS', 28: 'SGNSS'}
 
   mask=obnum.nonzero()
 
@@ -78,7 +65,7 @@ def plot_dfsabs_bar(infile) :
 
   plt.xticks(strpos_mask, barlabs,rotation=90)
 
-  plt.ylabel('DFS',fontsize=10)
+  plt.ylabel('Absolute DFS',fontsize=10)
 
   # Hide the right and top spines
   ax.spines['right'].set_visible(False)
@@ -96,6 +83,51 @@ def plot_dfsabs_bar(infile) :
 
   plt.savefig('dfsabs_bar.png',format='png')
 
+#############################################################################################
+def plot_dfsrel_bar(infile,dfsobsstr) :
+
+  obindex,obnum,obdfsabs= np.loadtxt(infile, unpack=True)
+
+  fig = plt.figure()
+  ax = fig.add_subplot(111)
+
+#EW   dfsobsstr= { 1: 'SYNOP-Z', 2: 'SYNOP-T2',  3: 'SYNOP-R2', 4: 'SYNOP-U10', 5: 'GNSS-ZTD',
+#EW                6: 'TEMP-U',  7: 'TEMP-T',    8: 'TEMP-Z',   9: 'TEMP-Q',   10: 'AIREP-T',
+#EW               11: 'AIREP-U', 12: 'SATOB-U', 13: 'BUOY-Z', 14: 'BUOY-U',   15: 'PILOT-Z',
+#EW               16: 'PILOT-U', 17: 'AMSU-A',  18: 'MHS',    19: 'MWHS',     20: 'ATMS',
+#EW               21: 'IASI',    22: 'CRIS',    23: 'SEVIRI', 24: 'SCATT-U',  25: 'RADAR-Z',
+#EW               26: 'RADAR-U', 27: 'TEMP_CLS', 28: 'SGNSS'}
+
+  mask=obnum.nonzero()
+
+  barlabs=[]
+  for i in range(len(mask[0])):
+    j=mask[0][i]
+    barlabs.append(dfsobsstr[j])
+
+  strpos_mask=np.arange(len(barlabs))
+  plt.bar(strpos_mask, obdfsabs[mask]/obnum[mask], align='center', alpha=0.5)
+
+  plt.xticks(strpos_mask, barlabs,rotation=90)
+
+  plt.ylabel('Relative DFS',fontsize=10)
+
+  # Hide the right and top spines
+  ax.spines['right'].set_visible(False)
+  ax.spines['top'].set_visible(False)
+
+  # Only show ticks on the left and bottom spines
+  ax.yaxis.set_ticks_position('left')
+  ax.xaxis.set_ticks_position('bottom')
+
+  ax.xaxis.set_visible(True)
+  ax.yaxis.set_visible(True)
+  ax.yaxis.grid() #horizontal lines
+
+  plt.tight_layout()
+
+  plt.savefig('dfsrel_bar.png',format='png')
+
 #############################################################################################
 def main(argv) :
 
@@ -108,8 +140,18 @@ def main(argv) :
 
   args = parser.parse_args()
 
-  plot_dfsabs_bar(args.ipath)
-  plot_dfsabs_pie(args.ipath)
+
+  dfsobsstr= ('SYNOP-Z','SYNOP-T2','SYNOP-R2','SYNOP-U10','GNSS-ZTD',
+              'TEMP-U','TEMP-T','TEMP-Z','TEMP-Q','AIREP-T',
+              'AIREP-U','SATOB-U','BUOY-Z','BUOY-U','PILOT-Z',
+              'PILOT-U','AMSU-A','MHS','MWHS','ATMS',
+              'IASI','CRIS','SEVIRI','SCATT-U','RADAR-Z',
+              'RADAR-U','TEMP_CLS','SGNSS')
+
+
+  plot_dfsabs_bar(args.ipath,dfsobsstr)
+  plot_dfsrel_bar(args.ipath,dfsobsstr)
+  plot_dfsabs_pie(args.ipath,dfsobsstr)
 
 if __name__ == "__main__":
     sys.exit(main(sys.argv))
diff --git a/src/dfscomp.F90 b/src/dfscomp.F90
@@ -49,7 +49,7 @@ PROGRAM DFSCOMP
  LOGICAL            :: LEX
  INTEGER            :: INUM,IERR
  INTEGER            :: NUN1,NUN2,NUN3
- INTEGER, PARAMETER :: NUMV=27
+ INTEGER, PARAMETER :: NUMV=28
  REAL               :: ZDFS(NUMV)
  INTEGER            :: NDFS(NUMV)
  INTEGER            :: ITOT,INOTC,ICHAN
@@ -182,6 +182,7 @@ PROGRAM DFSCOMP
 !    25 : RADAR_Z
 !    26 : RADAR_U
 !    27 : TEMP_CLS
+!    28 : SGNSS
 
 
   IF (IOT1 .EQ. 1 ) THEN
@@ -266,6 +267,10 @@ PROGRAM DFSCOMP
       IF (IOV1.EQ.195) IIND=26
   ENDIF
 
+  ! SGNSS
+  IF (IOT1 .EQ. 19) IIND=28
+
+
   IF (IIND .EQ. 0 ) THEN
      INOTC = INOTC + 1
      ! This should not happen