From 51aed288a3ca6cf60a617b6dfacd422da956e613 Mon Sep 17 00:00:00 2001 From: Howard Soh Date: Wed, 22 Feb 2023 16:44:30 -0700 Subject: [PATCH] #2428 Initial release --- data/wrappers/read_tmp_point_nc.py | 30 +++ data/wrappers/write_tmp_point_nc.py | 59 ++++++ scripts/python/met_point_obs_nc.py | 276 ++++++++++++++++++++++++++++ 3 files changed, 365 insertions(+) create mode 100644 data/wrappers/read_tmp_point_nc.py create mode 100644 data/wrappers/write_tmp_point_nc.py create mode 100644 scripts/python/met_point_obs_nc.py diff --git a/data/wrappers/read_tmp_point_nc.py b/data/wrappers/read_tmp_point_nc.py new file mode 100644 index 0000000000..8afaf16684 --- /dev/null +++ b/data/wrappers/read_tmp_point_nc.py @@ -0,0 +1,30 @@ +######################################################################## +# +# Reads temporary file into memory. +# +# usage: /path/to/python read_tmp_dataplane.py dataplane.tmp +# +######################################################################## + +import os +import sys + +met_base_dir = os.getenv('MET_BASE',None) +if met_base_dir is not None: + sys.path.append(os.path.join(met_base_dir, 'python')) + +# add share/met/python directory to system path to find met_point_obs +sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), + os.pardir, 'python'))) +from met_point_obs import met_point_obs +from met_point_obs_nc import nc_point_obs + +netcdf_filename = sys.argv[1] + +# read NetCDF file +print('{p} reading{f}'.format(p=met_point_obs.get_prompt(), f=netcdf_filename)) +point_obs_data = nc_point_obs() +point_obs_data.read_data(netcdf_filename) + +met_point_data = point_obs_data.get_point_data() +met_point_data['met_point_data'] = point_obs_data diff --git a/data/wrappers/write_tmp_point_nc.py b/data/wrappers/write_tmp_point_nc.py new file mode 100644 index 0000000000..1dbda19807 --- /dev/null +++ b/data/wrappers/write_tmp_point_nc.py @@ -0,0 +1,59 @@ +######################################################################## +# +# Adapted from a script provided by George McCabe +# Adapted by Randy Bullock +# +# usage: /path/to/python write_tmp_point.py \ +# tmp_output_filename .py +# +######################################################################## + +import os +import sys +import importlib.util + +met_base_dir = os.getenv('MET_BASE',None) +if met_base_dir is not None: + sys.path.append(os.path.join(met_base_dir, 'python')) + +# add share/met/python directory to system path to find met_point_obs +sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), + os.pardir, 'python'))) + +from met_point_obs import met_point_obs +from met_point_obs_nc import nc_point_obs + +PROMPT = met_point_obs.get_prompt() +print("{p} Python Script:\t".format(p=PROMPT) + repr(sys.argv[0])) +print("{p} User Command:\t".format(p=PROMPT) + repr(' '.join(sys.argv[2:]))) +print("{p} Temporary File:\t".format(p=PROMPT) + repr(sys.argv[1])) + +tmp_filename = sys.argv[1] +pyembed_module_name = sys.argv[2] +sys.argv = sys.argv[2:] + +# append user script dir to system path +pyembed_dir, pyembed_file = os.path.split(pyembed_module_name) +if pyembed_dir: + sys.path.insert(0, pyembed_dir) + +if not pyembed_module_name.endswith('.py'): + pyembed_module_name += '.py' + +user_base = os.path.basename(pyembed_module_name).replace('.py','') + +spec = importlib.util.spec_from_file_location(user_base, pyembed_module_name) +met_in = importlib.util.module_from_spec(spec) +spec.loader.exec_module(met_in) + +if hasattr(met_in, 'point_obs_data'): + met_in.point_obs_data.save_ncfile(tmp_filename) +else: + if hasattr(met_in.met_point_data, 'point_obs_data'): + met_in.met_point_data['point_obs_data'].save_ncfile(tmp_filename) + else: + tmp_point_obs = nc_point_obs() + tmp_point_obs.put_data(met_in.met_point_data) + tmp_point_obs.save_ncfile(tmp_filename) + +#print('{p} writing {f}'.format(p=PROMPT, f=tmp_filename)) diff --git a/scripts/python/met_point_obs_nc.py b/scripts/python/met_point_obs_nc.py new file mode 100644 index 0000000000..cb86be5011 --- /dev/null +++ b/scripts/python/met_point_obs_nc.py @@ -0,0 +1,276 @@ +#!/usr/bin/env python3 + +''' +Separated from read_met_point_obs on Feb 09, 2023 + +@author: hsoh + +This script reads the MET point observation NetCDF file like MET tools do. +''' + +import os +import sys +from datetime import datetime +import numpy as np +import netCDF4 as nc + +from met_point_obs import met_point_obs, base_met_point_obs, get_prompt + +DO_PRINT_DATA = False +ARG_PRINT_DATA = 'show_data' + +# Note: caller should import netCDF4 +# the argements nc_group(dataset) and nc_var should not be None +class nc_tools(): + + met_missing = -99999999. + + @staticmethod + def get_num_array(nc_group, var_name): + nc_var = nc_group.variables.get(var_name, None) + return [] if nc_var is None else nc_var[:] + + @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(nc_group, var_name): + nc_var = nc_group.variables.get(var_name, None) + return [] if nc_var is None else nc_tools.get_ncbyte_array_to_str(nc_var) + + +class nc_point_obs(met_point_obs): + + # args should be string, list, or dictionary + def get_nc_filename(self, args): + nc_filename = None + if isinstance(args, dict): + nc_filename = args.get('nc_name',None) + elif isinstance(args, list): + nc_filename = args[0] + elif args != ARG_PRINT_DATA: + nc_filename = args + + return nc_filename + + def read_data(self, nc_filename): + if nc_filename is None: + self.log_error_msg("The input NetCDF filename is missing") + elif not os.path.exists(nc_filename): + self.log_error_msg(f"input NetCDF file ({nc_filename}) does not exist") + 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('obs_unit', None) + if nc_var: + self.obs_var_unit = nc_var[:] + nc_var = dataset.variables.get('obs_desc', None) + if nc_var: + self.obs_var_desc = nc_var[:] + + 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') + + def save_ncfile(self, nc_filename): + met_data = self.get_point_data() + with nc.Dataset(nc_filename, 'w') as nc_dataset: + self.set_nc_data(nc_dataset) + return met_data + + def set_nc_data(self, nc_dataset): + return nc_point_obs.write_nc_data(nc_dataset, self) + + @staticmethod + def write_nc_file(nc_filename, point_obs): + with nc.Dataset(nc_filename, 'w') as nc_dataset: + nc_point_obs.set_nc_data(nc_dataset, point_obs) + + @staticmethod + def write_nc_data(nc_dataset, point_obs): + do_nothing = False + if 0 == point_obs.nhdr: + do_nothing = True + base_met_point_obs.info_msg("the header is empty") + if 0 == point_obs.nobs: + do_nothing = True + base_met_point_obs.info_msg("the observation data is empty") + if do_nothing: + print() + return + + # Set global attributes + nc_dataset.MET_Obs_version = "1.02" ; + nc_dataset.use_var_id = "true" if point_obs.use_var_id else "false" + + # Create dimensions + nc_dataset.createDimension('mxstr', 16) + nc_dataset.createDimension('mxstr2', 40) + nc_dataset.createDimension('mxstr3', 80) + nc_dataset.createDimension('nhdr', point_obs.nhdr) + nc_dataset.createDimension('nobs', point_obs.nobs) + #npbhdr = len(point_obs.hdr_prpt_typ) + if 0 < point_obs.npbhdr: + nc_dataset.createDimension('npbhdr', point_obs.npbhdr) + nc_dataset.createDimension('nhdr_typ', point_obs.nhdr_typ) + nc_dataset.createDimension('nhdr_sid', point_obs.nhdr_sid) + nc_dataset.createDimension('nhdr_vld', point_obs.nhdr_vld) + nc_dataset.createDimension('nobs_qty', point_obs.nobs_qty) + nc_dataset.createDimension('obs_var_num', point_obs.nobs_var) + + type_for_string = 'S1' # np.byte + dims_hdr = ('nhdr',) + dims_obs = ('nobs',) + + # Create header and observation variables + var_hdr_typ = nc_dataset.createVariable('hdr_typ', np.int32, dims_hdr, fill_value=-9999) + var_hdr_sid = nc_dataset.createVariable('hdr_sid', np.int32, dims_hdr, fill_value=-9999) + var_hdr_vld = nc_dataset.createVariable('hdr_vld', np.int32, dims_hdr, fill_value=-9999) + var_hdr_lat = nc_dataset.createVariable('hdr_lat', np.float32, dims_hdr, fill_value=-9999.) + var_hdr_lon = nc_dataset.createVariable('hdr_lon', np.float32, dims_hdr, fill_value=-9999.) + var_hdr_elv = nc_dataset.createVariable('hdr_elv', np.float32, dims_hdr, fill_value=-9999.) + + var_obs_qty = nc_dataset.createVariable('obs_qty', np.int32, dims_obs, fill_value=-9999) + var_obs_hid = nc_dataset.createVariable('obs_hid', np.int32, dims_obs, fill_value=-9999) + var_obs_vid = nc_dataset.createVariable('obs_vid', np.int32, dims_obs, fill_value=-9999) + var_obs_lvl = nc_dataset.createVariable('obs_lvl', np.float32, dims_obs, fill_value=-9999.) + var_obs_hgt = nc_dataset.createVariable('obs_hgt', np.float32, dims_obs, fill_value=-9999.) + var_obs_val = nc_dataset.createVariable('obs_val', np.float32, dims_obs, fill_value=-9999.) + + if 0 == point_obs.npbhdr: + var_hdr_prpt_typ = None + var_hdr_irpt_typ = None + var_hdr_inst_typ = None + else: + dims_npbhdr = ('npbhdr',) + var_hdr_prpt_typ = nc_dataset.createVariable('hdr_prpt_typ', np.int32, dims_npbhdr, fill_value=-9999.) + var_hdr_irpt_typ = nc_dataset.createVariable('hdr_irpt_typ', np.int32, dims_npbhdr, fill_value=-9999.) + var_hdr_inst_typ = nc_dataset.createVariable('hdr_inst_typ', np.int32, dims_npbhdr, fill_value=-9999.) + + var_hdr_typ_table = nc_dataset.createVariable('hdr_typ_table', type_for_string, ('nhdr_typ','mxstr2')) + var_hdr_sid_table = nc_dataset.createVariable('hdr_sid_table', type_for_string, ('nhdr_sid','mxstr2')) + var_hdr_vld_table = nc_dataset.createVariable('hdr_vld_table', type_for_string, ('nhdr_vld','mxstr')) + var_obs_qty_table = nc_dataset.createVariable('obs_qty_table', type_for_string, ('nobs_qty','mxstr')) + var_obs_var_table = nc_dataset.createVariable('obs_var', type_for_string, ('obs_var_num','mxstr2')) + var_obs_var_unit = nc_dataset.createVariable('obs_unit', type_for_string, ('obs_var_num','mxstr2')) + var_obs_var_desc = nc_dataset.createVariable('obs_desc', type_for_string, ('obs_var_num','mxstr3')) + + # Set variables + var_hdr_typ[:] = point_obs.hdr_typ[:] + var_hdr_sid[:] = point_obs.hdr_sid[:] + var_hdr_vld[:] = point_obs.hdr_vld[:] + var_hdr_lat[:] = point_obs.hdr_lat[:] + var_hdr_lon[:] = point_obs.hdr_lon[:] + var_hdr_elv[:] = point_obs.hdr_elv[:] + for i in range(0, point_obs.nhdr_typ): + for j in range(0, len(point_obs.hdr_typ_table[i])): + var_hdr_typ_table[i,j] = point_obs.hdr_typ_table[i][j] + for i in range(0, point_obs.nhdr_sid): + for j in range(0, len(point_obs.hdr_sid_table[i])): + var_hdr_sid_table[i,j] = point_obs.hdr_sid_table[i][j] + for i in range(0, point_obs.nhdr_vld): + for j in range(0, len(point_obs.hdr_vld_table[i])): + var_hdr_vld_table[i,j] = point_obs.hdr_vld_table[i][j] + if 0 < point_obs.npbhdr: + var_hdr_prpt_typ[:] = point_obs.hdr_prpt_typ[:] + var_hdr_irpt_typ[:] = point_obs.hdr_irpt_typ[:] + var_hdr_inst_typ[:] = point_obs.hdr_inst_typ[:] + + var_obs_qty[:] = point_obs.obs_qty[:] + var_obs_hid[:] = point_obs.obs_hid[:] + var_obs_vid[:] = point_obs.obs_vid[:] + var_obs_lvl[:] = point_obs.obs_lvl[:] + var_obs_hgt[:] = point_obs.obs_hgt[:] + var_obs_val[:] = point_obs.obs_val[:] + for i in range(0, point_obs.nobs_var): + for j in range(0, len(point_obs.obs_var_table[i])): + var_obs_var_table[i,j] = point_obs.obs_var_table[i][j] + var_obs_var_unit[i] = "" if i >= len(point_obs.obs_var_unit) else point_obs.obs_var_unit[i] + var_obs_var_desc[i] = "" if i >= len(point_obs.obs_var_desc) else point_obs.obs_var_desc[i] + for i in range(0, point_obs.nobs_qty): + for j in range(0, len(point_obs.obs_qty_table[i])): + var_obs_qty_table[i,j] = point_obs.obs_qty_table[i][j] + + # Set variable attributes + var_hdr_typ.long_name = "index of message type" + var_hdr_sid.long_name = "index of station identification" + var_hdr_vld.long_name = "index of valid time" + var_hdr_lat.long_name = "latitude" + var_hdr_lat.units = "degrees_north" + var_hdr_lon.long_name = "longitude" + var_hdr_lon.units = "degrees_east" + var_hdr_elv.long_name = "elevation" + var_hdr_elv.units = "meters above sea level (msl)" + + var_obs_qty.long_name = "index of quality flag" + var_obs_hid.long_name = "index of matching header data" + var_obs_vid.long_name = "index of BUFR variable corresponding to the observation type" + var_obs_lvl.long_name = "pressure level (hPa) or accumulation interval (sec)" + var_obs_hgt.long_name = "height in meters above sea level (msl)" + var_obs_val.long_name = "observation value" + var_hdr_typ_table.long_name = "message type" + var_hdr_sid_table.long_name = "station identification" + var_hdr_vld_table.long_name = "valid time" + var_hdr_vld_table.units = "YYYYMMDD_HHMMSS UTC" + var_obs_qty_table.long_name = "quality flag" + var_obs_var_table.long_name = "variable names" + var_obs_var_unit.long_name = "variable units" + var_obs_var_desc.long_name = "variable descriptions" + + +def main(argv): + if len(argv) != 1 and argv[1] != ARG_PRINT_DATA: + netcdf_filename = argv[1] + tmp_nc_name = 'tmp_met_point.nc' + point_obs_data = nc_point_obs() + point_obs_data.read_data(point_obs_data.get_nc_filename(netcdf_filename)) + met_point_data = point_obs_data.save_ncfile(tmp_nc_name) + print(f'{get_prompt()} saved met_point_data to {tmp_nc_name}') + met_point_data['met_point_data'] = point_obs_data + + if DO_PRINT_DATA or ARG_PRINT_DATA == argv[-1]: + met_point_obs.print_point_data(met_point_data) + +if __name__ == '__main__': + start_time = datetime.now() + main(sys.argv) + run_time = datetime.now() - start_time + print(f'{get_prompt()} Done python script {sys.argv[0]} took {run_time}')