Feature 1844 python embedding

met/Make-include

@@ -25,6 +25,7 @@ MET_CPPFLAGS = -I${top_builddir}/src/basic/vx_cal \
 	     -I${top_builddir}/src/libcode/vx_nc_util \
 	     -I${top_builddir}/src/libcode/vx_pb_util \
 	     -I${top_builddir}/src/libcode/vx_plot_util \
+	     -I${top_builddir}/src/libcode/vx_pointdata_python \
 	     -I${top_builddir}/src/libcode/vx_ps \
 	     -I${top_builddir}/src/libcode/vx_pxm \
 	     -I${top_builddir}/src/libcode/vx_render \
@@ -67,6 +68,7 @@ MET_LDFLAGS = -L${top_builddir}/src/basic/vx_cal \
 	     -L${top_builddir}/src/libcode/vx_nc_util \
 	     -L${top_builddir}/src/libcode/vx_pb_util \
 	     -L${top_builddir}/src/libcode/vx_plot_util \
+	     -L${top_builddir}/src/libcode/vx_pointdata_python \
 	     -L${top_builddir}/src/libcode/vx_ps \
 	     -L${top_builddir}/src/libcode/vx_pxm \
 	     -L${top_builddir}/src/libcode/vx_render \

met/configure.ac

@@ -1036,7 +1036,7 @@ if test "x$ENABLE_PYTHON" = "xyes"; then
     AC_DEFINE([ENABLE_PYTHON], [], ["build python embedding"])
     AC_MSG_NOTICE([python embedding will be compiled])
     CPPFLAGS="${CPPFLAGS} -DWITH_PYTHON"
-    PYTHON_LIBS="-lvx_data2d_python -lvx_python3_utils ${MET_PYTHON_LD}"
+    PYTHON_LIBS="-lvx_data2d_python -lvx_pointdata_python -lvx_python3_utils ${MET_PYTHON_LD}"
 else
     AC_MSG_NOTICE([python embedding will not be compiled])
     PYTHON_LIBS=
