diff --git a/met/Make-include b/met/Make-include index dfa6bd5faa..42f1d5f3b3 100644 --- a/met/Make-include +++ b/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 \ diff --git a/met/configure.ac b/met/configure.ac index 7008748079..df6a93a428 100644 --- a/met/configure.ac +++ b/met/configure.ac @@ -1043,7 +1043,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= @@ -1240,6 +1240,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 diff --git a/met/docs/Users_Guide/appendixF.rst b/met/docs/Users_Guide/appendixF.rst index 3a4046e9ea..586c264c9e 100644 --- a/met/docs/Users_Guide/appendixF.rst +++ b/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 `_ 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 + } diff --git a/met/docs/Users_Guide/ensemble-stat.rst b/met/docs/Users_Guide/ensemble-stat.rst index 5c6277b179..fd1dbf535a 100644 --- a/met/docs/Users_Guide/ensemble-stat.rst +++ b/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. diff --git a/met/docs/Users_Guide/plotting.rst b/met/docs/Users_Guide/plotting.rst index 7f748e7885..a96dca605a 100644 --- a/met/docs/Users_Guide/plotting.rst +++ b/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 ` for more details about Python embedding in MET. + + plot_point_obs configuration file ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/met/docs/Users_Guide/point-stat.rst b/met/docs/Users_Guide/point-stat.rst index 7a9966511c..ce7b74e905 100644 --- a/met/docs/Users_Guide/point-stat.rst +++ b/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. diff --git a/met/docs/Users_Guide/reformat_point.rst b/met/docs/Users_Guide/reformat_point.rst index 999393b72d..9e41ca80dc 100644 --- a/met/docs/Users_Guide/reformat_point.rst +++ b/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= [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 ` for more details about Python embedding in MET. + + point2grid output ~~~~~~~~~~~~~~~~~ diff --git a/met/scripts/python/Makefile.am b/met/scripts/python/Makefile.am index 78dc7d88bc..368da3da7f 100644 --- a/met/scripts/python/Makefile.am +++ b/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} diff --git a/met/scripts/python/met_point_obs.py b/met/scripts/python/met_point_obs.py new file mode 100755 index 0000000000..c306839226 --- /dev/null +++ b/met/scripts/python/met_point_obs.py @@ -0,0 +1,189 @@ +''' +Created on Nov 10, 2021 + +@author: hsoh + +- This is the base class and the customized script should extend the met_point_obs. +- The customized script (for example "custom_reader") must implement "def read_data(self, args)" + which fills the array (list) variables at __init__(). +- The args can be 1) single string argument, 2) the list of arguments, and 3) the dictionary of arguments. +- The variable "met_point_data" must be set for MET tools +- The customized script is expected to include following codes: + + # prepare arguments for the customized script + args = {'input', sys.argv[1]} # or args = [] + point_obs_data = custom_reader() + point_obs_data.read_data() + met_point_data = point_obs_data.get_point_data() +''' + +from abc import ABC, abstractmethod +import numpy as np + +COUNT_SHOW = 30 + +MET_PYTHON_OBS_ARGS = "MET_POINT_PYTHON_ARGS" + +class met_point_obs(ABC): + ''' + classdocs + ''' + + python_prefix = 'PYTHON_POINT_RAW' + + def __init__(self, use_var_id=True): + ''' + Constructor + ''' + 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 + + @abstractmethod + def read_data(self, args): + # args should be list or dictionary + 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_subset=True): + print(' === MET point data by python embedding ===') + if print_subset: + 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']) + else: + 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(' === MET point data by python embedding ===') + + +# This is a sample drived class +class sample_met_point_obs(met_point_obs): + + #@abstractmethod + def read_data(self, arg_map={}): + self.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 ]) + self.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 ]) + self.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 ]) + self.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. ]) + self.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. ]) + self.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. ]) + + self.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 ]) + self.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 ]) + self.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 ]) + self.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. ]) + self.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. ]) + self.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. ]) + + self.hdr_typ_table = [ "ADPSFC" ] + self.hdr_sid_table = [ "001", "002", "003", "004" ] + self.hdr_vld_table = [ + "20120409_115000", "20120409_115500", "20120409_120100", "20120409_120500", "20120409_121000", + "20120409_120000" ] + self.obs_var_table = [ "TMP", "RH" ] + self.obs_qty_table = [ "NA" ] + + +def main(): + args = {} # or args = [] + point_obs_data = sample_met_point_obs() + point_obs_data.read_data(args) + met_point_data = point_obs_data.get_point_data() + + point_obs_data.print_point_data(met_point_data, print_subset=False) + +if __name__ == '__main__': + main() + print('Done python scripot') diff --git a/met/scripts/python/read_met_point_obs.py b/met/scripts/python/read_met_point_obs.py new file mode 100755 index 0000000000..15667541ec --- /dev/null +++ b/met/scripts/python/read_met_point_obs.py @@ -0,0 +1,125 @@ +''' +Created on Nov 10, 2021 + +@author: hsoh + +This script reads the MET point observation NetCDF file like MET tools do. +''' + +import os +import sys +import time +import numpy as np +import netCDF4 as nc + +from met_point_obs import met_point_obs, sample_met_point_obs + +DO_PRINT_DATA = False +ARG_PRINT_DATA = 'show_data' + +class nc_point_obs(met_point_obs): + + def read_data(self, args): + # args should be list or dictionaryu + if isinstance(args, dict): + nc_filename = args.get('nc_name',None) + elif isinstance(args, list): + nc_filename = args[0] + else: + nc_filename = args + if nc_filename is None: + print("==ERROR== The input NetCDF filename is missing") + elif not os.path.exists(nc_filename): + print("==ERROR== input NetCDF file ({f}) does not exist".format(f=nc_filename)) + else: + dataset = nc.Dataset(nc_filename, 'r') + # Header + self.hdr_typ = dataset['hdr_typ'][:] + self.hdr_sid = dataset['hdr_sid'][:] + self.hdr_vld = dataset['hdr_vld'][:] + self.hdr_lat = dataset['hdr_lat'][:] + self.hdr_lon = dataset['hdr_lon'][:] + self.hdr_elv = dataset['hdr_elv'][:] + self.hdr_typ_table = nc_tools.get_string_array(dataset, 'hdr_typ_table') + self.hdr_sid_table = nc_tools.get_string_array(dataset, 'hdr_sid_table') + self.hdr_vld_table = nc_tools.get_string_array(dataset, 'hdr_vld_table') + + nc_var = dataset.variables.get('hdr_prpt_typ', None) + if nc_var: + self.hdr_prpt_typ = nc_var[:] + nc_var = dataset.variables.get('hdr_irpt_typ', None) + if nc_var: + self.hdr_irpt_typ = nc_var[:] + nc_var = dataset.variables.get('hdr_inst_typ', None) + if nc_var: + self.hdr_inst_typ =nc_var[:] + + #Observation data + self.hdr_sid = dataset['hdr_sid'][:] + self.obs_qty = np.array(dataset['obs_qty'][:]) + self.obs_hid = np.array(dataset['obs_hid'][:]) + self.obs_lvl = np.array(dataset['obs_lvl'][:]) + self.obs_hgt = np.array(dataset['obs_hgt'][:]) + self.obs_val = np.array(dataset['obs_val'][:]) + nc_var = dataset.variables.get('obs_vid', None) + if nc_var is None: + self.use_var_id = False + nc_var = dataset.variables.get('obs_gc', None) + else: + self.obs_var_table = nc_tools.get_string_array(dataset, 'obs_var') + if nc_var: + self.obs_vid = np.array(nc_var[:]) + + self.obs_qty_table = nc_tools.get_string_array(dataset, 'obs_qty_table') + + +class nc_tools(): + + @staticmethod + def get_num_array(dataset, var_name): + nc_var = dataset.variables.get(var_name, None) + return nc_var[:].filled(nc_var._FillValue if '_FillValue' in nc_var.ncattrs() else -9999) if nc_var else [] + + @staticmethod + def get_ncbyte_array_to_str(nc_var): + nc_str_data = nc_var[:] + if nc_var.datatype.name == 'bytes8': + nc_str_data = [ str(s.compressed(),"utf-8") for s in nc_var[:]] + return nc_str_data + + @staticmethod + def get_string_array(dataset, var_name): + nc_var = dataset.variables.get(var_name, None) + return nc_tools.get_ncbyte_array_to_str(nc_var) if nc_var else [] + + +perf_start_time = time.time() +perf_start_counter = time.perf_counter_ns() + +point_obs_data = None +if len(sys.argv) == 1: + point_obs_data = sample_met_point_obs() + point_obs_data.read_data([]) +else: + netcdf_filename = sys.argv[1] + args = [ netcdf_filename ] + #args = { 'nc_name': netcdf_filename } + point_obs_data = nc_point_obs() + point_obs_data.read_data(args) + + if ARG_PRINT_DATA == sys.argv[-1]: + DO_PRINT_DATA = True + +if point_obs_data is not None: + met_point_data = point_obs_data.get_point_data() + met_point_data['met_point_data'] = point_obs_data + + if DO_PRINT_DATA: + met_point_obs.print_point_data(met_point_data) + +perf_end_time = time.time() +perf_end_counter = time.perf_counter_ns() +perf_duration = perf_end_time - perf_start_time +perf_duration_counter = (perf_end_counter - perf_start_counter) / 1000000000 + +print('Done python script {s} Took walltime: {t1} & perf: {t2} seconds'.format(s=sys.argv[0], t1=perf_duration, t2=perf_duration_counter)) diff --git a/met/src/libcode/Makefile.am b/met/src/libcode/Makefile.am index 76817f7749..8ed5d82c73 100644 --- a/met/src/libcode/Makefile.am +++ b/met/src/libcode/Makefile.am @@ -27,6 +27,7 @@ SUBDIRS = vx_grid \ if ENABLE_PYTHON SUBDIRS += vx_python3_utils SUBDIRS += vx_data2d_python + SUBDIRS += vx_pointdata_python endif diff --git a/met/src/libcode/vx_nc_obs/Makefile.am b/met/src/libcode/vx_nc_obs/Makefile.am index ba2fb8c910..22f0871fab 100644 --- a/met/src/libcode/vx_nc_obs/Makefile.am +++ b/met/src/libcode/vx_nc_obs/Makefile.am @@ -12,6 +12,7 @@ include ${top_srcdir}/Make-include noinst_LIBRARIES = libvx_nc_obs.a libvx_nc_obs_a_SOURCES = \ + met_point_data.cc met_point_data.h \ nc_obs_util.cc nc_obs_util.h \ nc_point_obs.cc nc_point_obs.h \ nc_point_obs_in.cc nc_point_obs_in.h \ diff --git a/met/src/libcode/vx_nc_obs/met_point_data.cc b/met/src/libcode/vx_nc_obs/met_point_data.cc new file mode 100644 index 0000000000..6cc15a7e81 --- /dev/null +++ b/met/src/libcode/vx_nc_obs/met_point_data.cc @@ -0,0 +1,437 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + + +//////////////////////////////////////////////////////////////////////// + + +using namespace std; + +#include +#include +#include +#include +#include +#include +#include + +#include "vx_log.h" +#include "is_bad_data.h" + +#include "met_point_data.h" +//#include "nc_point_obs.h" + +//////////////////////////////////////////////////////////////////////// + + + // + // Code for class MetPointData + // + + +//////////////////////////////////////////////////////////////////////// + +MetPointData::MetPointData() { + // Derived class should set obs_data + init_from_scratch(); +} + +//////////////////////////////////////////////////////////////////////// + +MetPointData::~MetPointData() { + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +void MetPointData::init_from_scratch() { + nobs = 0; + nhdr = 0; + qty_length = 0; + + use_var_id = false; + use_arr_vars = false; +} + + +//////////////////////////////////////////////////////////////////////// + +//void MetPointData::allocate() { +// obs_data.allocate(); +//} + +//////////////////////////////////////////////////////////////////////// + +void MetPointData::clear() { + if (obs_data) obs_data->clear(); + header_data.clear(); +} + +//////////////////////////////////////////////////////////////////////// + +bool MetPointData::get_header(int header_offset, float hdr_arr[HDR_ARRAY_LEN], + ConcatString &hdr_typ_str, ConcatString &hdr_sid_str, + ConcatString &hdr_vld_str) { + int hdr_idx; + + // Read the corresponding header array for this observation + + hdr_arr[0] = header_data.lat_array[header_offset]; + hdr_arr[1] = header_data.lon_array[header_offset]; + hdr_arr[2] = header_data.elv_array[header_offset]; + + // Read the corresponding header type for this observation + hdr_idx = use_arr_vars ? header_offset : header_data.typ_idx_array[header_offset]; + hdr_typ_str = header_data.typ_array[hdr_idx]; + + // Read the corresponding header Station ID for this observation + hdr_idx = use_arr_vars ? header_offset : header_data.sid_idx_array[header_offset]; + hdr_sid_str = header_data.sid_array[hdr_idx]; + + // Read the corresponding valid time for this observation + hdr_idx = use_arr_vars ? header_offset : header_data.vld_idx_array[header_offset]; + hdr_vld_str = header_data.vld_array[hdr_idx]; + + return true; +} + +//////////////////////////////////////////////////////////////////////// + +bool MetPointData::get_header_type(int header_offset, int hdr_typ_arr[HDR_TYPE_ARR_LEN]) { + int hdr_idx; + // Read the header integer types + hdr_typ_arr[0] = (header_data.prpt_typ_array.n() > header_offset ? + header_data.prpt_typ_array[header_offset] : bad_data_int); + hdr_typ_arr[1] = (header_data.irpt_typ_array.n() > header_offset ? + header_data.irpt_typ_array[header_offset] : bad_data_int); + hdr_typ_arr[2] = (header_data.inst_typ_array.n() > header_offset ? + header_data.inst_typ_array[header_offset] : bad_data_int); + + return true; +} + +//////////////////////////////////////////////////////////////////////// + +bool MetPointData::get_lats(float *hdr_lats) { + for (int idx=0; idxobs_cnt = obs_cnt; +} + + + + +//////////////////////////////////////////////////////////////////////// + + + // + // Code for class MetPointDataPython + // + + +//////////////////////////////////////////////////////////////////////// + +MetPointDataPython::MetPointDataPython() { + obs_data = new MetPointObsData(); + init_from_scratch(); +} + +//////////////////////////////////////////////////////////////////////// + +MetPointDataPython::MetPointDataPython(MetPointDataPython &d) { + init_from_scratch(); + //obs_data = d.get_point_obs_data(); + //header_data = d.get_header_data(); + obs_data->assign(*d.get_point_obs_data()); + header_data.assign(*d.get_header_data()); +cout << " DEBUG HS MetPointData(MetPointData &d) is called \n"; +cout << " DEBUG HS MetPointData(MetPointData &d) &header_data.lat_array=" << &(header_data.lat_array) << "\n"; +cout << " DEBUG HS MetPointData(MetPointData &d) header_data.lat_array.n()=" << header_data.lat_array.n() << "\n"; +cout << " DEBUG HS MetPointData(MetPointData &d) header_data.lon_array.n()=" << header_data.lon_array.n() << "\n"; +cout << " DEBUG HS MetPointData(MetPointData &d) header_data.elv_array.n()=" << header_data.elv_array.n() << "\n"; +cout << " DEBUG HS MetPointData(MetPointData &d) header_data.typ_idx_array.n()=" << header_data.typ_idx_array.n() << "\n"; +cout << " DEBUG HS MetPointData(MetPointData &d) header_data.sid_idx_array.n()=" << header_data.sid_idx_array.n() << "\n"; +cout << " DEBUG HS MetPointData(MetPointData &d) header_data.vld_idx_array.n()=" << header_data.vld_idx_array.n() << "\n"; +} + + +//////////////////////////////////////////////////////////////////////// + +MetPointDataPython::~MetPointDataPython() { + clear(); +} + +//////////////////////////////////////////////////////////////////////// + +//void MetPointDataPython::init_from_scratch() { +// nobs = 0; +// nhdr = 0; +// qty_length = 0; +// +// use_var_id = false; +// use_arr_vars = false; +//} + + +//////////////////////////////////////////////////////////////////////// + +void MetPointDataPython::allocate(int obs_cnt) { + set_obs_cnt(obs_cnt); + obs_data->allocate(); +} + +//////////////////////////////////////////////////////////////////////// + + + +/////////////////////////////////////////////////////////////////////////////// +// struct MetPointObsData + +MetPointObsData::MetPointObsData(): + obs_cnt(0), + obs_ids((int *)0), + obs_hids((int *)0), + obs_qids((int *)0), + obs_lvls((float *)0), + obs_hgts((float *)0), + obs_vals((float *)0), + obs_arr((float *)0), + is_obs_array(false) +{ +} + + +/////////////////////////////////////////////////////////////////////////////// + +void MetPointObsData::allocate() { + if (is_obs_array) obs_arr = new float[obs_cnt*OBS_ARRAY_LEN]; // nobs * 5 + else { + obs_ids = new int[obs_cnt]; // grib_code or var_id + obs_hids = new int[obs_cnt]; + obs_qids = new int[obs_cnt]; + obs_lvls = new float[obs_cnt]; + obs_hgts = new float[obs_cnt]; + obs_vals = new float[obs_cnt]; + } +} + +/////////////////////////////////////////////////////////////////////////////// + +void MetPointObsData::assign(MetPointObsData &o) { + obs_cnt = o.obs_cnt; + is_obs_array = o.is_obs_array; + + clear(); + allocate(); + if (is_obs_array) + for (int idx=0; idx + +#include "nc_utils.h" + + +//////////////////////////////////////////////////////////////////////// + +struct MetPointHeader { + int typ_len; + int sid_len; + int vld_len; + int strl_len; + int strll_len; + int min_vld_time; + int max_vld_time; + int hdr_count; + int hdr_type_count; + + StringArray typ_array; + StringArray sid_array; + StringArray vld_array; // string array for valid time + IntArray vld_num_array; // number array for valid time + IntArray typ_idx_array; + IntArray sid_idx_array; + IntArray vld_idx_array; // index to vld_array/vld_num_array + NumArray lat_array; + NumArray lon_array; + NumArray elv_array; + IntArray prpt_typ_array; + IntArray irpt_typ_array; + IntArray inst_typ_array; + + MetPointHeader(); + void assign(MetPointHeader &h); + void clear(); + void reset_counters(); +}; + +//////////////////////////////////////////////////////////////////////// + +struct MetPointObsData { + int obs_cnt; + bool is_obs_array; + + int *obs_ids; // grib_code or var_id + int *obs_hids; + int *obs_qids; + float *obs_lvls; + float *obs_hgts; + float *obs_vals; + float *obs_arr; // nobs * 5 + StringArray var_names; + StringArray qty_names; + + MetPointObsData(); + void allocate(); + void assign(MetPointObsData &d); + void clear(); + void clear_numbers(); + void clear_strings(); + bool fill_obs_buf(int buf_size, int offset, float *obs_arr, int *qty_idx_arr); + float get_obs_val(int index); +}; + + +class MetPointData { + + protected: + + int nobs; + int nhdr; + int qty_length; + bool use_var_id; + bool use_arr_vars; + + MetPointHeader header_data; + MetPointObsData *obs_data = 0; + + void init_from_scratch(); + + public: + + MetPointData(); + ~MetPointData(); + + //void close(); + void clear(); + + int get_buf_size(); + int get_hdr_cnt(); + int get_grib_code_or_var_index(const float obs_arr[OBS_ARRAY_LEN]); + MetPointHeader *get_header_data(); + bool get_header(int header_offset, float hdr_arr[HDR_ARRAY_LEN], + ConcatString &hdr_typ_str, ConcatString &hdr_sid_str, + ConcatString &hdr_vld_str); + int get_header_offset(const float obs_arr[OBS_ARRAY_LEN]); + bool get_header_type(int header_offset, int hdr_typ_arr[HDR_TYPE_ARR_LEN]); + bool get_lats(float *hdr_lats); + bool get_lons(float *hdr_lons); + int get_obs_cnt(); + MetPointObsData *get_point_obs_data(); + StringArray get_qty_data(); + StringArray get_var_names(); + + bool is_same_obs_values(const float obs_arr1[OBS_ARRAY_LEN], const float obs_arr2[OBS_ARRAY_LEN]); + bool is_using_var_id(); + + void set_grib_code_or_var_index(float obs_arr[OBS_ARRAY_LEN], int grib_code); + void set_hdr_cnt(int hdr_cnt); + void set_obs_cnt(int obs_cnt); + // variables + + // data + +}; // MetPointDataBase + + + +class MetPointDataPython : public MetPointData { + + protected: + + public: + + MetPointDataPython(); + MetPointDataPython(MetPointDataPython &d); + ~MetPointDataPython(); + + void allocate(int obs_cnt); + void set_use_var_id(bool new_use_var_id); + + // variables + + // data + +}; // MetPointData + +//////////////////////////////////////////////////////////////////////// + +inline MetPointHeader *MetPointData::get_header_data() { return &header_data; } +inline int MetPointData::get_buf_size() { return OBS_BUFFER_SIZE; } +inline int MetPointData::get_grib_code_or_var_index(const float obs_arr[OBS_ARRAY_LEN]) { return obs_arr[1]; }; +inline int MetPointData::get_hdr_cnt() { return nhdr; } +inline int MetPointData::get_header_offset(const float obs_arr[OBS_ARRAY_LEN]) { return obs_arr[0]; }; +inline int MetPointData::get_obs_cnt() { return nobs; } +inline MetPointObsData *MetPointData::get_point_obs_data() { return obs_data; } +inline StringArray MetPointData::get_qty_data() { return obs_data->qty_names; } +inline StringArray MetPointData::get_var_names() { return obs_data->var_names; } +inline bool MetPointData::is_using_var_id() { return use_var_id; } +inline void MetPointData::set_grib_code_or_var_index(float obs_arr[OBS_ARRAY_LEN], int grib_code) { obs_arr[1] = grib_code; } +inline void MetPointDataPython::set_use_var_id(bool new_use_var_id) { use_var_id = new_use_var_id; } + +//////////////////////////////////////////////////////////////////////// + + +#endif /* __MET_POINT_DATA_H__ */ + + +//////////////////////////////////////////////////////////////////////// + diff --git a/met/src/libcode/vx_nc_obs/nc_obs_util.cc b/met/src/libcode/vx_nc_obs/nc_obs_util.cc index 70ea844d49..39a59266d2 100644 --- a/met/src/libcode/vx_nc_obs/nc_obs_util.cc +++ b/met/src/libcode/vx_nc_obs/nc_obs_util.cc @@ -68,78 +68,12 @@ void NcDataBuffer::reset_counters() { /////////////////////////////////////////////////////////////////////////////// // struct NcPointObsData -NcPointObsData::NcPointObsData(): - obs_cnt(0), - obs_ids((int *)0), - obs_hids((int *)0), - obs_qids((int *)0), - obs_lvls((float *)0), - obs_hgts((float *)0), - obs_vals((float *)0), - obs_arr((float *)0), - is_obs_array(false) +NcPointObsData::NcPointObsData() { } /////////////////////////////////////////////////////////////////////////////// -void NcPointObsData::clear() { - obs_cnt = 0; - is_obs_array = false; - - clear_numbers(); - clear_strings(); -} - -/////////////////////////////////////////////////////////////////////////////// - -void NcPointObsData::clear_numbers() { - if (0 != obs_ids) { - delete [] obs_ids; - obs_ids = (int *)0; - } - if (0 != obs_hids) { - delete [] obs_hids; - obs_hids = (int *)0; - } - if (0 != obs_qids) { - delete [] obs_qids; - obs_qids = (int *)0; - } - if (0 != obs_lvls) { - delete [] obs_lvls; - obs_lvls = (float *)0; - } - if (0 != obs_hgts) { - delete [] obs_hgts; - obs_hgts = (float *)0; - } - if (0 != obs_vals) { - delete [] obs_vals; - obs_vals = (float *)0; - } - if (0 != obs_arr) { - delete [] obs_arr; - obs_arr = (float *)0; - } -} - -/////////////////////////////////////////////////////////////////////////////// - -void NcPointObsData::clear_strings() { - var_names.clear(); - qty_names.clear(); -} - -/////////////////////////////////////////////////////////////////////////////// - -float NcPointObsData::get_obs_val(int index) { - float obs_val = (is_obs_array ? obs_arr[index*obs_cnt+4] : obs_vals[index]); - return obs_val; -} - -/////////////////////////////////////////////////////////////////////////////// - bool NcPointObsData::read_obs_data_numbers(NetcdfObsVars obs_vars, bool stop) { bool succeed = true; const char* method_name = "NcPointObsData::read_obs_data_numbers()"; @@ -434,7 +368,7 @@ void NetcdfObsVars::create_pb_hdrs (NcFile *f_out, const int hdr_count) { /////////////////////////////////////////////////////////////////////////////// -void NetcdfObsVars::create_table_vars (NcFile *f_out, NcHeaderData &hdr_data, +void NetcdfObsVars::create_table_vars (NcFile *f_out, MetPointHeader &hdr_data, NcDataBuffer &data_buffer) { const string method_name = " create_table_vars()"; @@ -603,7 +537,7 @@ void NetcdfObsVars::read_dims_vars(NcFile *f_in) { //////////////////////////////////////////////////////////////////////// -void NetcdfObsVars::read_header_data(NcHeaderData &hdr_data) { +void NetcdfObsVars::read_header_data(MetPointHeader &hdr_data) { bool is_valid_obs_nc = true; long nhdr_count = get_dim_size(&hdr_dim); int strl_len = get_dim_size(&strl_dim); @@ -970,7 +904,7 @@ bool NetcdfObsVars::read_obs_data(int buf_size, int offset, //////////////////////////////////////////////////////////////////////// -void NetcdfObsVars::read_pb_hdr_data(NcHeaderData &hdr_data) { +void NetcdfObsVars::read_pb_hdr_data(MetPointHeader &hdr_data) { const char *method_name = "get_nc_pb_hdr_data() -> "; int pb_hdr_count = get_dim_size(&pb_hdr_dim); @@ -1166,7 +1100,7 @@ void NetcdfObsVars::write_header_to_nc(NcDataBuffer &data_buf, /////////////////////////////////////////////////////////////////////////////// -void NetcdfObsVars::write_table_vars (NcHeaderData &hdr_data, NcDataBuffer &data_buffer) +void NetcdfObsVars::write_table_vars (MetPointHeader &hdr_data, NcDataBuffer &data_buffer) { mlog << Debug(5) << "write_table_vars() is called. valid hdr_typ_tbl_var: " << !IS_INVALID_NC(hdr_typ_tbl_var) << "\n"; @@ -1261,51 +1195,6 @@ void NetcdfObsVars::write_obs_var_descriptions(StringArray &descriptions) { //////////////////////////////////////////////////////////////////////// // End of NetcdfObsVars -/////////////////////////////////////////////////////////////////////////////// - -// struct NcHeaderData - -NcHeaderData::NcHeaderData() -{ - reset_counters(); -} - -/////////////////////////////////////////////////////////////////////////////// - -void NcHeaderData::clear() { - reset_counters(); - - typ_array.clear(); - sid_array.clear(); - vld_array.clear(); - vld_num_array.clear(); - typ_idx_array.clear(); - sid_idx_array.clear(); - vld_idx_array.clear(); - lat_array.clear(); - lon_array.clear(); - elv_array.clear(); - prpt_typ_array.clear(); - irpt_typ_array.clear(); - inst_typ_array.clear(); -} - -/////////////////////////////////////////////////////////////////////////////// - -void NcHeaderData::reset_counters() { - valid_point_obs = false; - typ_len = 0; - sid_len = 0; - vld_len = 0; - strl_len = 0; - strll_len = 0; - hdr_count = 0; - - min_vld_time = -1; - max_vld_time = -1; -} - - /////////////////////////////////////////////////////////////////////////////// long count_nc_headers(vector< Observation > &observations) diff --git a/met/src/libcode/vx_nc_obs/nc_obs_util.h b/met/src/libcode/vx_nc_obs/nc_obs_util.h index 4985639c09..d3d07f3007 100644 --- a/met/src/libcode/vx_nc_obs/nc_obs_util.h +++ b/met/src/libcode/vx_nc_obs/nc_obs_util.h @@ -17,6 +17,7 @@ using namespace netCDF; #include "nc_summary.h" +#include "met_point_data.h" //////////////////////////////////////////////////////////////////////// @@ -78,39 +79,6 @@ struct NcDataBuffer { //////////////////////////////////////////////////////////////////////// -struct NcHeaderData { - bool valid_point_obs; - int typ_len; - int sid_len; - int vld_len; - int strl_len; - int strll_len; - int min_vld_time; - int max_vld_time; - int hdr_count; - int hdr_type_count; - - StringArray typ_array; - StringArray sid_array; - StringArray vld_array; - IntArray vld_num_array; - IntArray typ_idx_array; - IntArray sid_idx_array; - IntArray vld_idx_array; - NumArray lat_array; - NumArray lon_array; - NumArray elv_array; - IntArray prpt_typ_array; - IntArray irpt_typ_array; - IntArray inst_typ_array; - - NcHeaderData(); - void clear(); - void reset_counters(); -}; - -//////////////////////////////////////////////////////////////////////// - struct NetcdfObsVars { bool attr_agl ; bool attr_pb2nc ; @@ -167,7 +135,7 @@ struct NetcdfObsVars { void create_hdr_vars (NcFile *f_out, const int hdr_count); void create_obs_vars (NcFile *f_out); void create_obs_name_vars (NcFile *f_out, const int var_count, const int unit_count); - void create_table_vars (NcFile *f_out, NcHeaderData &hdr_data, NcDataBuffer &data_buffer); + void create_table_vars (NcFile *f_out, MetPointHeader &hdr_data, NcDataBuffer &data_buffer); void create_pb_hdrs (NcFile *f_out, const int hdr_count); NcDim create_var_obs_var (NcFile *f_out, int var_count); @@ -175,41 +143,25 @@ struct NetcdfObsVars { int get_obs_index(); void read_dims_vars(NcFile *f_in); - void read_header_data(NcHeaderData &hdr_data); + void read_header_data(MetPointHeader &hdr_data); bool read_obs_data(int buf_size, int offset, int qty_len, float *obs_arr, int *qty_idx_arr, char *obs_qty_buf); - void read_pb_hdr_data(NcHeaderData &hdr_data); + void read_pb_hdr_data(MetPointHeader &hdr_data); void write_header_to_nc(NcDataBuffer &data_buf, const int buf_size, const bool is_pb = false); void write_obs_buffer(NcDataBuffer &data_buffer, const int buf_size); void write_obs_var_names(StringArray &obs_names); void write_obs_var_units(StringArray &units); void write_obs_var_descriptions(StringArray &descriptions); - void write_table_vars(NcHeaderData &hdr_data, NcDataBuffer &data_buffer); + void write_table_vars(MetPointHeader &hdr_data, NcDataBuffer &data_buffer); }; // NetcdfObsVars //////////////////////////////////////////////////////////////////////// -struct NcPointObsData { - int obs_cnt; - bool is_obs_array; - - int *obs_ids; // grib_code or var_id - int *obs_hids; - int *obs_qids; - float *obs_lvls; - float *obs_hgts; - float *obs_vals; - float *obs_arr; // nobs * 5 - StringArray var_names; - StringArray qty_names; +struct NcPointObsData : public MetPointObsData { NcPointObsData(); - void clear(); - void clear_numbers(); - void clear_strings(); - float get_obs_val(int index); bool read_obs_data_numbers(NetcdfObsVars obs_vars, bool stop=true); bool read_obs_data_table_lookups(NetcdfObsVars obs_vars, bool stop=true); }; diff --git a/met/src/libcode/vx_nc_obs/nc_point_obs.cc b/met/src/libcode/vx_nc_obs/nc_point_obs.cc index b665e9d2e9..0869f3a31e 100644 --- a/met/src/libcode/vx_nc_obs/nc_point_obs.cc +++ b/met/src/libcode/vx_nc_obs/nc_point_obs.cc @@ -36,6 +36,7 @@ using namespace std; //////////////////////////////////////////////////////////////////////// MetNcPointObs::MetNcPointObs() { + obs_data = new NcPointObsData(); init_from_scratch(); } @@ -48,104 +49,21 @@ MetNcPointObs::~MetNcPointObs() { //////////////////////////////////////////////////////////////////////// void MetNcPointObs::init_from_scratch() { - obs_nc = (NcFile *) 0; + MetPointData::init_from_scratch(); - nobs = 0; - nhdr = 0; - qty_length = 0; keep_nc = false; - use_var_id = false; - use_arr_vars = false; - - return; + obs_nc = (NcFile *) 0; } //////////////////////////////////////////////////////////////////////// void MetNcPointObs::close() { + MetPointData::clear(); + if ( !keep_nc && obs_nc ) { delete obs_nc; obs_nc = (NcFile *) 0; } - - obs_data.clear(); - header_data.clear(); - return; -} - -//////////////////////////////////////////////////////////////////////// - -int MetNcPointObs::get_qty_length() { - qty_length = get_nc_string_length(&obs_vars.obs_qty_tbl_var); - return qty_length; -} - -//////////////////////////////////////////////////////////////////////// - -bool MetNcPointObs::get_header(int header_offset, float hdr_arr[HDR_ARRAY_LEN], - ConcatString &hdr_typ_str, ConcatString &hdr_sid_str, ConcatString &hdr_vld_str) { - int hdr_idx; - - // Read the corresponding header array for this observation - hdr_arr[0] = header_data.lat_array[header_offset]; - hdr_arr[1] = header_data.lon_array[header_offset]; - hdr_arr[2] = header_data.elv_array[header_offset]; - - // Read the corresponding header type for this observation - hdr_idx = use_arr_vars ? header_offset : header_data.typ_idx_array[header_offset]; - hdr_typ_str = header_data.typ_array[hdr_idx]; - - // Read the corresponding header Station ID for this observation - hdr_idx = use_arr_vars ? header_offset : header_data.sid_idx_array[header_offset]; - hdr_sid_str = header_data.sid_array[hdr_idx]; - - // Read the corresponding valid time for this observation - hdr_idx = use_arr_vars ? header_offset : header_data.vld_idx_array[header_offset]; - hdr_vld_str = header_data.vld_array[hdr_idx]; - - return true; -} - -//////////////////////////////////////////////////////////////////////// - -bool MetNcPointObs::get_header_type(int header_offset, int hdr_typ_arr[HDR_TYPE_ARR_LEN]) { - int hdr_idx; - // Read the header integer types - hdr_typ_arr[0] = (header_data.prpt_typ_array.n() > header_offset ? - header_data.prpt_typ_array[header_offset] : bad_data_int); - hdr_typ_arr[1] = (header_data.irpt_typ_array.n() > header_offset ? - header_data.irpt_typ_array[header_offset] : bad_data_int); - hdr_typ_arr[2] = (header_data.inst_typ_array.n() > header_offset ? - header_data.inst_typ_array[header_offset] : bad_data_int); - - return true; -} - -//////////////////////////////////////////////////////////////////////// - -bool MetNcPointObs::get_lats(float *hdr_lats) { - for (int idx=0; idxread_obs_data_numbers(obs_vars); } return status; } @@ -104,7 +104,7 @@ bool MetNcPointObsIn::read_obs_data_numbers() { bool MetNcPointObsIn::read_obs_data_table_lookups() { bool status = false; if( IS_VALID_NC_P(obs_nc) ) { - status = obs_data.read_obs_data_table_lookups(obs_vars); + status = ((NcPointObsData *)obs_data)->read_obs_data_table_lookups(obs_vars); } return status; } diff --git a/met/src/libcode/vx_pointdata_python/Makefile.am b/met/src/libcode/vx_pointdata_python/Makefile.am new file mode 100644 index 0000000000..c48c07e2df --- /dev/null +++ b/met/src/libcode/vx_pointdata_python/Makefile.am @@ -0,0 +1,19 @@ +## @start 1 +## Makefile.am -- Process this file with automake to produce Makefile.in +## @end 1 + +MAINTAINERCLEANFILES = Makefile.in + +# Include the project definitions + +include ${top_srcdir}/Make-include + +# The library + +noinst_LIBRARIES = libvx_pointdata_python.a +libvx_pointdata_python_a_SOURCES = \ + pointdata_python.h pointdata_python.cc \ + pointdata_from_array.h pointdata_from_array.cc \ + python_pointdata.h python_pointdata.cc +libvx_pointdata_python_a_CPPFLAGS = ${MET_CPPFLAGS} -I../vx_python2_utils ${MET_PYTHON_CC} $(MET_PYTHON_LD) + diff --git a/met/src/libcode/vx_pointdata_python/pointdata_from_array.cc b/met/src/libcode/vx_pointdata_python/pointdata_from_array.cc new file mode 100644 index 0000000000..b3d6cd84fe --- /dev/null +++ b/met/src/libcode/vx_pointdata_python/pointdata_from_array.cc @@ -0,0 +1,659 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + + +#include "string.h" + +#include "vx_python3_utils.h" +#include "vx_statistics.h" +#include "check_endian.h" + +#include "pointdata_from_array.h" + +//////////////////////////////////////////////////////////////////////// + +static const int api_delug_level = 11; + +//////////////////////////////////////////////////////////////////////// + + +template +void load_numpy (void * buf, + const int n, + const int data_endian, + void (*shuf)(void *), + float * out) +{ + +bool need_swap = (shuf != 0) && (native_endian != data_endian); + +int j, x, y, r, c; +T * u = (T *) buf; +T value; + +for (j=0; j +void load_numpy (void * buf, + const int n, + const int data_endian, + void (*shuf)(void *), + int * out) +{ + +const char *method_name = "load_numpy(int *) "; +bool need_swap = (shuf != 0) && (native_endian != data_endian); + +int j; +T * u = (T *) buf; +T value; + +for (j=0; j +void load_numpy_int (void * buf, + const int n, + const int data_endian, + void (*shuf)(void *), + IntArray *out) +{ + +const char *method_name = "load_numpy_int(IntArray *) "; +bool need_swap = (shuf != 0) && (native_endian != data_endian); + +int j; +T * u = (T *) buf; +T value; + +out->extend(n); + +for (j=0; jadd((int)value); + + mlog << Debug(api_delug_level) << method_name << " [" << j << "] value=" << value << "\n"; +} // for j + + + +return; + +} + + +//////////////////////////////////////////////////////////////////////// + + +template +void load_numpy_num (void * buf, + const int n, + const int data_endian, + void (*shuf)(void *), + NumArray *out) +{ + +const char *method_name = "load_numpy_num(NumArray *) "; +bool need_swap = (shuf != 0) && (native_endian != data_endian); + +int j; +T * u = (T *) buf; +T value; + +out->extend(n); + +for (j=0; jadd((float)value); + + mlog << Debug(api_delug_level) << method_name << "[" << j << "] value=" << value << "\n"; +} // for j + + + +return; + +} + + +//////////////////////////////////////////////////////////////////////// + + +bool pointdata_from_np_array(Python3_Numpy & np, float * data_out) +{ + +const char *method_name = "pointdata_from_np_array(float) -> "; + + // + // make sure it's a 1D array + // + +if ( np.n_dims() != 1 ) { + + mlog << Error << "\n" << method_name + << "data array is not 1-dimensional! ... " + << "(dim = " << (np.n_dims()) << ")\n\n"; + + exit ( 1 ); + + +} + +int n = np.dim(0); + + + // + // load the data + // + +const ConcatString dtype = np.dtype(); + + // 1 byte integers + + if ( dtype == "|i1" ) load_numpy (np.buffer(), n, little_endian, 0, data_out); +else if ( dtype == "|u1" ) load_numpy (np.buffer(), n, little_endian, 0, data_out); + + // 2 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_2, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_2, data_out); + +else if ( dtype == ">i2" ) load_numpy (np.buffer(), n, big_endian, shuffle_2, data_out); +else if ( dtype == ">u2" ) load_numpy (np.buffer(), n, big_endian, shuffle_2, data_out); + + // 4 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); + +else if ( dtype == ">i4" ) load_numpy (np.buffer(), n, big_endian, shuffle_4, data_out); +else if ( dtype == ">u4" ) load_numpy (np.buffer(), n, big_endian, shuffle_4, data_out); + + // 8 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); + +else if ( dtype == ">i8" ) load_numpy (np.buffer(), n, big_endian, shuffle_8, data_out); +else if ( dtype == ">u8" ) load_numpy (np.buffer(), n, big_endian, shuffle_8, data_out); + + // single precision floats + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); +else if ( dtype == ">f4" ) load_numpy (np.buffer(), n, big_endian, shuffle_4, data_out); + + // double precision floats + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); +else if ( dtype == ">f8" ) load_numpy (np.buffer(), n, big_endian, shuffle_8, data_out); + + // + // nope ... the only other numerical data type for numpy arrays + // is single or double precision complex numbers, and + // we're not supporting those at this time + // + +else { + + mlog << Error << "\n" << method_name + << "unsupported data type \"" << dtype << "\"\n\n"; + + exit ( 1 ); + +} + + //////////////////// + + + // + // done + // + +return ( true ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +bool pointdata_from_np_array(Python3_Numpy & np, int * data_out) +{ + +const char *method_name = "pointdata_from_np_array(int) -> "; + + // + // make sure it's a 1D array + // + +if ( np.n_dims() != 1 ) { + + mlog << Error << "\n" << method_name + << "data array is not 1-dimensional! ... " + << "(dim = " << (np.n_dims()) << ")\n\n"; + + exit ( 1 ); + + +} + +int n = np.dim(0); + + + // + // load the data + // + +const ConcatString dtype = np.dtype(); + + // 1 byte integers + + if ( dtype == "|i1" ) load_numpy (np.buffer(), n, little_endian, 0, data_out); +else if ( dtype == "|u1" ) load_numpy (np.buffer(), n, little_endian, 0, data_out); + + // 2 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_2, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_2, data_out); + +else if ( dtype == ">i2" ) load_numpy (np.buffer(), n, big_endian, shuffle_2, data_out); +else if ( dtype == ">u2" ) load_numpy (np.buffer(), n, big_endian, shuffle_2, data_out); + + // 4 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); + +else if ( dtype == ">i4" ) load_numpy (np.buffer(), n, big_endian, shuffle_4, data_out); +else if ( dtype == ">u4" ) load_numpy (np.buffer(), n, big_endian, shuffle_4, data_out); + + // 8 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); + +else if ( dtype == ">i8" ) load_numpy (np.buffer(), n, big_endian, shuffle_8, data_out); +else if ( dtype == ">u8" ) load_numpy (np.buffer(), n, big_endian, shuffle_8, data_out); + + // single precision floats + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); +else if ( dtype == ">f4" ) load_numpy (np.buffer(), n, big_endian, shuffle_4, data_out); + + // double precision floats + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); +else if ( dtype == ">f8" ) load_numpy (np.buffer(), n, big_endian, shuffle_8, data_out); + + // + // nope ... the only other numerical data type for numpy arrays + // is single or double precision complex numbers, and + // we're not supporting those at this time + // + +else { + + mlog << Error << "\n" << method_name + << "unsupported data type \"" << dtype << "\"\n\n"; + + exit ( 1 ); + +} + + //////////////////// + + + // + // done + // + +return ( true ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +bool pointdata_from_np_array(Python3_Numpy & np, IntArray *data_out) +{ + +const char *method_name = "pointdata_from_np_array(IntArray) -> "; + + // + // make sure it's a 1D array + // + +if ( np.n_dims() != 1 ) { + + mlog << Error << "\n" << method_name + << "data array is not 1-dimensional! ... " + << "(dim = " << (np.n_dims()) << ")\n\n"; + + exit ( 1 ); + + +} + +int n = np.dim(0); + + + // + // load the data + // + +IntArray data; +const ConcatString dtype = np.dtype(); + + // 1 byte integers + + if ( dtype == "|i1" ) load_numpy_int (np.buffer(), n, little_endian, 0, data_out); +else if ( dtype == "|u1" ) load_numpy_int (np.buffer(), n, little_endian, 0, data_out); + + // 2 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_2, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_2, data_out); + +else if ( dtype == ">i2" ) load_numpy_int (np.buffer(), n, big_endian, shuffle_2, data_out); +else if ( dtype == ">u2" ) load_numpy_int (np.buffer(), n, big_endian, shuffle_2, data_out); + + // 4 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); + +else if ( dtype == ">i4" ) load_numpy_int (np.buffer(), n, big_endian, shuffle_4, data_out); +else if ( dtype == ">u4" ) load_numpy_int (np.buffer(), n, big_endian, shuffle_4, data_out); + + // 8 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); + +else if ( dtype == ">i8" ) load_numpy_int (np.buffer(), n, big_endian, shuffle_8, data_out); +else if ( dtype == ">u8" ) load_numpy_int (np.buffer(), n, big_endian, shuffle_8, data_out); + + // single precision floats + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); +else if ( dtype == ">f4" ) load_numpy_int (np.buffer(), n, big_endian, shuffle_4, data_out); + + // double precision floats + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); +else if ( dtype == ">f8" ) load_numpy_int (np.buffer(), n, big_endian, shuffle_8, data_out); + + // + // nope ... the only other numerical data type for numpy arrays + // is single or double precision complex numbers, and + // we're not supporting those at this time + // + +else { + + mlog << Error << "\n" << method_name + << "unsupported data type \"" << dtype << "\"\n\n"; + + exit ( 1 ); + +} + + //////////////////// + + + // + // done + // + +return ( true ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +bool pointdata_from_np_array(Python3_Numpy & np, NumArray *data_out) +{ + +const char *method_name = "pointdata_from_np_array(NumArray) -> "; + + // + // make sure it's a 1D array + // + +if ( np.n_dims() != 1 ) { + + mlog << Error << "\n" << method_name + << "data array is not 1-dimensional! ... " + << "(dim = " << (np.n_dims()) << ")\n\n"; + + exit ( 1 ); + + +} + +int n = np.dim(0); + + + // + // load the data + // + +NumArray data; +const ConcatString dtype = np.dtype(); + + // 1 byte integers + + if ( dtype == "|i1" ) load_numpy_num (np.buffer(), n, little_endian, 0, data_out); +else if ( dtype == "|u1" ) load_numpy_num (np.buffer(), n, little_endian, 0, data_out); + + // 2 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_2, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_2, data_out); + +else if ( dtype == ">i2" ) load_numpy_num (np.buffer(), n, big_endian, shuffle_2, data_out); +else if ( dtype == ">u2" ) load_numpy_num (np.buffer(), n, big_endian, shuffle_2, data_out); + + // 4 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); + +else if ( dtype == ">i4" ) load_numpy_num (np.buffer(), n, big_endian, shuffle_4, data_out); +else if ( dtype == ">u4" ) load_numpy_num (np.buffer(), n, big_endian, shuffle_4, data_out); + + // 8 byte integers + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); + +else if ( dtype == ">i8" ) load_numpy_num (np.buffer(), n, big_endian, shuffle_8, data_out); +else if ( dtype == ">u8" ) load_numpy_num (np.buffer(), n, big_endian, shuffle_8, data_out); + + // single precision floats + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_4, data_out); +else if ( dtype == ">f4" ) load_numpy_num (np.buffer(), n, big_endian, shuffle_4, data_out); + + // double precision floats + +else if ( dtype == " (np.buffer(), n, little_endian, shuffle_8, data_out); +else if ( dtype == ">f8" ) load_numpy_num (np.buffer(), n, big_endian, shuffle_8, data_out); + + // + // nope ... the only other numerical data type for numpy arrays + // is single or double precision complex numbers, and + // we're not supporting those at this time + // + +else { + + mlog << Error << "\n" << method_name + << "unsupported data type \"" << dtype << "\"\n\n"; + + exit ( 1 ); + +} + + //////////////////// + + + // + // done + // + +return ( true ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +bool pointdata_from_str_array(PyObject *data_array, StringArray *data_out) +{ + +//const char *method_name = "pointdata_from_str_array(StringArray) -> "; + +StringArray a = pyobject_as_string_array(data_array); +data_out->clear(); +data_out->add(a); + + //////////////////// + + + // + // done + // + +return ( true ); + +} + + +//////////////////////////////////////////////////////////////////////// + + // + // we just grab the numpy array and the attributes dictionary + // + // from the xarray DataArray object, and then hand them + // + // off to pointdata_from_numpy_array + // + +bool pointdata_from_xarray(PyObject * data_array, float *data_out) +{ + +Python3_Numpy np; +PyObject *numpy_array = PyObject_GetAttrString(data_array, data_attr_name); + + ///////////////////// + +np.set(numpy_array); + +bool status = pointdata_from_np_array(np, data_out); + + // + // done + // + +return ( status ); + +} + + +//////////////////////////////////////////////////////////////////////// + + // + // we just grab the numpy array and the attributes dictionary + // + // from the xarray DataArray object, and then hand them + // + // off to pointdata_from_numpy_array + // + +bool pointdata_from_xarray(PyObject * data_array, int *data_out) +{ + +Python3_Numpy np; +PyObject *numpy_array = PyObject_GetAttrString(data_array, data_attr_name); + + + ///////////////////// + + +np.set(numpy_array); + +bool status = pointdata_from_np_array(np, data_out); + + // + // done + // + +return ( status ); + +} + + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_pointdata_python/pointdata_from_array.h b/met/src/libcode/vx_pointdata_python/pointdata_from_array.h new file mode 100644 index 0000000000..bec9380705 --- /dev/null +++ b/met/src/libcode/vx_pointdata_python/pointdata_from_array.h @@ -0,0 +1,59 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + + +#ifndef __MET_POINTDATA_FROM_ARRAY_H__ +#define __MET_POINTDATA_FROM_ARRAY_H__ + + +//////////////////////////////////////////////////////////////////////// + + +#include "met_point_data.h" + +#include "python3_dict.h" +#include "python3_numpy.h" + + +extern "C" { + +#include "Python.h" + +} + + +//////////////////////////////////////////////////////////////////////// + +static const char data_attr_name [] = "values"; + +//////////////////////////////////////////////////////////////////////// + + +extern bool pointdata_from_np_array(Python3_Numpy & np, int * data_out); +extern bool pointdata_from_np_array(Python3_Numpy & np, float * data_out); +extern bool pointdata_from_np_array(Python3_Numpy & np, IntArray *data_out); +extern bool pointdata_from_np_array(Python3_Numpy & np, NumArray *data_out); +extern bool pointdata_from_str_array(PyObject *data_array, StringArray *data_out); + + +extern bool pointdata_from_xarray(PyObject *data_array, int *data_out); +extern bool pointdata_from_xarray(PyObject *data_array, float *data_out); +extern bool pointdata_from_xarray(PyObject *data_array, IntArray &data_out); +extern bool pointdata_from_xarray(PyObject *data_array, NumArray &data_out); +extern bool pointdata_from_xarray(PyObject *data_array, StringArray &data_out); + +//////////////////////////////////////////////////////////////////////// + + +#endif /* __MET_POINTDATA_FROM_NUMPY_ARRAY_H__ */ + + +//////////////////////////////////////////////////////////////////////// + diff --git a/met/src/libcode/vx_pointdata_python/pointdata_python.cc b/met/src/libcode/vx_pointdata_python/pointdata_python.cc new file mode 100644 index 0000000000..846966a6f1 --- /dev/null +++ b/met/src/libcode/vx_pointdata_python/pointdata_python.cc @@ -0,0 +1,290 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + + +//////////////////////////////////////////////////////////////////////// + + +using namespace std; + +#include +#include +#include +#include + +#include "pointdata_python.h" +#include "pointdata_from_array.h" +#include "vx_python3_utils.h" + +#include "vx_math.h" +#include "vx_log.h" + + +//////////////////////////////////////////////////////////////////////// + + + // + // Code for class MetPythonPointDataFile + // + + +//////////////////////////////////////////////////////////////////////// + + +MetPythonPointDataFile::MetPythonPointDataFile() + +{ + +python_init_from_scratch(); + +} + + +//////////////////////////////////////////////////////////////////////// + + +MetPythonPointDataFile::~MetPythonPointDataFile() + +{ + +close(); + +} + + +//////////////////////////////////////////////////////////////////////// + + +MetPythonPointDataFile::MetPythonPointDataFile(const MetPythonPointDataFile &) + +{ + +mlog << Error << "\nMetPythonPointDataFile::MetPythonPointDataFile(const MetPythonPointDataFile &) -> " + << "should never be called!\n\n"; + +exit ( 1 ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +MetPythonPointDataFile & MetPythonPointDataFile::operator=(const MetPythonPointDataFile &) + +{ + +mlog << Error << "\nMetPythonPointDataFile::operator=(const MetPythonPointDataFile &) -> " + << "should never be called!\n\n"; + +exit ( 1 ); + +return ( * this ); + +} + + +//////////////////////////////////////////////////////////////////////// + + +void MetPythonPointDataFile::python_init_from_scratch() + +{ + +PythonCommand.clear(); +//met_data.clear(); + +close(); + +return; + +} + + +//////////////////////////////////////////////////////////////////////// + + +void MetPythonPointDataFile::close() + +{ + +met_data.clear(); + + // + // Don't reset the Type field + // Don't reset the PythonCommand + // + +return; + +} + + +//////////////////////////////////////////////////////////////////////// + +/* +void MetPythonPointDataFile::set_type(const GrdFileType t) + +{ + +Type = t; + +return; + +} +*/ + +//////////////////////////////////////////////////////////////////////// + + +bool MetPythonPointDataFile::open(const char * cur_command, bool use_xarray) + +{ + +close(); + +ConcatString full_path, file_name; +int i, file_argc; +char **file_argv = (char **) 0; // allocated +StringArray sa; +const char *method_name = "MetPythonPointDataFile::open() "; + + // + // Store the PythonCommand that is being run + // + +PythonCommand = cur_command; + + // + // parse and store argc and argv + // + +sa = PythonCommand.split(" "); + +file_argc = sa.n_elements(); + +if ( file_argc > 0 ) { + file_argv = new char * [ file_argc ]; + char a_var_name[512+1]; + + for ( i=0; idump(out, depth + 1); + +} else { + + out << prefix << "No Grid!\n"; + +} +*/ + // + // done + // + +out.flush(); + +return; + +} + + +//////////////////////////////////////////////////////////////////////// +/* + +bool MetPythonPointDataFile::data_ok(int x, int y) const + +{ + +//const double value = get(x, y); + +//return ( !is_bad_data(value) ); +return true; + +} + + +//////////////////////////////////////////////////////////////////////// + + +void MetPythonPointDataFile::data_minmax(double & data_min, double & data_max) const + +{ + +//Plane.data_range(data_min, data_max); + +return; + +} +*/ + +//////////////////////////////////////////////////////////////////////// + + +MetPointDataPython *MetPythonPointDataFile::get_met_point_data() + +{ + +return &met_data; + +} + + +//////////////////////////////////////////////////////////////////////// + diff --git a/met/src/libcode/vx_pointdata_python/pointdata_python.h b/met/src/libcode/vx_pointdata_python/pointdata_python.h new file mode 100644 index 0000000000..22b32c7f52 --- /dev/null +++ b/met/src/libcode/vx_pointdata_python/pointdata_python.h @@ -0,0 +1,112 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + + +//////////////////////////////////////////////////////////////////////// + + +#ifndef __MET_VX_POINTDATA_PYTHON_H__ +#define __MET_VX_POINTDATA_PYTHON_H__ + + +//////////////////////////////////////////////////////////////////////// + + +#include "data_plane.h" +#include "data_class.h" +#include "two_to_one.h" + +#include "python_pointdata.h" // the main access point +#include "pointdata_from_array.h" // takes an NumPy array or xarray DataArray + +#include "global_python.h" + + +//////////////////////////////////////////////////////////////////////// + + +//class MetPythonPointDataFile : public Met2dDataFile { +class MetPythonPointDataFile { + + private: + + void python_init_from_scratch(); + + MetPythonPointDataFile(const MetPythonPointDataFile &); + MetPythonPointDataFile & operator=(const MetPythonPointDataFile &); + + ConcatString PythonCommand; + + MetPointDataPython met_data; + + //VarInfoPython VInfo; + + //GrdFileType Type; // FileType_Python_Xarray or FileType_Python_Numpy + + public: + + MetPythonPointDataFile(); + ~MetPythonPointDataFile(); + + + // + // set stuff + // + + //void set_type(const GrdFileType); + + // + // get stuff + // + +// GrdFileType file_type() const; + +// double operator () (int x, int y) const; + + //double get (int x, int y) const; + +// bool data_ok (int x, int y) const; + +// void data_minmax (double & data_min, double & data_max) const; + + // + // do stuff + // + + bool open(const char * cur_command, bool use_xarray); + + void close(); + + + void dump (ostream &, int depth = 0) const; + + //MetPointData get_met_point_data(); + MetPointDataPython *get_met_point_data(); + + //int data_plane_array(VarInfo &, DataPlaneArray &); + + //bool met_point_data(MetPointData &); + +}; + + +//////////////////////////////////////////////////////////////////////// + + +//inline double MetPythonPointDataFile::operator () (int x, int y) const { return ( get(x, y) ); } +//inline GrdFileType MetPythonPointDataFile::file_type () const { return ( Type ); } + + +//////////////////////////////////////////////////////////////////////// + + +#endif /* __MET_VX_POINTDATA_PYTHON_H__ */ + + +//////////////////////////////////////////////////////////////////////// + diff --git a/met/src/libcode/vx_pointdata_python/python_pointdata.cc b/met/src/libcode/vx_pointdata_python/python_pointdata.cc new file mode 100644 index 0000000000..43196acc5f --- /dev/null +++ b/met/src/libcode/vx_pointdata_python/python_pointdata.cc @@ -0,0 +1,588 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + + +#include "vx_python3_utils.h" +#include "python_pointdata.h" +#include "pointdata_from_array.h" +#include "vx_util.h" + + +#include "global_python.h" +#include "wchar_argv.h" + +//////////////////////////////////////////////////////////////////////// + +//////////////////////////////////////////////////////////////////////// + + +static bool straight_python_point_data(const char * script_name, + int script_argc, char ** script_argv, + const bool use_xarray, MetPointDataPython &met_pd_out); + +//////////////////////////////////////////////////////////////////////// + +template +static void set_array_from_python(PyObject *python_data, const char *python_key, T *out, bool required=true) { + const char *method_name = "set_array_from_python(T *) -> "; + PyObject *numpy_array_obj = PyDict_GetItemString (python_data, python_key); + if (numpy_array_obj) { + Python3_Numpy np; + np.set(numpy_array_obj); + pointdata_from_np_array(np, out); + mlog << Debug(7) << method_name + << "get the point data for " << python_key << " from python object\n"; + } + else { + if (required) { + mlog << Error << "\n" << method_name + << "error getting the point data by the key (" << python_key << ") from python object\n\n"; + exit (1); + } + else mlog << Debug(3) << method_name + << "not exists the point data (" << python_key << ") from python object\n"; + } +} + +/* +static void set_array_from_python(PyObject *python_data, const char *python_key, int *out, bool required=true) { + const char *method_name = "set_array_from_python(T *) -> "; + PyObject *numpy_array_obj = PyDict_GetItemString (python_data, python_key); + if (numpy_array_obj) { + Python3_Numpy np; +mlog << Debug(7) << method_name + << "get the point data for " << python_key << " from python object - before\n"; + np.set(numpy_array_obj); +mlog << Debug(7) << method_name + << "get the point data for " << python_key << " from python object - after np.set\n"; + pointdata_from_np_array(np, out); + mlog << Debug(7) << method_name + << "get the point data for " << python_key << " from python object\n"; + } + else { + if (required) { + mlog << Error << "\n" << method_name + << "error getting the point data by the key (" << python_key << ") from python object\n\n"; + exit (1); + } + else mlog << Debug(3) << method_name + << "not exists the point data (" << python_key << ") from python object\n"; + } +} +*/ + +/* +static void set_array_from_python(PyObject *python_data, const char *python_key, int *out, bool required=true) { + const char *method_name = "set_array_from_python(int *) -> "; + PyObject *numpy_array_obj = PyDict_GetItemString (python_data, python_key); + if (numpy_array_obj) { + Python3_Numpy np; + np.set(numpy_array_obj); + pointdata_from_np_array(np, out); + } + else if (required) { + mlog << Error << "\nset_array_from_python(int *) -> " + << "error getting member (" << python_key << ") from python object\"\n\n"; + exit (1); + } + else mlog << Debug(3) << method_name + << "not exists the point data (" << python_key << ") from python object\n"; +} + +//////////////////////////////////////////////////////////////////////// + +static void set_array_from_python(PyObject *python_data, const char *python_key, float *out, bool required=true) { + const char *method_name = "set_array_from_python(float *) -> "; + PyObject *numpy_array_obj = PyDict_GetItemString (python_data, python_key); + if (numpy_array_obj) { + Python3_Numpy np; + np.set(numpy_array_obj); + pointdata_from_np_array(np, out); + mlog << Debug(7) << method_name + << "get the point data for " << python_key << " from python object\n"; + } + else if (required) { + mlog << Error << "\n" << method_name + << "error getting member (" << python_key << ") from python object\"\n\n"; + exit (1); + } + else mlog << Debug(3) << method_name + << "not exists the point data (" << python_key << ") from python object\n"; +} +*/ +//////////////////////////////////////////////////////////////////////// + +/* +static void set_met_array_from_python(PyObject *python_data, const char *python_key, IntArray *out, bool required=true) { + const char *method_name = "set_met_array_from_python(IntArray *) -> "; + PyObject *numpy_array_obj = PyDict_GetItemString (python_data, python_key); + if (numpy_array_obj) { + Python3_Numpy np; + np.set(numpy_array_obj); + pointdata_from_np_array(np, out); + mlog << Debug(7) << method_name + << "get the point data for " << python_key << " from python object\n"; + } + else if (required) { + mlog << Error << "\n" << method_name + << "error getting member (" << python_key << ") from python object\"\n\n"; + exit (1); + } + else mlog << Debug(3) << method_name + << "not exists the point data (" << python_key << ") from python object\n"; +} +*/ +/* +//////////////////////////////////////////////////////////////////////// + +static void set_met_array_from_python(PyObject *python_data, const char *python_key, NumArray &out) { + PyObject *numpy_array_obj = PyDict_GetItemString (python_data, python_key); + if (numpy_array_obj) { + Python3_Numpy np; + np.set(numpy_array_obj); + pointdata_from_np_array(np, out); + } + else { + mlog << Error << "\nset_array_from_python(NumArray) -> " + << "error getting member (" << python_key << ") from python object\"\n\n"; + exit (1); + } +} +*/ +//////////////////////////////////////////////////////////////////////// + +static void set_str_array_from_python(PyObject *python_data, const char *python_key, StringArray *out) { + const char *method_name = "set_met_array_from_python(StringArray *) -> "; + PyObject *str_array_obj = PyDict_GetItemString (python_data, python_key); + if (str_array_obj) { + pointdata_from_str_array(str_array_obj, out); + mlog << Debug(7) << method_name + << "get the point data for " << python_key << " from python object\n"; + } + else { + mlog << Error << "\n" << method_name + << "error getting member (" << python_key << ") from python object\"\n\n"; + exit (1); + } +} + + +//////////////////////////////////////////////////////////////////////// + + +bool python_point_data(const char * script_name, int script_argc, char ** script_argv, + const bool use_xarray, MetPointDataPython &met_pd_out) + +{ + +bool status = straight_python_point_data(script_name, script_argc, script_argv, + use_xarray, met_pd_out); + +return ( status ); + +} + +//////////////////////////////////////////////////////////////////////// + + +bool straight_python_point_data(const char * script_name, int script_argc, char ** script_argv, + const bool use_xarray, MetPointDataPython &met_pd_out) +{ + +int int_value; +PyObject * module_obj = 0; +PyObject * module_dict_obj = 0; +PyObject * python_key = 0; +PyObject * python_value = 0; +PyObject * numpy_array_obj = 0; +PyObject * python_met_point_data = 0; +ConcatString cs, user_dir, user_base; +const char *method_name = "straight_python_point_data -> "; +const char *method_name_s = "straight_python_point_data()"; + +mlog << Debug(3) << "Running user's python script (" + << script_name << ").\n"; + +cs = script_name; +user_dir = cs.dirname(); +user_base = cs.basename(); + +Wchar_Argv wa; + +wa.set(script_argc, script_argv); + + // + // if the global python object has already been initialized, + // we need to reload the module + // + +bool do_reload = GP.is_initialized; + +GP.initialize(); + + // + // start up the python interpreter + // + +if ( PyErr_Occurred() ) { + + PyErr_Print(); + + mlog << Warning << "\n" << method_name + << "an error occurred initializing python\n\n"; + + return ( false ); + +} + + // + // set the arguments + // + +run_python_string("import os"); +run_python_string("import sys"); + +ConcatString command; + +command << cs_erase + << "sys.path.append(\"" + << user_dir + << "\")"; + +run_python_string(command.text()); + +if ( script_argc > 0 ) { + + PySys_SetArgv (wa.wargc(), wa.wargv()); + +} + + // + // import the python script as a module + // + +module_obj = PyImport_ImportModule (user_base.c_str()); + + // + // if needed, reload the module + // + +if ( do_reload ) { + + module_obj = PyImport_ReloadModule (module_obj); + +} + +if ( PyErr_Occurred() ) { + + PyErr_Print(); + + mlog << Warning << "\n" << method_name + << "an error occurred importing module \"" + << script_name << "\"\n\n"; + + return ( false ); + +} + +if ( ! module_obj ) { + + mlog << Warning << "\n" << method_name + << "error running python script \"" + << script_name << "\"\n\n"; + + return ( false ); + +} + + // + // get the namespace for the module (as a dictionary) + // + +module_dict_obj = PyModule_GetDict (module_obj); + + // + // get handles to the objects of interest from the module_dict + // + +python_met_point_data = PyDict_GetItemString (module_dict_obj, python_key_point_data); + +python_value = PyDict_GetItemString (python_met_point_data, python_use_var_id); + +bool use_var_id = pyobject_as_bool(python_value); +met_pd_out.set_use_var_id(use_var_id); +mlog << Debug(9) << method_name << "use_var_id: \"" << use_var_id << "\" from python. is_using_var_id(): " << met_pd_out.is_using_var_id() << "\n"; + +python_value = PyDict_GetItemString (python_met_point_data, python_key_nhdr); + +int_value = pyobject_as_int(python_value); +met_pd_out.set_hdr_cnt(int_value); + +python_value = PyDict_GetItemString (python_met_point_data, python_key_nobs); +int_value = pyobject_as_int(python_value); + +met_pd_out.allocate(int_value); + +MetPointObsData *obs_data = met_pd_out.get_point_obs_data(); +MetPointHeader *header_data = met_pd_out.get_header_data(); + +if ( use_xarray ) { + + PyObject * data_array_obj = 0; + + // look up the data array variable name from the dictionary + + data_array_obj = PyDict_GetItemString (python_met_point_data, numpy_array_hdr_typ); + + if ( ! data_array_obj ) { + + mlog << Warning << "\n" << method_name + << "trouble reading data from \"" + << script_name << "\"\n\n"; + + return ( false ); + } + + //pointdata_from_xarray(data_array_obj, &header_data.hdr_typ); + +} else { // numpy array & dict + + // look up the data array variable name from the dictionary + +/* + set_met_array_from_python(python_met_point_data, numpy_array_hdr_typ, header_data->typ_idx_array); + set_met_array_from_python(python_met_point_data, numpy_array_hdr_sid, header_data->sid_idx_array); + set_met_array_from_python(python_met_point_data, numpy_array_hdr_vld, header_data->vld_idx_array); + //IntArray vld_num_array; // number array for valid time + + set_met_array_from_python(python_met_point_data, numpy_array_hdr_lat, header_data->lat_array); + set_met_array_from_python(python_met_point_data, numpy_array_hdr_lon, header_data->lon_array); + set_met_array_from_python(python_met_point_data, numpy_array_hdr_elv, header_data->elv_array); +*/ + set_array_from_python(python_met_point_data, numpy_array_hdr_typ, &header_data->typ_idx_array); + set_array_from_python(python_met_point_data, numpy_array_hdr_sid, &header_data->sid_idx_array); + set_array_from_python(python_met_point_data, numpy_array_hdr_vld, &header_data->vld_idx_array); + set_array_from_python(python_met_point_data, numpy_array_hdr_lat, &header_data->lat_array); + set_array_from_python(python_met_point_data, numpy_array_hdr_lon, &header_data->lon_array); + set_array_from_python(python_met_point_data, numpy_array_hdr_elv, &header_data->elv_array); + + set_str_array_from_python(python_met_point_data, numpy_array_hdr_typ_table, &header_data->typ_array); + set_str_array_from_python(python_met_point_data, numpy_array_hdr_sid_table, &header_data->sid_array); + set_str_array_from_python(python_met_point_data, numpy_array_hdr_vld_table, &header_data->vld_array); + set_array_from_python(python_met_point_data, numpy_array_prpt_typ_table, &header_data->prpt_typ_array, false); + set_array_from_python(python_met_point_data, numpy_array_irpt_typ_table, &header_data->irpt_typ_array, false); + set_array_from_python(python_met_point_data, numpy_array_inst_typ_table, &header_data->inst_typ_array, false); + + +// if ( !numpy_array_obj ) { +// +// mlog << Warning << "\n" << method_name +// << "trouble reading data from \"" +// << script_name << "\"\n\n"; +// +//// return ( false ); +// } +/* +*/ + set_array_from_python(python_met_point_data, numpy_array_obs_qty, obs_data->obs_qids); + set_array_from_python(python_met_point_data, numpy_array_obs_hid, obs_data->obs_hids); + set_array_from_python(python_met_point_data, numpy_array_obs_vid, obs_data->obs_ids); + set_array_from_python(python_met_point_data, numpy_array_obs_lvl, obs_data->obs_lvls); + set_array_from_python(python_met_point_data, numpy_array_obs_hgt, obs_data->obs_hgts); + set_array_from_python(python_met_point_data, numpy_array_obs_val, obs_data->obs_vals); + + set_str_array_from_python(python_met_point_data, numpy_array_obs_qty_table, &obs_data->qty_names); + set_str_array_from_python(python_met_point_data, numpy_array_obs_var_table, &obs_data->var_names); + + if(mlog.verbosity_level()>=point_data_debug_level) print_met_data(obs_data, header_data, method_name_s); + +} + + // + // done + // + +return ( true ); + +} + + +//////////////////////////////////////////////////////////////////////// + +void print_met_data(MetPointObsData *obs_data, MetPointHeader *header_data, + const char *caller, int debug_level) { + int log_count, count; + const int min_count = 20; + const char *method_name = "print_met_data() "; + + mlog << Debug(debug_level) << "\n" << method_name << "by " << caller << "\n" + << " obs_data.obs_cnt = " << obs_data->obs_cnt << "\n" + << " header_data.hdr_count = " << header_data->hdr_count << " type=" + << header_data->typ_idx_array.n() << ", sid=" + << header_data->sid_idx_array.n() << ", valid=" + << header_data->vld_idx_array.n() << ", lat=" + << header_data->lat_array.n() << ", lon=" + << header_data->lon_array.n() << ", elv=" + << header_data->elv_array.n() << ", message_type=" + << header_data->typ_array.n() << ", station_id=" + << header_data->sid_array.n() << ", valid_time=" + << header_data->vld_array.n() << ", prpt=" + << header_data->prpt_typ_array.n() << ", irpt=" + << header_data->irpt_typ_array.n() << ", inst=" + << header_data->inst_typ_array.n() << " &header_data=" << header_data + << "\n"; + + log_count = (header_data->hdr_count > min_count) ? min_count : header_data->hdr_count; + mlog << Debug(debug_level) << method_name + << "header_data: message_type,station_id,time_time,lat,lon.elv\n"; + for (int idx=0; idxsid_idx_array[idx] << ", " + << header_data->vld_idx_array[idx] << ", " + << header_data->lat_array[idx] << ", " + << header_data->lon_array[idx] << ", " + << header_data->elv_array[idx] << "\n"; + } + if (header_data->hdr_count > log_count) { + log_count = header_data->hdr_count - min_count; + if (log_count < min_count) log_count = min_count; + else mlog << Debug(debug_level) + << " header_data[...] = ...\n"; + for (int idx=log_count; idxhdr_count; idx++) { + mlog << Debug(debug_level) + << " header_data[" << idx << "] = " + << header_data->typ_idx_array[idx] << ", " + << header_data->sid_idx_array[idx] << ", " + << header_data->vld_idx_array[idx] << ", " + << header_data->lat_array[idx] << ", " + << header_data->lon_array[idx] << ", " + << header_data->elv_array[idx] << "\n"; + } + } + if (header_data->typ_array.n() > 0) { + mlog << Debug(debug_level) << "\n"; + count = header_data->typ_array.n(); + for (int idx=0; idxsid_array.n() > 0) { + mlog << Debug(debug_level) << "\n"; + count = header_data->sid_array.n(); + log_count = (count > min_count) ? min_count : count; + for (int idx=0; idx log_count) { + log_count = count - min_count; + if (log_count < min_count) log_count = min_count; + else mlog << Debug(debug_level) + << " station_id[...] = ...\n"; + for (int idx=log_count; idxvld_array.n() > 0) { + mlog << Debug(debug_level) << "\n"; + count = header_data->vld_array.n(); + log_count = (count > min_count) ? min_count : count; + for (int idx=0; idx log_count) { + log_count = count - min_count; + if (log_count < min_count) log_count = min_count; + else mlog << Debug(debug_level) + << " valid time[...] = ...\n"; + for (int idx=log_count; idxprpt_typ_array.n() > 0) { + mlog << Debug(debug_level) << "\n"; + for (int idx=0; idxprpt_typ_array.n(); idx++) + mlog << Debug(debug_level) + << " prpt_type[" << idx << "] = " << header_data->prpt_typ_array[idx] << "\n"; + } + if (header_data->irpt_typ_array.n() > 0) { + mlog << Debug(debug_level) << "\n"; + for (int idx=0; idxirpt_typ_array.n(); idx++) + mlog << Debug(debug_level) + << " irpt_type[" << idx << "] = " << header_data->irpt_typ_array[idx] << "\n"; + } + if (header_data->inst_typ_array.n() > 0) { + mlog << Debug(debug_level) << "\n"; + for (int idx=0; idxinst_typ_array.n(); idx++) + mlog << Debug(debug_level) + << " inst_type[" << idx << "] = " << header_data->inst_typ_array[idx] << "\n"; + } + + log_count = (obs_data->obs_cnt > min_count) ? min_count : obs_data->obs_cnt; + mlog << Debug(debug_level) << "\n" << method_name + << "obs_data: hid,vid.level,height,value,qty\n"; + for (int idx=0; idxobs_ids[idx] << ", " + << obs_data->obs_lvls[idx] << ", " + << obs_data->obs_hgts[idx] << ", " + << obs_data->obs_vals[idx] << ", " + << obs_data->obs_qids[idx] << "\n"; + } + if (obs_data->obs_cnt > log_count) { + log_count = obs_data->obs_cnt - min_count; + if (log_count < min_count) log_count = min_count; + else mlog << Debug(debug_level) + << " obs_data[...] = ...\n"; + for (int idx=log_count; idxobs_cnt; idx++) { + mlog << Debug(debug_level) + << " obs_data[" << idx << "] = " + << obs_data->obs_hids[idx] << ", " + << obs_data->obs_ids[idx] << ", " + << obs_data->obs_lvls[idx] << ", " + << obs_data->obs_hgts[idx] << ", " + << obs_data->obs_vals[idx] << ", " + << obs_data->obs_qids[idx] << "\n"; + } + } + + if (obs_data->var_names.n() > 0) { + mlog << Debug(debug_level) << "\n"; + for (int idx=0; idxvar_names.n(); idx++) + mlog << Debug(debug_level) + << " var_names[" << idx << "] = " << obs_data->var_names[idx] << "\n"; + } + else mlog << Debug(debug_level) + << " var_names is empty!!!\n"; + if (obs_data->qty_names.n() > 0) { + mlog << Debug(debug_level) << "\n"; + for (int idx=0; idxqty_names.n(); idx++) + mlog << Debug(debug_level) + << " qty_names[" << idx << "] = " << obs_data->qty_names[idx] << "\n"; + } + else mlog << Debug(debug_level) + << " qty_names is empty!!!\n"; + + mlog << Debug(debug_level) << "Done " << method_name << "by " << caller << "\n\n"; + +} + +//////////////////////////////////////////////////////////////////////// diff --git a/met/src/libcode/vx_pointdata_python/python_pointdata.h b/met/src/libcode/vx_pointdata_python/python_pointdata.h new file mode 100644 index 0000000000..5be19111f0 --- /dev/null +++ b/met/src/libcode/vx_pointdata_python/python_pointdata.h @@ -0,0 +1,80 @@ +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* +// ** Copyright UCAR (c) 1992 - 2021 +// ** University Corporation for Atmospheric Research (UCAR) +// ** National Center for Atmospheric Research (NCAR) +// ** Research Applications Lab (RAL) +// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA +// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* + +//////////////////////////////////////////////////////////////////////// + + +#ifndef __PYTHON_POINTDATA__ +#define __PYTHON_POINTDATA__ + + +//////////////////////////////////////////////////////////////////////// + + +#include "met_point_data.h" + + +extern "C" { + +#include "Python.h" + +} + + +//////////////////////////////////////////////////////////////////////// + +static const char python_key_point_data [] = "met_point_data"; + +static const char python_key_nhdr [] = "nhdr"; +//static const char python_key_npbhdr [] = "npbhdr"; +static const char python_use_var_id [] = "use_var_id"; +static const char numpy_array_hdr_typ [] = "hdr_typ"; // message type IDs +static const char numpy_array_hdr_sid [] = "hdr_sid"; // station IDs +static const char numpy_array_hdr_vld [] = "hdr_vld"; // valid time IDs +static const char numpy_array_hdr_lat [] = "hdr_lat"; +static const char numpy_array_hdr_lon [] = "hdr_lon"; +static const char numpy_array_hdr_elv [] = "hdr_elv"; +static const char numpy_array_hdr_typ_table [] = "hdr_typ_table"; // message type list +static const char numpy_array_hdr_sid_table [] = "hdr_sid_table"; // station ID list +static const char numpy_array_hdr_vld_table [] = "hdr_vld_table"; // valid time list +static const char numpy_array_prpt_typ_table[] = "prpt_typ_table"; +static const char numpy_array_irpt_typ_table[] = "irpt_typ_table"; +static const char numpy_array_inst_typ_table[] = "inst_typ_table"; + +static const char python_key_nobs [] = "nobs"; +static const char numpy_array_obs_qty [] = "obs_qty"; // quality_id +static const char numpy_array_obs_hid [] = "obs_hid"; // header id +static const char numpy_array_obs_vid [] = "obs_vid"; // variable id or grib code +static const char numpy_array_obs_lvl [] = "obs_lvl"; +static const char numpy_array_obs_hgt [] = "obs_hgt"; +static const char numpy_array_obs_val [] = "obs_val"; +static const char numpy_array_obs_qty_table [] = "obs_qty_table"; +static const char numpy_array_obs_var_table [] = "obs_var_table"; // variable names or grib codes as string + +static const int point_data_debug_level = 10; + + +//////////////////////////////////////////////////////////////////////// + + +extern bool python_point_data(const char * script_name, int script_argc, char ** script_argv, + const bool use_xarray, MetPointDataPython &met_pd_out); +//extern bool python_point_data(const char *python_command, const bool use_xarray, +// MetPointData & po_out); +extern void print_met_data(MetPointObsData *obs_data, MetPointHeader *header_data, + const char *caller, int debug_level=point_data_debug_level); + + +//////////////////////////////////////////////////////////////////////// + + +#endif /* __PYTHON_POINTDATA__ */ + + +//////////////////////////////////////////////////////////////////////// + diff --git a/met/src/libcode/vx_python3_utils/python3_util.cc b/met/src/libcode/vx_python3_utils/python3_util.cc index fe56b6f6b2..ab17acd7b8 100644 --- a/met/src/libcode/vx_python3_utils/python3_util.cc +++ b/met/src/libcode/vx_python3_utils/python3_util.cc @@ -12,7 +12,6 @@ using namespace std; #include #include "vx_log.h" -#include "concat_string.h" #include "vx_math.h" #include "python3_util.h" @@ -125,6 +124,18 @@ return ( k ); //////////////////////////////////////////////////////////////////////// +bool pyobject_as_bool (PyObject * obj) + +{ + +return ( 1 == PyObject_IsTrue(obj) ); + +} + + +//////////////////////////////////////////////////////////////////////// + + double pyobject_as_double (PyObject * obj) { @@ -215,6 +226,29 @@ return ( s ); //////////////////////////////////////////////////////////////////////// +StringArray pyobject_as_string_array (PyObject * obj) + +{ + +StringArray a; +ConcatString s; +PyObject *item = 0; + +int size = PyList_Size (obj); +for (int idx=0; idxget_hdr_cnt(); + int obs_count = met_point_obs->get_obs_cnt(); mlog << Debug(2) << "Searching " << (obs_count) << " observations from " << (hdr_count) << " header messages.\n"; - int buf_size = ((obs_count > DEF_NC_BUFFER_SIZE) ? DEF_NC_BUFFER_SIZE : (obs_count)); + const int buf_size = ((obs_count > DEF_NC_BUFFER_SIZE) ? DEF_NC_BUFFER_SIZE : (obs_count)); int obs_qty_idx_block[buf_size]; float obs_arr_block[buf_size][OBS_ARRAY_LEN]; @@ -1107,19 +1146,25 @@ void process_point_obs(int i_nc) { ConcatString obs_qty_str; ConcatString var_name; StringArray var_names; - StringArray obs_qty_array = nc_point_obs.get_qty_data(); - if(use_var_id) var_names = nc_point_obs.get_var_names(); + StringArray obs_qty_array = met_point_obs->get_qty_data(); + if(use_var_id) var_names = met_point_obs->get_var_names(); for(int i_start=0; i_start DEF_NC_BUFFER_SIZE) - ? DEF_NC_BUFFER_SIZE : (obs_count-i_start); + int buf_size2 = (obs_count-i_start); + if (buf_size2 > buf_size) buf_size2 = buf_size; + +#ifdef WITH_PYTHON + if (use_python) + status = met_point_obs->get_point_obs_data()->fill_obs_buf( + buf_size2, i_start, (float *)obs_arr_block, obs_qty_idx_block); + else +#endif + status = nc_point_obs.read_obs_data(buf_size2, i_start, (float *)obs_arr_block, + obs_qty_idx_block, (char *)0); + if (!status) exit(1); - if (!nc_point_obs.read_obs_data(buf_size, i_start, (float *)obs_arr_block, - obs_qty_idx_block, (char *)0)) { - exit(1); - } // Process each observation in the file - for(int i_offset=0; i_offsetget_header_offset(obs_arr); // Range check the header offset if(headerOffset < 0 || headerOffset >= hdr_count) { @@ -1144,16 +1189,21 @@ void process_point_obs(int i_nc) { // Read the corresponding header array for this observation // - the corresponding header type, header Station ID, and valid time - nc_point_obs.get_header(headerOffset, hdr_arr, hdr_typ_str, - hdr_sid_str, hdr_vld_str); +#ifdef WITH_PYTHON + if (use_python) + met_point_obs->get_header(headerOffset, hdr_arr, hdr_typ_str, hdr_sid_str, hdr_vld_str); + else +#endif + nc_point_obs.get_header(headerOffset, hdr_arr, hdr_typ_str, + hdr_sid_str, hdr_vld_str); // Read the header integer types - nc_point_obs.get_header_type(headerOffset, hdr_typ_arr); + met_point_obs->get_header_type(headerOffset, hdr_typ_arr); // Convert string to a unixtime hdr_ut = timestring_to_unix(hdr_vld_str.c_str()); - int grib_code = nc_point_obs.get_grib_code_or_var_index(obs_arr); + int grib_code = met_point_obs->get_grib_code_or_var_index(obs_arr); if (use_var_id && grib_code < var_names.n()) { var_name = var_names[grib_code]; } @@ -1176,6 +1226,10 @@ void process_point_obs(int i_nc) { } // end for i_start // Deallocate and clean up +#ifdef WITH_PYTHON + if (use_python) met_point_file.close(); + else +#endif nc_point_obs.close(); return; diff --git a/met/src/tools/core/point_stat/point_stat.cc b/met/src/tools/core/point_stat/point_stat.cc index f62f3ae511..48ffd21c6e 100644 --- a/met/src/tools/core/point_stat/point_stat.cc +++ b/met/src/tools/core/point_stat/point_stat.cc @@ -127,6 +127,11 @@ using namespace std; #include "nc_obs_util.h" #include "nc_point_obs_in.h" +#ifdef WITH_PYTHON +#include "data2d_nc_met.h" +#include "pointdata_python.h" +#endif + //////////////////////////////////////////////////////////////////////// #define BUFFER_SIZE (DEF_NC_BUFFER_SIZE/2) @@ -672,49 +677,88 @@ void process_obs_file(int i_nc) { // Open the observation file as a NetCDF file. // The observation file must be in NetCDF format as the // output of the PB2NC or ASCII2NC tool. + bool status; + bool use_var_id = true; + bool use_arr_vars = false; + bool use_python = false; MetNcPointObsIn nc_point_obs; - if( !nc_point_obs.open(obs_file[i_nc].c_str()) ) { - nc_point_obs.close(); + MetPointData *met_point_obs = 0; +#ifdef WITH_PYTHON + MetPythonPointDataFile met_point_file; + string python_command = obs_file[i_nc]; + bool use_xarray = (0 == python_command.find(conf_val_python_xarray)); + use_python = use_xarray || (0 == python_command.find(conf_val_python_numpy)); + if (use_python) { + int offset = python_command.find("="); + if (offset == std::string::npos) { + mlog << Error << "\n" << method_name + << "trouble parsing the python command " << python_command << ".\n\n"; + exit(1); + } - mlog << Warning << "\n" << method_name - << "can't open observation netCDF file: " - << obs_file[i_nc] << "\n\n"; - return; - } + if(!met_point_file.open(python_command.substr(offset+1).c_str(), use_xarray)) { + met_point_file.close(); + mlog << Error << "\n" << method_name + << "trouble getting point observation file from python command " + << python_command << ".\n\n"; + exit(1); + } - nc_point_obs.read_dim_headers(); - nc_point_obs.check_nc(obs_file[i_nc].c_str(), method_name); - nc_point_obs.read_obs_data_table_lookups(); + met_point_obs = met_point_file.get_met_point_data(); + } + else { +#endif + if( !nc_point_obs.open(obs_file[i_nc].c_str()) ) { + nc_point_obs.close(); + + mlog << Warning << "\n" << method_name + << "can't open observation netCDF file: " + << obs_file[i_nc] << "\n\n"; + return; + } + + nc_point_obs.read_dim_headers(); + nc_point_obs.check_nc(obs_file[i_nc].c_str(), method_name); + nc_point_obs.read_obs_data_table_lookups(); + met_point_obs = (MetPointData *)&nc_point_obs; + use_var_id = nc_point_obs.is_using_var_id(); + use_arr_vars = nc_point_obs.is_using_obs_arr(); +#ifdef WITH_PYTHON + } +#endif - bool use_var_id = nc_point_obs.is_using_var_id(); - int hdr_count = nc_point_obs.get_hdr_cnt(); - int obs_count = nc_point_obs.get_obs_cnt(); + int hdr_count = met_point_obs->get_hdr_cnt(); + int obs_count = met_point_obs->get_obs_cnt(); mlog << Debug(2) << "Searching " << obs_count << " observations from " << hdr_count << " messages.\n"; ConcatString var_name(""); - bool use_arr_vars = nc_point_obs.is_using_obs_arr(); StringArray var_names; - StringArray obs_qty_array = nc_point_obs.get_qty_data(); - if( use_var_id ) var_names = nc_point_obs.get_var_names(); + StringArray obs_qty_array = met_point_obs->get_qty_data(); + if( use_var_id ) var_names = met_point_obs->get_var_names(); - int buf_size = ((obs_count > BUFFER_SIZE) ? BUFFER_SIZE : (obs_count)); + const int buf_size = ((obs_count > BUFFER_SIZE) ? BUFFER_SIZE : (obs_count)); int obs_qty_idx_block[buf_size]; float obs_arr_block[buf_size][OBS_ARRAY_LEN]; // Process each observation in the file int str_length, block_size; - for(int i_block_start_idx=0; i_block_start_idx BUFFER_SIZE) block_size = BUFFER_SIZE; - - if (!nc_point_obs.read_obs_data(block_size, i_block_start_idx, - (float *)obs_arr_block, - obs_qty_idx_block, (char *)0)) { - exit(1); - } + if (block_size > buf_size) block_size = buf_size; + +#ifdef WITH_PYTHON + if (use_python) + status = met_point_obs->get_point_obs_data()->fill_obs_buf( + block_size, i_block_start_idx, (float *)obs_arr_block, obs_qty_idx_block); + else +#endif + status = nc_point_obs.read_obs_data(block_size, i_block_start_idx, + (float *)obs_arr_block, + obs_qty_idx_block, (char *)0); + if (!status) exit(1); int hdr_idx; for(int i_block_idx=0; i_block_idxget_header_offset(obs_arr); // Range check the header offset if(headerOffset < 0 || headerOffset >= hdr_count) { @@ -741,11 +785,16 @@ void process_obs_file(int i_nc) { // Read the corresponding header array for this observation // - the corresponding header type, header Station ID, and valid time - nc_point_obs.get_header(headerOffset, hdr_arr, hdr_typ_str, - hdr_sid_str, hdr_vld_str); +#ifdef WITH_PYTHON + if (use_python) + met_point_obs->get_header(headerOffset, hdr_arr, hdr_typ_str, hdr_sid_str, hdr_vld_str); + else +#endif + nc_point_obs.get_header(headerOffset, hdr_arr, hdr_typ_str, + hdr_sid_str, hdr_vld_str); // Store the variable name - int org_grib_code = nc_point_obs.get_grib_code_or_var_index(obs_arr); + int org_grib_code = met_point_obs->get_grib_code_or_var_index(obs_arr); int grib_code = org_grib_code; if (use_var_id && grib_code < var_names.n()) { var_name = var_names[grib_code]; @@ -774,7 +823,7 @@ void process_obs_file(int i_nc) { // and at the same vertical level. if(vflag && is_vgrd) { - if(!nc_point_obs.is_same_obs_values(obs_arr, prev_obs_arr)) { + if(!met_point_obs->is_same_obs_values(obs_arr, prev_obs_arr)) { mlog << Error << "\n" << method_name << "for observation index " << i_obs << ", when computing VL1L2 and/or VAL1L2 vector winds " @@ -802,12 +851,16 @@ void process_obs_file(int i_nc) { grid, var_name.c_str()); } - nc_point_obs.set_grib_code_or_var_index(obs_arr, org_grib_code); + met_point_obs->set_grib_code_or_var_index(obs_arr, org_grib_code); } } // end for i_block_start_idx // Deallocate and clean up +#ifdef WITH_PYTHON + if (use_python) met_point_file.close(); + else +#endif nc_point_obs.close(); return; diff --git a/met/src/tools/other/plot_point_obs/plot_point_obs.cc b/met/src/tools/other/plot_point_obs/plot_point_obs.cc index 1ba1b9938b..243f4d0c29 100644 --- a/met/src/tools/other/plot_point_obs/plot_point_obs.cc +++ b/met/src/tools/other/plot_point_obs/plot_point_obs.cc @@ -45,26 +45,11 @@ using namespace std; #include "plot_point_obs.h" -#include "nc_utils.h" -#include "vx_log.h" -#include "data_plane.h" #include "data_plane_plot.h" -#include "write_netcdf.h" -#include "vx_data2d.h" -#include "vx_data2d_factory.h" -#include "vx_data2d_grib.h" -#include "vx_data2d_nc_met.h" -#include "vx_data2d_nc_pinterp.h" -#include "vx_util.h" -#include "vx_cal.h" -#include "vx_grid.h" -#include "vx_math.h" -#include "vx_color.h" #include "vx_ps.h" #include "vx_pxm.h" #include "vx_render.h" #include "vx_plot_util.h" -#include "nc_obs_util.h" #include "nc_point_obs_in.h" //////////////////////////////////////////////////////////////////////// @@ -149,42 +134,75 @@ int main(int argc, char *argv[]) { void process_point_obs(const char *point_obs_filename) { int h, v, i_obs; - const char *method_name = "\nprocess_point_obs() -> "; - const char *method_name_s = "\nprocess_point_obs() "; + const char *method_name = "process_point_obs() -> "; + const char *method_name_s = "process_point_obs() "; // Open the netCDF point observation file mlog << Debug(1) << "Reading point observation file: " << point_obs_filename << "\n"; + bool status; + bool use_var_id = true; + bool use_obs_arr = false; + + bool use_python = false; MetNcPointObsIn nc_point_obs; - if(!nc_point_obs.open(point_obs_filename)) { - nc_point_obs.close(); + MetPointData *met_point_obs = 0; +#ifdef WITH_PYTHON + MetPythonPointDataFile met_point_file; + string python_command = point_obs_filename; + bool use_xarray = (0 == python_command.find(conf_val_python_xarray)); + use_python = use_xarray || (0 == python_command.find(conf_val_python_numpy)); + if (use_python) { + int offset = python_command.find("="); + if (offset == std::string::npos) { + mlog << Error << "\n" << method_name + << "trouble parsing the python command " << python_command << ".\n\n"; + exit(1); + } - mlog << Error << "\n" << method_name - << "trouble opening point observation file " - << point_obs_filename << ".\n\n"; + if(!met_point_file.open(python_command.substr(offset+1).c_str(), use_xarray)) { + met_point_file.close(); + mlog << Error << "\n" << method_name + << "trouble getting point observation file from python command " + << python_command << ".\n\n"; + exit(1); + } - exit(1); + met_point_obs = met_point_file.get_met_point_data(); } + else +#endif + { + if(!nc_point_obs.open(point_obs_filename)) { + nc_point_obs.close(); - // Read the dimensions and variables - nc_point_obs.read_dim_headers(); - nc_point_obs.check_nc(point_obs_filename, method_name_s); // exit if missing dims/vars - nc_point_obs.read_obs_data(); + mlog << Error << "\n" << method_name + << "trouble opening point observation file " + << point_obs_filename << ".\n\n"; - bool use_var_id = nc_point_obs.is_using_var_id(); + exit(1); + } - // Print warning about ineffective command line arguments - if(use_var_id && ivar.n() != 0) { - mlog << Warning << "\n" << method_name << "-gc option is ignored!\n\n"; - } - if(!use_var_id && svar.n() != 0) { - mlog << Warning << "\n" << method_name << "-obs_var option is ignored!\n\n"; - } + // Read the dimensions and variables + nc_point_obs.read_dim_headers(); + nc_point_obs.check_nc(point_obs_filename, method_name_s); // exit if missing dims/vars + nc_point_obs.read_obs_data(); + use_obs_arr = nc_point_obs.is_using_obs_arr(); - long nhdr_count = nc_point_obs.get_hdr_cnt(); - long nobs_count = nc_point_obs.get_obs_cnt(); + // Print warning about ineffective command line arguments + if(use_var_id && ivar.n() != 0) { + mlog << Warning << "\n" << method_name << "-gc option is ignored!\n\n"; + } + if(!use_var_id && svar.n() != 0) { + mlog << Warning << "\n" << method_name << "-obs_var option is ignored!\n\n"; + } + met_point_obs = (MetPointData *)&nc_point_obs; + } + + long nhdr_count = met_point_obs->get_hdr_cnt(); + long nobs_count = met_point_obs->get_obs_cnt(); mlog << Debug(2) << "Processing " << nobs_count << " observations at " << nhdr_count << " locations.\n"; @@ -192,8 +210,6 @@ void process_point_obs(const char *point_obs_filename) { // Get the corresponding header: // message type, staton_id, valid_time, and lat/lon/elv - bool use_obs_arr = nc_point_obs.is_using_obs_arr(); - int buf_size = (nobs_count > DEF_NC_BUFFER_SIZE) ? DEF_NC_BUFFER_SIZE : nobs_count; // Allocate space to store the data @@ -202,8 +218,18 @@ void process_point_obs(const char *point_obs_filename) { int obs_qty_block[buf_size]; float obs_arr_block[buf_size][OBS_ARRAY_LEN]; - if(use_var_id) var_list = nc_point_obs.get_var_names(); - qty_list = nc_point_obs.get_qty_data(); + use_var_id = met_point_obs->is_using_var_id(); + if(use_var_id) var_list = met_point_obs->get_var_names(); + qty_list = met_point_obs->get_qty_data(); + + if (0 == qty_list.n()) { + mlog << Error << "\n" << method_name << "Missing quality marks\n\n"; + exit (1); + } + if (use_var_id && (0 == var_list.n())) { + mlog << Debug(1) << method_name << "Missing variable names\ns\n"; + exit (1); + } ConcatString hdr_typ_str; ConcatString hdr_sid_str; @@ -211,16 +237,21 @@ void process_point_obs(const char *point_obs_filename) { ConcatString obs_qty_str; for(int i_start=0; i_start DEF_NC_BUFFER_SIZE) ? - DEF_NC_BUFFER_SIZE : (nobs_count-i_start); - - if (!nc_point_obs.read_obs_data(buf_size, i_start, (float *)obs_arr_block, - obs_qty_block, (char *)0)) { - exit(1); - } + int buf_size2 = ((nobs_count-i_start) > DEF_NC_BUFFER_SIZE) ? + DEF_NC_BUFFER_SIZE : (nobs_count-i_start); + +#ifdef WITH_PYTHON + if (use_python) + status = met_point_obs->get_point_obs_data()->fill_obs_buf( + buf_size2, i_start, (float *)obs_arr_block, obs_qty_block); + else +#endif + status = nc_point_obs.read_obs_data(buf_size2, i_start, (float *)obs_arr_block, + obs_qty_block, (char *)0); + if (!status) exit(1); int typ_idx, sid_idx, vld_idx; - for(int i_offset=0; i_offsetget_header(h, hdr_arr, hdr_typ_str, hdr_sid_str, hdr_vld_str); + else +#endif + nc_point_obs.get_header(h, hdr_arr, hdr_typ_str, hdr_sid_str, hdr_vld_str); // Store data in an observation object Observation cur_obs( @@ -263,12 +299,15 @@ void process_point_obs(const char *point_obs_filename) { } // end for i_start // Clean up +#ifdef WITH_PYTHON + if (use_python) met_point_file.close(); + else +#endif nc_point_obs.close(); return; } - //////////////////////////////////////////////////////////////////////// void create_plot() { @@ -470,7 +509,7 @@ void create_plot() { if(use_flate) plot.end_flate(); mlog << Debug(2) - << "Finished plotting " << plot_count << " locations.\n" + << "Finished plotting " << plot_count << " locations for " << ps_file << ".\n" << "Skipped " << skip_count << " locations off the grid.\n"; plot.showpage(); diff --git a/met/src/tools/other/plot_point_obs/plot_point_obs.h b/met/src/tools/other/plot_point_obs/plot_point_obs.h index 7fec0e19f0..0bf33b0679 100644 --- a/met/src/tools/other/plot_point_obs/plot_point_obs.h +++ b/met/src/tools/other/plot_point_obs/plot_point_obs.h @@ -43,12 +43,13 @@ using namespace netCDF; #include "plot_point_obs_conf_info.h" -#include "vx_data2d_factory.h" -#include "vx_grid.h" #include "vx_util.h" #include "vx_stat_out.h" #include "vx_gsl_prob.h" #include "nc_utils.h" +#ifdef WITH_PYTHON +#include "pointdata_python.h" +#endif //////////////////////////////////////////////////////////////////////// // diff --git a/met/src/tools/other/point2grid/point2grid.cc b/met/src/tools/other/point2grid/point2grid.cc index 7b840f08d9..f806023c11 100644 --- a/met/src/tools/other/point2grid/point2grid.cc +++ b/met/src/tools/other/point2grid/point2grid.cc @@ -43,6 +43,11 @@ using namespace std; #include "point2grid_conf_info.h" +#ifdef WITH_PYTHON +#include "data2d_nc_met.h" +#include "pointdata_python.h" +#endif + //////////////////////////////////////////////////////////////////////// static ConcatString program_name; @@ -53,6 +58,7 @@ static const int TYPE_OBS = 1; // MET Point Obs NetCDF (from xxx2nc) static const int TYPE_NCCF = 2; // CF NetCDF with time and lat/lon variables static const int TYPE_GOES = 5; static const int TYPE_GOES_ADP = 6; +static const int TYPE_PYTHON = 7; // MET Point Obs NetCDF from PYTHON static const InterpMthd DefaultInterpMthd = InterpMthd_UW_Mean; static const int DefaultInterpWdth = 2; @@ -120,14 +126,19 @@ static NcDim lon_dim ; static void process_command_line(int, char **); static void process_data_file(); static void process_point_file(NcFile *nc_in, MetConfig &config, - VarInfo *, const Grid fr_grid, const Grid to_grid); + VarInfo *, const Grid to_grid); +#ifdef WITH_PYTHON +static void process_point_python(string python_command, MetConfig &config, + VarInfo *vinfo, const Grid to_grid, bool use_xarray); +#endif static void process_point_nccf_file(NcFile *nc_in, MetConfig &config, - VarInfo *, Met2dDataFile *fr_mtddf, const Grid to_grid); + VarInfo *, Met2dDataFile *fr_mtddf, + const Grid to_grid); static void open_nc(const Grid &grid, const ConcatString run_cs); static void write_nc(const DataPlane &dp, const Grid &grid, const VarInfo *vinfo, const char *vname); static void write_nc_int(const DataPlane &dp, const Grid &grid, - const VarInfo *vinfo, const char *vname); + const VarInfo *vinfo, const char *vname); static void close_nc(); static void usage(); static void set_field(const StringArray &); @@ -265,11 +276,26 @@ void process_command_line(int argc, char **argv) { OutputFilename = cline[2]; // Check if the input file - if ( !file_exists(InputFilename.c_str()) ) { - mlog << Error << "\n" << method_name - << "file does not exist \"" << InputFilename << "\"\n\n"; - exit(1); + bool use_python = false; +#ifdef WITH_PYTHON + string python_command = InputFilename; + bool use_xarray = (0 == python_command.find(conf_val_python_xarray)); + use_python = use_xarray || (0 == python_command.find(conf_val_python_numpy)); + if (use_python) { + int offset = python_command.find("="); + if (offset == std::string::npos) { + mlog << Error << "\n" << method_name + << "trouble parsing the python command " << python_command << ".\n\n"; + exit(1); + } } + else +#endif + if ( !file_exists(InputFilename.c_str()) ) { + mlog << Error << "\n" << method_name + << "file does not exist \"" << InputFilename << "\"\n\n"; + exit(1); + } // Check for at least one configuration string if(FieldSA.n() < 1) { @@ -347,18 +373,42 @@ void process_data_file() { // Open the input file mlog << Debug(1) << "Reading data file: " << InputFilename << "\n"; - nc_in = open_ncfile(InputFilename.c_str()); - - // Get the obs type before opening NetCDF - int obs_type = get_obs_type(nc_in); - bool goes_data = (obs_type == TYPE_GOES || obs_type == TYPE_GOES_ADP); - - if (obs_type == TYPE_NCCF) setenv(nc_att_met_point_nccf, "yes", 1); - - // Read the input data file + bool goes_data = false; + bool use_python = false; + int obs_type = TYPE_UNKNOWN; Met2dDataFileFactory m_factory; Met2dDataFile *fr_mtddf = (Met2dDataFile *) 0; - fr_mtddf = m_factory.new_met_2d_data_file(InputFilename.c_str(), ftype); +#ifdef WITH_PYTHON + string python_command = InputFilename; + MetPointData *met_point_obs = 0; + bool use_xarray = (0 == python_command.find(conf_val_python_xarray)); + use_python = use_xarray || (0 == python_command.find(conf_val_python_numpy)); + if (use_python) { + int offset = python_command.find("="); + if (offset == std::string::npos) { + mlog << Error << "\n" << method_name + << "trouble parsing the python command " << python_command << ".\n\n"; + exit(1); + } + + python_command = python_command.substr(offset+1); + obs_type = TYPE_PYTHON; + fr_mtddf = new MetNcMetDataFile(); + } + else +#endif + { + nc_in = open_ncfile(InputFilename.c_str()); + + // Get the obs type before opening NetCDF + obs_type = get_obs_type(nc_in); + goes_data = (obs_type == TYPE_GOES || obs_type == TYPE_GOES_ADP); + + if (obs_type == TYPE_NCCF) setenv(nc_att_met_point_nccf, "yes", 1); + + // Read the input data file + fr_mtddf = m_factory.new_met_2d_data_file(InputFilename.c_str(), ftype); + } if(!fr_mtddf) { mlog << Error << "\n" << method_name @@ -384,7 +434,11 @@ void process_data_file() { // For python types read the first field to set the grid // Determine the "from" grid +#ifdef WITH_PYTHON + if (!use_python) fr_grid = fr_mtddf->grid(); +#else fr_grid = fr_mtddf->grid(); +#endif // Determine the "to" grid to_grid = parse_vx_grid(RGInfo, &fr_grid, &fr_grid); @@ -410,12 +464,17 @@ void process_data_file() { process_goes_file(nc_in, config, vinfo, fr_grid, to_grid); } else if (TYPE_OBS == obs_type) { - process_point_file(nc_in, config, vinfo, fr_grid, to_grid); + process_point_file(nc_in, config, vinfo, to_grid); } else if (TYPE_NCCF == obs_type) { process_point_nccf_file(nc_in, config, vinfo, fr_mtddf, to_grid); unsetenv(nc_att_met_point_nccf); } +#ifdef WITH_PYTHON + else if (TYPE_PYTHON == obs_type) { + process_point_python(python_command, config, vinfo, to_grid, use_xarray); + } +#endif else { mlog << Error << "\n" << method_name << "Please check the input file. Only supports GOES, MET point obs, " @@ -599,8 +658,8 @@ IntArray prepare_qc_array(const IntArray qc_flags, StringArray qc_tables) { //////////////////////////////////////////////////////////////////////// -void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, - const Grid fr_grid, const Grid to_grid) { +void process_point_met_data(MetPointData *met_point_obs, MetConfig &config, VarInfo *vinfo, + const Grid to_grid) { int nhdr, nobs; int nx, ny, var_count, to_count, var_count2; int idx, hdr_idx; @@ -616,8 +675,8 @@ void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, bool has_prob_thresh = !prob_cat_thresh.check(bad_data_double); unixtime requested_valid_time, valid_time; - static const char *method_name = "process_point_file() -> "; - static const char *method_name_s = "process_point_file()"; + static const char *method_name = "process_point_met_data() -> "; + static const char *method_name_s = "process_point_met_data()"; // Check for at least one configuration string if(FieldSA.n() < 1) { @@ -626,35 +685,29 @@ void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, usage(); } - MetNcPointObsIn nc_point_obs; - nc_point_obs.set_netcdf(nc_in, true); - // Read the dimensions and variables - nc_point_obs.read_dim_headers(); - nc_point_obs.check_nc(GET_NC_NAME_P(nc_in).c_str(), method_name_s); // exit if missing dims/vars - // Read all obs data to compute the cell mapping - nc_point_obs.read_obs_data(); - NcHeaderData header_data = nc_point_obs.get_header_data(); - NcPointObsData obs_data = nc_point_obs.get_point_obs_data(); +// MetNcPointObsIn nc_point_obs; + MetPointHeader *header_data = met_point_obs->get_header_data(); + MetPointObsData *obs_data = met_point_obs->get_point_obs_data(); - nhdr = nc_point_obs.get_hdr_cnt(); - nobs = nc_point_obs.get_obs_cnt(); + nhdr = met_point_obs->get_hdr_cnt(); + nobs = met_point_obs->get_obs_cnt(); bool empty_input = (nhdr == 0 && nobs == 0); - bool use_var_id = nc_point_obs.is_using_var_id(); + bool use_var_id = met_point_obs->is_using_var_id(); float *hdr_lats = new float[nhdr]; float *hdr_lons = new float[nhdr]; IntArray var_index_array; IntArray valid_time_array; - StringArray qc_tables = nc_point_obs.get_qty_data(); - StringArray var_names = nc_point_obs.get_var_names(); - StringArray hdr_valid_times = header_data.vld_array; + StringArray qc_tables = met_point_obs->get_qty_data(); + StringArray var_names = met_point_obs->get_var_names(); + StringArray hdr_valid_times = header_data->vld_array; hdr_valid_times.sort(); - nc_point_obs.get_lats(hdr_lats); - nc_point_obs.get_lons(hdr_lons); + met_point_obs->get_lats(hdr_lats); + met_point_obs->get_lons(hdr_lons); // Check the message types - prepare_message_types(header_data.typ_array); + prepare_message_types(header_data->typ_array); // Check and read obs_vid and obs_var if exists bool success_to_read = true; @@ -732,7 +785,7 @@ void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, else { bool not_found_grib_code = true; for (idx=0; idxobs_ids[idx]) { not_found_grib_code = false; break; } @@ -757,10 +810,10 @@ void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, log_msg << "GRIB codes: "; IntArray grib_codes; for (idx=0; idxobs_ids[idx])) { + grib_codes.add(obs_data->obs_ids[idx]); if (0 < idx) log_msg << ", "; - log_msg << obs_data.obs_ids[idx]; + log_msg << obs_data->obs_ids[idx]; } } } @@ -836,29 +889,29 @@ void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, var_count = var_count2 = to_count = 0; filtered_by_time = filtered_by_msg_type = filtered_by_qc = 0; for (idx=0; idx < nobs; idx++) { - if (var_idx_or_gc == obs_data.obs_ids[idx]) { + if (var_idx_or_gc == obs_data->obs_ids[idx]) { var_count2++; - hdr_idx = obs_data.obs_hids[idx]; + hdr_idx = obs_data->obs_hids[idx]; if (0 < valid_time_array.n() && - !valid_time_array.has(header_data.vld_idx_array[hdr_idx])) { + !valid_time_array.has(header_data->vld_idx_array[hdr_idx])) { filtered_by_time++; continue; } - if(!keep_message_type(header_data.typ_idx_array[hdr_idx])) { + if(!keep_message_type(header_data->typ_idx_array[hdr_idx])) { filtered_by_msg_type++; continue; } // Filter by QC flag - if (has_qc_flags && !qc_idx_array.has(obs_data.obs_qids[idx])) { + if (has_qc_flags && !qc_idx_array.has(obs_data->obs_qids[idx])) { filtered_by_qc++; continue; } var_index_array.add(idx); var_count++; - if (is_eq(obs_data.obs_vals[idx], 0.)) obs_count_zero_from++; + if (is_eq(obs_data->obs_vals[idx], 0.)) obs_count_zero_from++; else obs_count_non_zero_from++; } } @@ -869,7 +922,7 @@ void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, } cellMapping = new IntArray[nx * ny]; if( get_grid_mapping(to_grid, cellMapping, var_index_array, - obs_data.obs_hids, hdr_lats, hdr_lons) ) { + obs_data->obs_hids, hdr_lats, hdr_lons) ) { int from_index; IntArray cellArray; NumArray dataArray; @@ -905,7 +958,7 @@ void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, dataArray.extend(cellArray.n()); for (int dIdx=0; dIdxget_obs_val(from_index); if (is_bad_data(data_value)) continue; if(mlog.verbosity_level() >= 4) { @@ -1073,13 +1126,108 @@ void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, delete [] hdr_lats; delete [] hdr_lons; + + return; +} + +//////////////////////////////////////////////////////////////////////// + +void process_point_file(NcFile *nc_in, MetConfig &config, VarInfo *vinfo, + const Grid to_grid) { + int nhdr, nobs; + int nx, ny, var_count, to_count, var_count2; + int idx, hdr_idx; + int var_idx_or_gc; + int filtered_by_time, filtered_by_msg_type, filtered_by_qc; + ConcatString vname, vname_cnt, vname_mask; + DataPlane fr_dp, to_dp; + DataPlane cnt_dp, mask_dp; + DataPlane prob_dp, prob_mask_dp; + NcVar var_obs_gc, var_obs_var; + + clock_t start_clock = clock(); + bool has_prob_thresh = !prob_cat_thresh.check(bad_data_double); + + unixtime requested_valid_time, valid_time; + static const char *method_name = "process_point_file() -> "; + static const char *method_name_s = "process_point_file()"; + + // Check for at least one configuration string + if(FieldSA.n() < 1) { + mlog << Error << "\n" << method_name + << "The -field option must be used at least once!\n\n"; + usage(); + } + + MetNcPointObsIn nc_point_obs; + nc_point_obs.set_netcdf(nc_in, true); + // Read the dimensions and variables + nc_point_obs.read_dim_headers(); + nc_point_obs.check_nc(GET_NC_NAME_P(nc_in).c_str(), method_name_s); // exit if missing dims/vars + // Read all obs data to compute the cell mapping + nc_point_obs.read_obs_data(); + process_point_met_data(&nc_point_obs, config, vinfo, to_grid); + nc_point_obs.close(); mlog << Debug(LEVEL_FOR_PERFORMANCE) << method_name << "took " << (clock()-start_clock)/double(CLOCKS_PER_SEC) << " seconds\n"; +} + +//////////////////////////////////////////////////////////////////////// + +#ifdef WITH_PYTHON + +void process_point_python(string python_command, MetConfig &config, VarInfo *vinfo, + const Grid to_grid, bool use_xarray) { + int nhdr, nobs; + int nx, ny, var_count, to_count, var_count2; + int idx, hdr_idx; + int var_idx_or_gc; + int filtered_by_time, filtered_by_msg_type, filtered_by_qc; + ConcatString vname, vname_cnt, vname_mask; + DataPlane fr_dp, to_dp; + DataPlane cnt_dp, mask_dp; + DataPlane prob_dp, prob_mask_dp; + NcVar var_obs_gc, var_obs_var; + + clock_t start_clock = clock(); + bool has_prob_thresh = !prob_cat_thresh.check(bad_data_double); + + unixtime requested_valid_time, valid_time; + static const char *method_name = "process_point_python() -> "; + static const char *method_name_s = "process_point_python()"; + + // Check for at least one configuration string + if(FieldSA.n() < 1) { + mlog << Error << "\n" << method_name + << "The -field option must be used at least once!\n\n"; + usage(); + } + + MetPythonPointDataFile met_point_file; + if(!met_point_file.open(python_command.c_str(), use_xarray)) { + met_point_file.close(); + + mlog << Error << "\n" << method_name + << "trouble getting point observation file from python command " + << python_command << "\n\n"; + + exit(1); + } + + MetPointData *met_point_obs = met_point_file.get_met_point_data(); + process_point_met_data(met_point_obs, config, vinfo, to_grid); + + met_point_file.close(); + + mlog << Debug(LEVEL_FOR_PERFORMANCE) << method_name << "took " + << (clock()-start_clock)/double(CLOCKS_PER_SEC) << " seconds\n"; + return; } +#endif //////////////////////////////////////////////////////////////////////// @@ -1097,7 +1245,7 @@ void process_point_nccf_file(NcFile *nc_in, MetConfig &config, bool opt_all_attrs = false; Grid fr_grid = fr_mtddf->grid(); int from_size = fr_grid.nx() * fr_grid.ny(); - static const char *method_name = "process_point_file_with_latlon() -> "; + static const char *method_name = "process_point_nccf_file() -> "; NcVar var_lat = get_nc_var_lat(nc_in); NcVar var_lon = get_nc_var_lon(nc_in); diff --git a/test/xml/unit_python.xml b/test/xml/unit_python.xml index 698f550f69..226d61869b 100644 --- a/test/xml/unit_python.xml +++ b/test/xml/unit_python.xml @@ -448,4 +448,111 @@ + + &MET_BIN;/point2grid + \ + 'PYTHON_NUMPY=&MET_BASE;/python/read_met_point_obs.py &OUTPUT_DIR;/pb2nc/ndas.20120409.t12z.prepbufr.tm00.nc' \ + G212 \ + &OUTPUT_DIR;/python/pb2nc_TMP.nc \ + -field 'name="TMP"; level="*"; valid_time="20120409_120000"; censor_thresh=[ <0 ]; censor_val=[0];' \ + -name TEMP \ + -v 1 + + + &OUTPUT_DIR;/python/pb2nc_TMP.nc + + + + + &MET_BIN;/plot_point_obs + + TO_GRID NONE + + \ + &OUTPUT_DIR;/pb2nc/ndas.20120409.t12z.prepbufr.tm00.nc \ + &OUTPUT_DIR;/python/nam_and_ndas.20120409.t12z.prepbufr_CONFIG.ps \ + -point_obs &OUTPUT_DIR;/ascii2nc/trmm_2012040912_3hr.nc \ + -plot_grid &DATA_DIR_MODEL;/grib2/nam/nam_2012040900_F012.grib2 \ + -config &CONFIG_DIR;/PlotPointObsConfig \ + -title "NAM 2012040900 F12 vs NDAS 500mb RH and TRMM 3h > 0" \ + -v 3 + + + &OUTPUT_DIR;/python/nam_and_ndas.20120409.t12z.prepbufr_CONFIG.ps + + + + + &MET_BIN;/plot_point_obs + + TO_GRID NONE + + \ + 'PYTHON_NUMPY=&MET_BASE;/python/read_met_point_obs.py &OUTPUT_DIR;/pb2nc/ndas.20120409.t12z.prepbufr.tm00.nc' \ + &OUTPUT_DIR;/python/nam_and_ndas.20120409.t12z.prepbufr_CONFIG.ps \ + -point_obs 'PYTHON_NUMPY=&MET_BASE;/python/read_met_point_obs.py &OUTPUT_DIR;/ascii2nc/trmm_2012040912_3hr.nc' \ + -plot_grid &DATA_DIR_MODEL;/grib2/nam/nam_2012040900_F012.grib2 \ + -config &CONFIG_DIR;/PlotPointObsConfig \ + -title "NAM 2012040900 F12 vs NDAS 500mb RH and TRMM 3h > 0" \ + -v 3 + + + &OUTPUT_DIR;/python/nam_and_ndas.20120409.t12z.prepbufr_CONFIG.ps + + + + + echo "&DATA_DIR_MODEL;/grib1/arw-fer-gep1/arw-fer-gep1_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-fer-gep5/arw-fer-gep5_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-sch-gep2/arw-sch-gep2_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-sch-gep6/arw-sch-gep6_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-tom-gep3/arw-tom-gep3_2012040912_F024.grib \ + &DATA_DIR_MODEL;/grib1/arw-tom-gep7/arw-tom-gep7_2012040912_F024.grib" \ + > &OUTPUT_DIR;/python/ensemble_stat/input_file_list; \ + &MET_BIN;/ensemble_stat + + DESC NA + OBS_ERROR_FLAG FALSE + SKIP_CONST FALSE + OUTPUT_PREFIX FILE_LIST + + \ + &OUTPUT_DIR;/python/ensemble_stat/input_file_list \ + &CONFIG_DIR;/EnsembleStatConfig \ + -grid_obs &DATA_DIR_OBS;/laps/laps_2012041012_F000.grib \ + -point_obs 'PYTHON_NUMPY=&MET_BASE;/python/read_met_point_obs.py &OUTPUT_DIR;/ascii2nc/gauge_2012041012_24hr.nc' \ + -outdir &OUTPUT_DIR;/python/ensemble_stat -v 1 + + + &OUTPUT_DIR;/python/ensemble_stat/ensemble_stat_FILE_LIST_20120410_120000V.stat + &OUTPUT_DIR;/python/ensemble_stat/ensemble_stat_FILE_LIST_20120410_120000V_ecnt.txt + &OUTPUT_DIR;/python/ensemble_stat/ensemble_stat_FILE_LIST_20120410_120000V_rhist.txt + &OUTPUT_DIR;/python/ensemble_stat/ensemble_stat_FILE_LIST_20120410_120000V_phist.txt + &OUTPUT_DIR;/python/ensemble_stat/ensemble_stat_FILE_LIST_20120410_120000V_orank.txt + &OUTPUT_DIR;/python/ensemble_stat/ensemble_stat_FILE_LIST_20120410_120000V_ssvar.txt + &OUTPUT_DIR;/python/ensemble_stat/ensemble_stat_FILE_LIST_20120410_120000V_ens.nc + &OUTPUT_DIR;/python/ensemble_stat/ensemble_stat_FILE_LIST_20120410_120000V_orank.nc + + + + + &MET_BIN;/point_stat + + BEG_DS -1800 + END_DS 1800 + OUTPUT_PREFIX GRIB1_NAM_GDAS_WINDS + CONFIG_DIR &CONFIG_DIR; + CLIMO_FILE "&DATA_DIR_MODEL;/grib1/gfs/gfs_2012040900_F012_gNam.grib" + + \ + &DATA_DIR_MODEL;/grib1/nam/nam_2012040900_F012.grib \ + 'PYTHON_NUMPY=&MET_BASE;/python/read_met_point_obs.py &OUTPUT_DIR;/pb2nc/gdas1.20120409.t12z.prepbufr.nc' \ + &CONFIG_DIR;/PointStatConfig_WINDS \ + -outdir &OUTPUT_DIR;/python -v 1 + + + &OUTPUT_DIR;/python/point_stat_GRIB1_NAM_GDAS_WINDS_120000L_20120409_120000V.stat + + +