radar intercomparison functions

config/processing/mals_emm_intercomp.txt

@@ -44,7 +44,7 @@ dempath STRARR 2
 # --------------------------------------------
 # Output data
 # -----------------------------------------------
-saveimgbasepath STRING /data/pyrad_products/
+saveimgbasepath STRING /data/pyrad_examples/
 saveimg         INT 1
 imgformat       STRING png
 

config/processing/paradiso_fvj_vol_intercomp_prod.txt

@@ -5,7 +5,7 @@
 # List of datasets to generate.
 # The detailed specification of each dataset is given below.
 dataSetList STRARR 1
-    DX50_PLD_intercomp
+    DX50_PLD_dBZ_avg_intercomp
 
 # ==========================================================================================
 #                 colocated gates
@@ -65,24 +65,51 @@ dataSetList STRARR 1
 # ==========================================================================================
 #                 intercomparison
 # ==========================================================================================
-DX50_PLD_intercomp STRUCT 8
-	type STRING INTERCOMP
-	datatype STRARR 2
+#DX50_PLD_intercomp STRUCT 8
+#	type STRING INTERCOMP
+#	datatype STRARR 2
+#		RADAR001:CFRADIAL:dBZc,dBZ_avg,SAVEVOL
+#        RADAR002:CFRADIAL:dBZc,dBZ_avg,SAVEVOL
+#    coloc_data_dir STRING DX50_PLD_DATA  # must be the same as product WRITE_INTERCOMP
+#    coloc_radars_name STRING DX50_PLD
+#    ele_tol FLOAT 0.5
+#    azi_tol FLOAT 0.5
+#    rng_tol FLOAT 100.
+#    products STRUCT 2        
+#        DX50_PLD_PLOT STRUCT 3
+#            type STRING PLOT_SCATTER_INTERCOMP
+#            voltype STRING dBZc
+#            step FLOAT 0.5
+#        DX50_PLD_DATA STRUCT 2
+#            type STRING WRITE_INTERCOMP
+#            voltype STRING dBZc
+
+DX50_PLD_dBZ_avg_intercomp STRUCT 12
+	type STRING INTERCOMP_TIME_AVG
+	datatype STRARR 6
 		RADAR001:CFRADIAL:dBZc,dBZ_avg,SAVEVOL
+        RADAR001:CFRADIAL:PhiDPc,PhiDP_avg,SAVEVOL
+        RADAR001:CFRADIAL:time_avg_flag,flag_avg,SAVEVOL
         RADAR002:CFRADIAL:dBZc,dBZ_avg,SAVEVOL
+        RADAR002:CFRADIAL:PhiDPc,PhiDP_avg,SAVEVOL
+        RADAR002:CFRADIAL:time_avg_flag,flag_avg,SAVEVOL
     coloc_data_dir STRING DX50_PLD_DATA  # must be the same as product WRITE_INTERCOMP
     coloc_radars_name STRING DX50_PLD
     ele_tol FLOAT 0.5
     azi_tol FLOAT 0.5
     rng_tol FLOAT 100.
-    products STRUCT 1        
+    clt_max INT 0
+    phi_excess_max INT 0
+    non_rain_max INT 100
+    phi_avg_max FLOAT 20.
+    products STRUCT 2        
         DX50_PLD_PLOT STRUCT 3
             type STRING PLOT_SCATTER_INTERCOMP
             voltype STRING dBZc
             step FLOAT 0.5
-#        DX50_PLD_DATA STRUCT 2
-#            type STRING WRITE_INTERCOMP
-#            voltype STRING dBZc
+        DX50_PLD_DATA STRUCT 2
+            type STRING WRITE_INTERCOMP_TIME_AVG
+            voltype STRING dBZc
 
 
 

doc/TODO_pyrad.txt

@@ -1,7 +1,6 @@
 Pending modifications for current code in Pyrad
 -----------------------------------------------
 - In function process_echo_id: Add condition on differential reflectivity at low reflectivity values. Allow user to select the various thresholds. Add a speckle filter
-- In selfconsistency functions: get tables as a function of elevation, allow user to define the tresholds
 - Functions for PhiDP and KDP processing consume too much memory up to the point that they cannot handle rad4alp PH data. It is probably due to the rolling window creating a monstrous matrix.
 We should try to reduce the memory consumption.
 
@@ -26,9 +25,6 @@ Things to check before making it operational:
     - Check if transmission losses have to be subtracted from received power
     - Check if conversion from solar flux to received power is done properly
 - Polarimetric data monitoring:
-    - Add metadata in the density and histogram plot
-    - Check what happens to data that is out of the range of the histogram
+    - Add metadata in the density and histogram plot    
 - Rainfall rate estimation:
-    - Distinguish between the various rainfall rate estimators (no single datatype RR)
-- Intercomparison:
-    - Add metadata in the plot
+    - Distinguish between the various rainfall rate estimators (no single datatype RR)

src/pyrad_proc/pyrad/flow/flow_control.py

@@ -521,14 +521,14 @@ def _get_datatype_list(cfg, radarnr='RADAR001'):
         proclevel, dataset = get_dataset_fields(datasetdescr)
         if 'datatype' in cfg[dataset]:
             if isinstance(cfg[dataset]['datatype'], str):
-                (radarnr_descr, datagroup, datatype, dataset_save,
+                (radarnr_descr, datagroup, datatype_aux, dataset_save,
                  product_save) = (
                     get_datatype_fields(cfg[dataset]['datatype']))
                 if datagroup != 'PROC' and radarnr_descr == radarnr:
                     if ((dataset_save is None) and (product_save is None)):
-                        datatypesdescr.add(datagroup+":"+datatype)
+                        datatypesdescr.add(datagroup+":"+datatype_aux)
                     else:
-                        datatypesdescr.add(datagroup + ":" + datatype +
+                        datatypesdescr.add(datagroup + ":" + datatype_aux +
                                            "," + dataset_save + "," +
                                            product_save)
             else:
@@ -540,7 +540,7 @@ def _get_datatype_list(cfg, radarnr='RADAR001'):
                         if ((dataset_save is None) and (product_save is None)):
                             datatypesdescr.add(datagroup+":"+datatype_aux)
                         else:
-                            datatypesdescr.add(datagroup + ":" + datatype +
+                            datatypesdescr.add(datagroup + ":" + datatype_aux +
                                                "," + dataset_save + "," +
                                                product_save)
 
@@ -605,8 +605,7 @@ def _get_masterfile_list(datatypesdescr, starttime, endtime, datacfg,
         radarnr, datagroup, datatype, dataset, product = get_datatype_fields(
             datatypedescr)
         if ((datagroup != 'COSMO') and (datagroup != 'RAD4ALPCOSMO') and
-                (datagroup != 'DEM') and (datagroup != 'RAD4ALPDEM') and
-                (datagroup != 'SAVED')):
+                (datagroup != 'DEM') and (datagroup != 'RAD4ALPDEM')):
             masterdatatypedescr = datatypedescr
             if scan_list is not None:
                 masterscan = scan_list[int(radarnr[5:8])-1][0]
@@ -617,7 +616,7 @@ def _get_masterfile_list(datatypesdescr, starttime, endtime, datacfg,
         for datatypedescr in datatypesdescr:
             radarnr, datagroup, datatype, dataset, product = (
                 get_datatype_fields(datatypedescr))
-            if ((datagroup == 'COSMO') or (datagroup == 'SAVED')):
+            if datagroup == 'COSMO':
                 masterdatatypedescr = radarnr+':RAINBOW:dBZ'
                 if scan_list is not None:
                     masterscan = scan_list[int(radarnr[5:8])-1][0]

src/pyrad_proc/pyrad/graph/plots.py

@@ -475,7 +475,7 @@ def plot_density(hist_obj, hist_type, field_name, ind_sweep, prdcfg,
 
 
 def plot_scatter(bins1, bins2, hist_2d, field_name1, field_name2, fname_list,
-                 prdcfg, metadata=None):
+                 prdcfg, metadata=None,  lin_regr=None, lin_regr_slope1=None):
     """
     2D histogram
 
@@ -523,6 +523,13 @@ def plot_scatter(bins1, bins2, hist_2d, field_name1, field_name2, fname_list,
 
     # plot reference
     plt.plot(bins1, bins2, 'k--')
+
+    # plot linear regression
+    if lin_regr is not None:
+        plt.plot(bins1, lin_regr[0]*bins1+lin_regr[1], 'r')
+    if lin_regr_slope1 is not None:
+        plt.plot(bins1, bins1+lin_regr_slope1, 'g')
+
     plt.autoscale(enable=True, axis='both', tight=True)
 
     plt.xlabel(labelx)

src/pyrad_proc/pyrad/io/__init__.py

@@ -53,6 +53,7 @@
 
     write_colocated_gates
     write_colocated_data
+    write_colocated_data_time_avg
     write_timeseries
     write_ts_polar_data
     write_monitoring_ts
@@ -109,6 +110,7 @@
 from .write_data import write_monitoring_ts
 from .write_data import write_sun_hits, write_sun_retrieval
 from .write_data import write_colocated_gates, write_colocated_data
+from .write_data import write_colocated_data_time_avg
 
 from .io_aux import get_save_dir, make_filename
 from .io_aux import get_datetime, get_dataset_fields

src/pyrad_proc/pyrad/io/read_data_other.py

@@ -12,6 +12,7 @@
     read_rad4alp_vis
     read_colocated_gates
     read_colocated_data
+    read_colocated_data_time_avg
     read_timeseries
     read_monitoring_ts
     read_sun_hits_multiple_days
@@ -287,6 +288,71 @@ def read_colocated_data(fname):
         return None, None, None, None, None, None, None, None
 
 
+def read_colocated_data_time_avg(fname):
+    """
+    Reads a csv files containing time averaged colocated data
+
+    Parameters
+    ----------
+    fname : str
+        path of time series file
+
+    Returns
+    -------
+    rad1_ele , rad1_azi, rad1_rng, rad1_val, rad2_ele, rad2_azi, rad2_rng,
+    rad2_val : tupple
+        A tupple with the data read. None otherwise
+
+    """
+    try:
+        with open(fname, 'r', newline='') as csvfile:
+            # first count the lines
+            reader = csv.DictReader(
+                row for row in csvfile if not row.startswith('#'))
+            nrows = sum(1 for row in reader)
+            rad1_ele = np.empty(nrows, dtype=float)
+            rad1_azi = np.empty(nrows, dtype=float)
+            rad1_rng = np.empty(nrows, dtype=float)
+            rad1_dBZavg = np.empty(nrows, dtype=float)
+            rad1_PhiDPavg = np.empty(nrows, dtype=float)
+            rad1_Flagavg = np.empty(nrows, dtype=float)
+            rad2_ele = np.empty(nrows, dtype=float)
+            rad2_azi = np.empty(nrows, dtype=float)
+            rad2_rng = np.empty(nrows, dtype=float)
+            rad2_dBZavg = np.empty(nrows, dtype=float)
+            rad2_PhiDPavg = np.empty(nrows, dtype=float)
+            rad2_Flagavg = np.empty(nrows, dtype=float)
+
+            # now read the data
+            csvfile.seek(0)
+            reader = csv.DictReader(
+                row for row in csvfile if not row.startswith('#'))
+            i = 0
+            for row in reader:
+                rad1_ele[i] = float(row['rad1_ele'])
+                rad1_azi[i] = float(row['rad1_azi'])
+                rad1_rng[i] = float(row['rad1_rng'])
+                rad1_dBZavg[i] = float(row['rad1_dBZavg'])
+                rad1_PhiDPavg[i] = float(row['rad1_PhiDPavg'])
+                rad1_Flagavg[i] = float(row['rad1_Flagavg'])
+                rad2_ele[i] = float(row['rad2_ele'])
+                rad2_azi[i] = float(row['rad2_azi'])
+                rad2_rng[i] = float(row['rad2_rng'])
+                rad2_dBZavg[i] = float(row['rad2_dBZavg'])
+                rad2_PhiDPavg[i] = float(row['rad2_PhiDPavg'])
+                rad2_Flagavg[i] = float(row['rad2_Flagavg'])
+                i += 1
+
+            return (rad1_ele, rad1_azi, rad1_rng,
+                    rad1_dBZavg, rad1_PhiDPavg, rad1_Flagavg,
+                    rad2_ele, rad2_azi, rad2_rng,
+                    rad2_dBZavg, rad2_PhiDPavg, rad2_Flagavg)
+    except EnvironmentError:
+        warn('Unable to read file '+fname)
+        return (None, None, None, None, None, None, None, None, None, None,
+                None, None)
+
+
 def read_timeseries(fname):
     """
     Reads a time series contained in a csv file

src/pyrad_proc/pyrad/io/write_data.py

@@ -11,6 +11,7 @@
     write_monitoring_ts
     write_colocated_gates
     write_colocated_data
+    write_colocated_data_time_avg
     write_sun_hits
     write_sun_retrieval
     generate_field_name_str
@@ -204,7 +205,7 @@ def write_colocated_gates(coloc_gates, fname):
 
 def write_colocated_data(coloc_data, fname):
     """
-    Writes the position of gates colocated with two radars
+    Writes the data of gates colocated with two radars
 
     Parameters
     ----------
@@ -264,6 +265,82 @@ def write_colocated_data(coloc_data, fname):
     return fname
 
 
+def write_colocated_data_time_avg(coloc_data, fname):
+    """
+    Writes the time averaged data of gates colocated with two radars
+
+    Parameters
+    ----------
+    coloc_data : dict
+        dictionary containing the colocated data parameters
+    fname : str
+        file name where to store the data
+
+    Returns
+    -------
+    fname : str
+        the name of the file where data has written
+
+    """
+    filelist = glob.glob(fname)
+    ngates = len(coloc_data['rad1_ele'])
+    if len(filelist) == 0:
+        with open(fname, 'w', newline='') as csvfile:
+            csvfile.write('# Colocated radar gates data file\n')
+            csvfile.write('# Comment lines are preceded by "#"\n')
+            csvfile.write('#\n')
+
+            fieldnames = [
+                'rad1_ele', 'rad1_azi', 'rad1_rng',
+                'rad1_dBZavg', 'rad1_PhiDPavg', 'rad1_Flagavg',
+                'rad2_ele', 'rad2_azi', 'rad2_rng',
+                'rad2_dBZavg', 'rad2_PhiDPavg', 'rad2_Flagavg']
+            writer = csv.DictWriter(csvfile, fieldnames)
+            writer.writeheader()
+            for i in range(ngates):
+                writer.writerow(
+                    {'rad1_ele': coloc_data['rad1_ele'][i],
+                     'rad1_azi': coloc_data['rad1_azi'][i],
+                     'rad1_rng': coloc_data['rad1_rng'][i],
+                     'rad1_dBZavg': coloc_data['rad1_dBZavg'][i],
+                     'rad1_PhiDPavg': coloc_data['rad1_PhiDPavg'][i],
+                     'rad1_dBZavg': coloc_data['rad1_dBZavg'][i],
+                     'rad1_Flagavg': coloc_data['rad1_Flagavg'][i],
+                     'rad2_ele': coloc_data['rad2_ele'][i],
+                     'rad2_azi': coloc_data['rad2_azi'][i],
+                     'rad2_rng': coloc_data['rad2_rng'][i],
+                     'rad2_dBZavg': coloc_data['rad2_dBZavg'][i],
+                     'rad2_PhiDPavg': coloc_data['rad2_PhiDPavg'][i],
+                     'rad2_Flagavg': coloc_data['rad2_Flagavg'][i]})
+            csvfile.close()
+    else:
+        with open(fname, 'a', newline='') as csvfile:
+            fieldnames = [
+                'rad1_ele', 'rad1_azi', 'rad1_rng',
+                'rad1_dBZavg', 'rad1_PhiDPavg', 'rad1_Flagavg',
+                'rad2_ele', 'rad2_azi', 'rad2_rng',
+                'rad2_dBZavg', 'rad2_PhiDPavg', 'rad2_Flagavg']
+            writer = csv.DictWriter(csvfile, fieldnames)
+            for i in range(ngates):
+                writer.writerow(
+                    {'rad1_ele': coloc_data['rad1_ele'][i],
+                     'rad1_azi': coloc_data['rad1_azi'][i],
+                     'rad1_rng': coloc_data['rad1_rng'][i],
+                     'rad1_dBZavg': coloc_data['rad1_dBZavg'][i],
+                     'rad1_PhiDPavg': coloc_data['rad1_PhiDPavg'][i],
+                     'rad1_dBZavg': coloc_data['rad1_dBZavg'][i],
+                     'rad1_Flagavg': coloc_data['rad1_Flagavg'][i],
+                     'rad2_ele': coloc_data['rad2_ele'][i],
+                     'rad2_azi': coloc_data['rad2_azi'][i],
+                     'rad2_rng': coloc_data['rad2_rng'][i],
+                     'rad2_dBZavg': coloc_data['rad2_dBZavg'][i],
+                     'rad2_PhiDPavg': coloc_data['rad2_PhiDPavg'][i],
+                     'rad2_Flagavg': coloc_data['rad2_Flagavg'][i]})
+            csvfile.close()
+
+    return fname
+
+
 def write_sun_hits(sun_hits, fname):
     """
     Writes sun hits data.

src/pyrad_proc/pyrad/proc/__init__.py

@@ -68,6 +68,7 @@
     process_time_avg_flag
     process_colocated_gates
     process_intercomp
+    process_intercomp_time_avg
 
 Retrievals
 ==========
@@ -108,6 +109,7 @@
 from .process_calib import process_time_avg, process_weighted_time_avg
 from .process_calib import process_time_avg_flag
 from .process_calib import process_colocated_gates, process_intercomp
+from .process_calib import process_intercomp_time_avg
 from .process_calib import process_sun_hits
 
 from .process_retrieve import process_signal_power, process_snr