preparation for implementation of Giangrande melting layer detection …

…algorithm

config/pyart/mch_config.py

@@ -148,6 +148,7 @@
 radar_estimated_rain_rate = 'radar_estimated_rain_rate'
 radar_echo_classification = 'radar_echo_classification'
 radar_echo_id = 'radar_echo_id'
+melting_layer = 'melting_layer'
 
 # attenuation
 specific_attenuation = 'specific_attenuation'
@@ -319,6 +320,7 @@
     'radar_estimated_rain_rate': radar_estimated_rain_rate,
     'radar_echo_classification': radar_echo_classification,
     'radar_echo_id': radar_echo_id,
+    'melting_layer': melting_layer,
     'specific_attenuation': specific_attenuation,
     'path_integrated_attenuation': path_integrated_attenuation,
     'specific_differential_attenuation': specific_differential_attenuation,
@@ -993,6 +995,15 @@
         'boundaries': [0.5, 1.5, 2.5, 3.5],
         'coordinates': 'elevation azimuth range'},
 
+    melting_layer: {
+        'units': 'legend',
+        'standard_name': 'melting_layer',
+        'long_name': 'Position of the range bin respect to the melting layer',
+        'labels': ['BELOW', 'ENTERING', 'INSIDE', 'EXITING', 'ABOVE'],
+        'ticks': [1, 2, 3, 4, 5],
+        'boundaries': [0.5, 1.5, 2.5, 3.5, 4.5, 5.5],
+        'coordinates': 'elevation azimuth range'},
+
     specific_attenuation: {
         'units': 'dB/km',
         'standard_name': 'specific_attenuation',
@@ -1948,6 +1959,7 @@ def spectrum_width_limit(container=None, selection=0):
 
     radar_echo_classification: 'pyart_LangRainbow12',
     radar_echo_id: 'pyart_LangRainbow12',
+    melting_layer: 'pyart_LangRainbow12',
 
     specific_attenuation: 'pyart_Carbone17',
     path_integrated_attenuation: 'pyart_Carbone17',
@@ -2065,8 +2077,9 @@ def spectrum_width_limit(container=None, selection=0):
     rain_rate: (0., 10.),
     radar_estimated_rain_rate: (0., 10.),
 
-    radar_echo_classification: (-0.5, 9.5),
+    radar_echo_classification: (0., 9.),
     radar_echo_id: (0, 3),
+    melting_layer: (0, 5),
 
     sun_hit_h: (0, 1),
     sun_hit_v: (0, 1),

src/pyrad_proc/pyrad/proc/__init__.py

@@ -34,6 +34,7 @@
     process_filter_visibility
     process_outlier_filter
     process_hydroclass
+    process_melting_layer
 
 Phase processing and attenuation correction
 ===========================================
@@ -124,7 +125,7 @@
 from .process_echoclass import process_echo_id, process_echo_filter
 from .process_echoclass import process_filter_snr, process_filter_visibility
 from .process_echoclass import process_outlier_filter, process_hydroclass
-from .process_echoclass import process_cdf
+from .process_echoclass import process_cdf, process_melting_layer
 
 from .process_phase import process_correct_phidp0
 from .process_phase import process_smooth_phidp_single_window

src/pyrad_proc/pyrad/proc/process_aux.py

@@ -118,6 +118,8 @@ def get_process_func(dataset_type, dsname):
         func_name = 'process_windshear'
     elif dataset_type == 'HYDROCLASS':
         func_name = 'process_hydroclass'
+    elif dataset_type == 'ML_DETECTION':
+        func_name = 'process_melting_layer'
     elif dataset_type == 'PHIDP0_ESTIMATE':
         func_name = 'process_estimate_phidp0'
     elif dataset_type == 'RHOHV_RAIN':

src/pyrad_proc/pyrad/proc/process_echoclass.py

@@ -14,6 +14,7 @@
     process_filter_visibility
     process_outlier_filter
     process_hydroclass
+    process_melting_layer
 
 """
 
@@ -806,3 +807,98 @@ def process_hydroclass(procstatus, dscfg, radar_list=None):
         new_dataset.add_field('radar_echo_classification', hydro)
 
         return new_dataset, ind_rad
+
+
+def process_melting_layer(procstatus, dscfg, radar_list=None):
+    """
+    Detects the melting layer
+
+    Parameters
+    ----------
+    procstatus : int
+        Processing status: 0 initializing, 1 processing volume,
+        2 post-processing
+    dscfg : dictionary of dictionaries
+        data set configuration. Accepted Configuration Keywords::
+
+        datatype : list of string. Dataset keyword
+            The input data types
+    radar_list : list of Radar objects
+        Optional. list of radar objects
+
+    Returns
+    -------
+    new_dataset : Radar
+        radar object
+    ind_rad : int
+        radar index
+
+    """
+    if procstatus != 1:
+        return None, None
+
+    if dscfg['ML_METHOD'] == 'GIANGRANDE':
+        temp_field = None
+        iso0_field = None
+        for datatypedescr in dscfg['datatype']:
+            radarnr, datagroup, datatype, dataset, product = (
+                get_datatype_fields(datatypedescr))
+            if datatype == 'dBZ':
+                refl_field = 'reflectivity'
+            if datatype == 'dBZc':
+                refl_field = 'corrected_reflectivity'
+            if datatype == 'ZDR':
+                zdr_field = 'differential_reflectivity'
+            if datatype == 'ZDRc':
+                zdr_field = 'corrected_differential_reflectivity'
+            if datatype == 'RhoHV':
+                rhv_field = 'cross_correlation_ratio'
+            if datatype == 'RhoHVc':
+                rhv_field = 'corrected_cross_correlation_ratio'
+            if datatype == 'TEMP':
+                temp_field = 'temperature'
+            if datatype == 'H_ISO0':
+                iso0_field = 'height_over_iso0'
+
+        ind_rad = int(radarnr[5:8])-1
+        if radar_list[ind_rad] is None:
+            warn('No valid radar')
+            return None, None
+        radar = radar_list[ind_rad]
+
+        if temp_field is None and iso0_field is None:
+            warn('iso0 or temperature fields needed to detect melting layer')
+            return None, None
+
+        if temp_field is not None and (temp_field not in radar.fields):
+            warn('Unable to detect melting layer. Missing temperature field')
+            return None, None
+
+        if iso0_field is not None and (iso0_field not in radar.fields):
+            warn('Unable to detect melting layer. ' +
+                 'Missing height over iso0 field')
+            return None, None
+
+        temp_ref = 'temperature'
+        if iso0_field is not None:
+            temp_ref = 'height_over_iso0'
+
+        if ((refl_field not in radar.fields) or
+                (zdr_field not in radar.fields) or
+                (rhv_field not in radar.fields)):
+            warn('Unable to detect melting layer. Missing data')
+            return None, None
+
+        # User defined variables here
+
+        ml = pyart.retrieve.melting_layer_giangrande(
+            radar, refl_field=refl_field, zdr_field=zdr_field,
+            rhv_field=rhv_field, temp_field=temp_field, iso0_field=iso0_field,
+            ml_field=None, temp_ref=temp_ref)
+
+        # prepare for exit
+        new_dataset = deepcopy(radar)
+        new_dataset.fields = dict()
+        new_dataset.add_field('melting_layer', ml)
+
+        return new_dataset, ind_rad