@@ -1233,6 +1233,7 @@ AC_CONFIG_FILES([Makefile
                  src/libcode/vx_summary/Makefile
                  src/libcode/vx_python3_utils/Makefile
                  src/libcode/vx_data2d_python/Makefile
+                 src/libcode/vx_pointdata_python/Makefile
                  src/tools/Makefile
                  src/tools/core/Makefile
                  src/tools/core/ensemble_stat/Makefile

met/docs/Users_Guide/appendixF.rst

@@ -257,3 +257,53 @@ The **read_ascii_mpr.py** sample script can be found in:
 • MET installation directory in *MET_BASE/python*.
 
 • `MET GitHub repository <https://github.com/dtcenter/MET>`_ in *met/scripts/python*.
+
+
+Python Embedding for Point Observations as input
+________________________________________________
+
+
+The point2grid, plot_point_obs, ensemble_stat, and point_stat tools use MET point observation NetCDF. They support the python embedding by the prefix 'PYTHON_NUMPY=" and followed by a python script name instead of the MET point observastion NetCDF filename. The customized python script is expected to extend MET_BASE/python/met_point_obs.py and to produce the python variable, **met_point_data**,  which is the dictionary of the MET point observation data. They are defined at MET_BASE/python/met_point_obs.py.
+
+
+.. _pyembed-point-obs-data:
+
+
+.. code-block:: none
+
+  met_point_data = {
+
+     'use_var_id':  Trur/False,     # obs_vid are variable index if True, otherwise GRIB codes
+
+     # Header data
+     'nhdr':        integer_value,  # number of headers
+     'pbhdr':       integer_value,  # number of PREPBUFR specific headers
+     'nhdr_typ':    integer_value,  # number of message types
+     'nhdr_sid':    integer_value,  # number of station IDs
+     'nhdr_vld':    integer_value,  # number of valid times
+     'hdr_typ':     nympy_integer_array,    # index of message type
+     'hdr_sid':     nympy_integer_array,    # index of station ID
+     'hdr_vld':     nympy_integer_array,    # index of valid time
+     'hdr_lat':     nympy_float_array,      # latitude
+     'hdr_lon':     nympy_float_array,      # longitude
+     'hdr_elv':     nympy_float_array,      # station elevation
+     'hdr_typ_table':   string_value,       # message types
+     'hdr_sid_table':   string_value,       # station IDs
+     'hdr_vld_table':   string_value,       # valid times "yyyymmdd_hhmmss"
+     'hdr_prpt_typ':    nympy_integer_array,   # optional
+     'hdr_irpt_typ':    nympy_integer_array,   # optional
+     'hdr_inst_typ':    nympy_integer_array,   # optional
+
+     # Observation data
+     'nobs':       integer_value,       # number of observation
+     'nobs_qty':   integer_value        # number of quality marks
+     'nobs_var':   integer_value        # number of variable names
+     'obs_qty':    nympy_integer_array, # index of quality mark
+     'obs_hid':    nympy_integer_array, # index of header
+     'obs_vid':    nympy_integer_array, # index of veriable or GRIB code
+     'obs_lvl':    nympy_float_array,   # pressure level
+     'obs_hgt':    nympy_float_array,   # height of observation data
+     'obs_val'     nympy_float_array,   # observatin value
+     'obs_qty_table':  string_array,    # quality marks
+     'obs_var_table':  string_array,    # variable names
+  }

met/docs/Users_Guide/ensemble-stat.rst

@@ -93,7 +93,7 @@ Optional arguments for ensemble_stat
 
 4. To produce ensemble statistics using gridded observations, use the **-grid_obs file** option to specify a gridded observation file. This option may be used multiple times if your observations are in several files.
 
-5. To produce ensemble statistics using point observations, use the **-point_obs file** option to specify a NetCDF point observation file. This option may be used multiple times if your observations are in several files.
+5. To produce ensemble statistics using point observations, use the **-point_obs file** option to specify a NetCDF point observation file. This option may be used multiple times if your observations are in several files. The python embedding will be activated if the **file** begines with 'PYTHON_NUMPY=" and followed by a python script name.
 
 6. To override the simple ensemble mean value of the input ensemble members for the ECNT, SSVAR, and ORANK line types, the **-ens_mean file** option specifies an ensemble mean model data file. This option replaces the **-ssvar_mean file** option from earlier versions of MET.
 

met/docs/Users_Guide/plotting.rst

@@ -8,6 +8,9 @@ __________________
 
 This section describes how to check your data files using plotting utilities. Point observations can be plotted using the Plot-Point-Obs utility. A single model level can be plotted using the plot_data_plane utility. For object based evaluations, the MODE objects can be plotted using plot_mode_field. Occasionally, a post-processing or timing error can lead to errors in MET. These tools can assist the user by showing the data to be verified to ensure that times and locations match up as expected. 
 
+MET version 10.1 adds support for the passing point observations to plot_point_obs using a Python scriptAn example of running plot_point_obs with Python embedding is included below.
+
+
 plot_point_obs usage
 ~~~~~~~~~~~~~~~~~~~~
 
@@ -66,6 +69,18 @@ An example of the plot_point_obs calling sequence is shown below:
 
 In this example, the Plot-Point-Obs tool will process the input sample_pb.nc file and write a postscript file containing a plot to a file named sample_pb.ps.
 
+This is an equivalent command with the python embedding. This is an example for the python embeddingt.
+
+.. code-block:: none
+
+  plot_point_obs 'PYTHON_NUMPY=MET_BASE/python/read_met_point_obs.py sample_pb.nc' sample_data.ps
+
+
+The user should replace the python script with the customized python script for the custom point observation data.
+
+Please refer to :numref:`Appendix F, Section %s <appendixF>` for more details about Python embedding in MET.
+
+
 plot_point_obs configuration file
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

met/docs/Users_Guide/point-stat.rst

@@ -286,7 +286,7 @@ Required arguments for point_stat
 Optional arguments for point_stat
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-4. The **-point_obs** file may be used to pass additional NetCDF point observation files to be used in the verification. 
+4. The **-point_obs** file may be used to pass additional NetCDF point observation files to be used in the verification. The python embedding will be activated if the **file** begines with 'PYTHON_NUMPY=" and followed by a python script name.
 
 5. The **-obs_valid_beg** time option in YYYYMMDD[_HH[MMSS]] format sets the beginning of the observation matching time window, overriding the configuration file setting.
 

met/docs/Users_Guide/reformat_point.rst

@@ -981,6 +981,8 @@ _______________
 
 The Point2Grid tool takes point observations from a NetCDF output file from one of the four previously mentioned MET tools (ascii2nc, madis2nc, pb2nc, lidar2nc) and creates a gridded NetCDF file. The other point observations are GOES-16/17 input files in NetCDF format (especially, Aerosol Optical Depth. Future development will include support for reading input files not produced from MET tools.
 
+MET version 10.1 adds support for the passing point observations to point2grid using a Python script with "input_filename" which begins with the "PYTHON_NUMPY=<python_script_name> [arguments]" instead of MET point observation NetCDF input. An example of running point2grid with Python embedding is included below.
+
 
 point2grid usage
 ~~~~~~~~~~~~~~~~
@@ -1013,6 +1015,8 @@ Required arguments for point2grid
 
 1. The **input_filename** argument indicates the name of the input NetCDF file to be processed. Currently, only NetCDF files produced from the ascii2nc, madis2nc, pb2nc, and lidar2nc are supported. And AOD dataset from GOES16/17 are supported, too. Support for additional file types will be added in future releases.
 
+The MET point observation NetCDF file name as **input_filename** argument is equivalent with "PYTHON_NUMPY=MET_BASE/python/read_met_point_obs.py netcdf_file name'.
+
 2. The **to_grid** argument defines the output grid as: (1) a named grid, (2) the path to a gridded data file, or (3) an explicit grid specification string.
 
 3. The **output_filename** argument is the name of the output NetCDF file to be written.
@@ -1065,6 +1069,21 @@ For the GOES-16 and GOES-17 data, the computing lat/long is time consuming. So t
 When processing GOES-16 data, the **-qc** option may also be used to specify the acceptable quality control flag values. The example above regrids the GOES-16 AOD values to NCEP Grid number 212 (which QC flags are high, medium, and low), writing to the output the maximum AOD value falling inside each grid box.
 
 
+Here is an example of processing the same set of observations but using Python embedding instead:
+
+
+.. code-block:: none
+
+		point2grid \
+		'PYTHON_NUMPY=MET_BASE/python/read_met_point_obs.py ascii2nc_edr_hourly.20130827.nc' \
+		G212 python_gridded_ascii_python.nc -config Point2GridConfig_edr \
+		-field 'name="200"; level="*"; valid_time="20130827_205959";' -method MAX -v 1
+
+The user should replace the python script with the customized python script for the custom point observation data. This is an example for the python embedding.
+
+Please refer to :numref:`Appendix F, Section %s <appendixF>` for more details about Python embedding in MET.
+
+
 point2grid output
 ~~~~~~~~~~~~~~~~~
 

met/scripts/python/Makefile.am

@@ -26,11 +26,13 @@
 pythonscriptsdir = $(pkgdatadir)/python
 
 pythonscripts_DATA = \
+	met_point_obs.py \
 	read_ascii_numpy.py \
 	read_ascii_numpy_grid.py \
 	read_ascii_xarray.py \
 	read_ascii_point.py \
-	read_ascii_mpr.py
+	read_ascii_mpr.py \
+	read_met_point_obs.py
 
 EXTRA_DIST = ${pythonscripts_DATA}
 

met/scripts/python/met_point_obs.py

@@ -0,0 +1,169 @@
+'''
+Created on Nov 10, 2021
+
+@author: hsoh
+'''
+
+
+import numpy as np
+
+COUNT_SHOW = 30
+
+MET_PYTHON_OBS_ARGS = "MET_POINT_PYTHON_ARGS"
+
+class met_point_obs(object):
+    '''
+    classdocs
+    '''
+    python_prefix = 'PYTHON_POINT_RAW'
+
+    def __init__(self, use_var_id=True):
+        '''
+        Constructor
+        '''
+        self.raw_data = None
+        self.use_var_id = use_var_id    # True if variable index, False if GRIB code
+
+        # Header
+        self.nhdr = 0
+        self.npbhdr = 0
+        self.nhdr_typ = 0   # type table
+        self.nhdr_sid = 0   #station_uid table
+        self.nhdr_vld = 0   # valid time strings
+        self.hdr_typ = []   # (nhdr)
+        self.hdr_sid = []   # (nhdr)
+        self.hdr_vld = []   # (nhdr)
+        self.hdr_lat = []   # (nhdr)
+        self.hdr_lon = []   # (nhdr)
+        self.hdr_elv = []   # (nhdr)
+        self.hdr_typ_table = [] # (nhdr_typ, mxstr2)
+        self.hdr_sid_table = [] # (nhdr_sid, mxstr2)
+        self.hdr_vld_table = [] # (nhdr_vld, mxstr)
+        self.hdr_prpt_typ  = [] # optional
+        self.hdr_irpt_typ  = [] # optional
+        self.hdr_inst_typ  = [] # optional
+
+        #Observation data        
+        self.nobs = 0
+        self.nobs_qty = 0
+        self.nobs_var = 0
+        self.obs_qty = []  # (nobs_qty )
+        self.obs_hid = []  # (nobs)
+        self.obs_vid = []  # (nobs) veriable index or GRIB code
+        self.obs_lvl = []  # (nobs)
+        self.obs_hgt = []  # (nobs)
+        self.obs_val = []  # (nobs)
+        self.obs_qty_table = []  # (nobs_qty, mxstr)
+        self.obs_var_table = []  # (nobs_var, mxstr2) optional if GRIB code at self.obs_vid
+
+    #def convert(self):
+    #    pass
+
+    def get_point_data(self):
+        if self.nhdr <= 0:
+            self.nhdr = len(self.hdr_lat)
+        if self.nobs <= 0:
+            self.nobs = len(self.obs_val)
+        if self.nhdr_typ <= 0:
+            self.nhdr_typ = len(self.hdr_typ_table)
+        if self.nhdr_sid <= 0:
+            self.nhdr_sid = len(self.hdr_sid_table)
+        if self.nhdr_vld <= 0:
+            self.nhdr_vld = len(self.hdr_vld_table)
+        if self.npbhdr <= 0:
+            self.npbhdr = len(self.hdr_prpt_typ)
+        if self.nobs_qty <= 0:
+            self.nobs_qty = len(self.obs_qty_table)
+        if self.nobs_var <= 0:
+            self.nobs_var = len(self.obs_var_table)
+        return self.__dict__
+
+    @staticmethod
+    def is_python_script(arg_value):
+        return arg_value.startswith(met_point_obs.python_prefix)
+
+    @staticmethod
+    def get_python_script(arg_value):
+        return arg_value[len(met_point_obs.python_prefix)+1:]
+
+    @staticmethod
+    def print_data(key, data_array, show_count=COUNT_SHOW):
+        if isinstance(data_array, list):
+            data_len = len(data_array)
+            if show_count >= data_len:
+                print("      {k:10s}: {v}".format(k=key, v= data_array))
+            else:
+                end_offset = int(show_count/2)
+                print("      {k:10s}: count={v}".format(k=key, v=data_len))
+                print("      {k:10s}[0:{o}] {v}".format(k=key, v=data_array[:end_offset], o=end_offset))
+                print("      {k:10s}[{s}:{e}]: {v}".format(k=key, v='...', s=end_offset+1, e=data_len-end_offset-1))
+                print("      {k:10s}[{s}:{e}]: {v}".format(k=key, v= data_array[-end_offset:], s=(data_len-end_offset), e=(data_len-1)))
+        else:
+            print("      {k:10s}: {v}".format(k=key, v= data_array))
+
+    @staticmethod
+    def print_point_data(met_point_data):
+        print(' === MET point data by python embedding ===')
+        met_point_obs.print_data('nhdr', met_point_data['nhdr'])
+        met_point_obs.print_data('nobs' , met_point_data['nobs'])
+        met_point_obs.print_data('use_var_id',met_point_data['use_var_id'])
+        met_point_obs.print_data('hdr_typ',met_point_data['hdr_typ'])
+        met_point_obs.print_data('obs_hid',met_point_data['obs_hid'])
+        met_point_obs.print_data('obs_vid',met_point_data['obs_vid'])
+        met_point_obs.print_data('obs_val',met_point_data['obs_val'])
+        met_point_obs.print_data('obs_var_table',met_point_data['obs_var_table'])
+        met_point_obs.print_data('obs_qty_table',met_point_data['obs_qty_table'])
+        met_point_obs.print_data('hdr_typ_table',met_point_data['hdr_typ_table'])
+        met_point_obs.print_data('hdr_sid_table',met_point_data['hdr_sid_table'])
+        met_point_obs.print_data('hdr_vld_table',met_point_data['hdr_vld_table'])
+        print(' === MET point data by python embedding ===')
+
+
+def get_sample_point_obs():
+    point_obs_data = met_point_obs()
+    point_obs_data.hdr_typ = np.array([ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ])
+    point_obs_data.hdr_sid = np.array([ 0, 0, 0, 0, 0, 1, 2, 3, 3, 1, 2, 2, 3, 0, 0, 0, 0, 0, 1, 2, 3, 3, 1, 2, 2, 3 ])
+    point_obs_data.hdr_vld = np.array([ 0, 1, 2, 3, 4, 4, 3, 4, 3, 4, 5, 4, 3, 0, 1, 2, 3, 4, 4, 3, 4, 3, 4, 5, 4, 3 ])
+    point_obs_data.hdr_lat = np.array([ 43., 43., 43., 43., 43., 43., 43., 43., 43., 46., 46., 46., 46., 43., 43., 43., 43., 43., 43., 43., 43., 43., 46., 46., 46., 46. ])
+    point_obs_data.hdr_lon = np.array([ -89., -89., -89., -89., -89., -89., -89., -89., -89., -92., -92., -92., -92., -89., -89., -89., -89., -89., -89., -89., -89., -89., -92., -92., -92., -92. ])
+    point_obs_data.hdr_elv = np.array([ 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220., 220. ])
+
+    point_obs_data.obs_qid = np.array([ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ])
+    point_obs_data.obs_hid = np.array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 25 ])
+    point_obs_data.obs_vid = np.array([ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ])
+    point_obs_data.obs_lvl = np.array([ 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000. ])
+    point_obs_data.obs_hgt = np.array([ 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2. ])
+    point_obs_data.obs_val = np.array([ 292., 292.5, 293., 293.5, 294., 294.5, 295., 295.5, 296., 292., 293.4, 293., 296., 294., 92., 92.5, 93., 93.5, 94., 94.5, 95., 95.5, 96., 92., 93.4, 93., 96., 94. ])
+
+    point_obs_data.hdr_typ_table = [ "ADPSFC" ]
+    point_obs_data.hdr_sid_table = [ "001", "002", "003", "004" ]
+    point_obs_data.hdr_vld_table = [
+            "20120409_115000", "20120409_115500", "20120409_120100", "20120409_120500", "20120409_121000",
+            "20120409_120000" ]
+    point_obs_data.obs_var_table = [ "TMP", "RH" ]
+    point_obs_data.obs_qty_table = [ "NA" ]
+    point_obs_data.nhdr = len(point_obs_data.hdr_lat)
+    point_obs_data.nobs = len(point_obs_data.obs_val)
+    #point_obs_data.met_point_data = point_obs_data
+    return point_obs_data
+
+
+def main():
+    met_point_data = get_sample_point_obs().get_point_data()
+
+    print('All',met_point_data)
+    print("      nhdr: ",met_point_data['nhdr'])
+    print("      nobs: ",met_point_data['nobs'])
+    print('use_var_id: ',met_point_data['use_var_id'])
+    print('hdr_typ: ',met_point_data['hdr_typ'])
+    print('obs_vid: ',met_point_data['obs_vid'])
+    print('obs_val: ',met_point_data['obs_val'])
+    print('obs_var_table: ',met_point_data['obs_var_table'])
+    print('obs_qty_table: ',met_point_data['obs_qty_table'])
+    print('hdr_typ_table: ',met_point_data['hdr_typ_table'])
+    print('hdr_sid_table: ',met_point_data['hdr_sid_table'])
+    print('hdr_vld_table: ',met_point_data['hdr_vld_table'])
+    print('Done python scripot')
+
+if __name__ == '__main__':
+    main()