# Copyright (C) 2015-2026 S[&]T, The Netherlands. # All rights reserved. # # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: # # 1. Redistributions of source code must retain the above copyright notice, # this list of conditions and the following disclaimer. # # 2. Redistributions in binary form must reproduce the above copyright # notice, this list of conditions and the following disclaimer in the # documentation and/or other materials provided with the distribution. # # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from # this software without specific prior written permission. # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF # SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS # INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN # CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. from __future__ import print_function from collections import OrderedDict import csv import datetime import glob import numpy import os import subprocess from functools import reduce try: from cStringIO import StringIO except ImportError: try: from StringIO import StringIO except ImportError: from io import StringIO from harp._harpc import ffi as _ffi __all__ = ["Error", "CLibraryError", "UnsupportedTypeError", "UnsupportedDimensionError", "NoDataError", "Variable", "Product", "get_encoding", "set_encoding", "version", "import_product", "import_product_metadata", "export_product", "concatenate", "execute_operations", "convert_unit", "to_dict"] class Error(Exception): """Exception base class for all HARP Python interface errors.""" pass class CLibraryError(Error): """Exception raised when an error occurs inside the HARP C library. Attributes: errno -- error code; if None, the error code will be retrieved from the HARP C library. strerror -- error message; if None, the error message will be retrieved from the HARP C library. """ def __init__(self, errno=None, strerror=None): if errno is None: errno = _lib.harp_errno if strerror is None: strerror = _decode_string(_ffi.string(_lib.harp_errno_to_string(errno))) super(CLibraryError, self).__init__(errno, strerror) self.errno = errno self.strerror = strerror def __str__(self): return self.strerror class UnsupportedTypeError(Error): """Exception raised when unsupported types are encountered, either on the Python or on the C side of the interface. """ pass class UnsupportedDimensionError(Error): """Exception raised when unsupported dimensions are encountered, either on the Python or on the C side of the interface. """ pass class NoDataError(Error): """Exception raised when the product returned from an import contains no variables, or variables without data. """ def __init__(self): super(NoDataError, self).__init__("product contains no variables, or variables without data") class Variable(object): """Python representation of a HARP variable. A variable consists of data (either a scalar or NumPy array), a list of dimension types that describe the dimensions of the data, and a number of optional attributes: physical unit, minimum valid value, maximum valid value, human-readable description, and enumeration name list. """ def __init__(self, data, dimension=[], unit=None, valid_min=None, valid_max=None, description=None, enum=None): self.data = data self.dimension = dimension if unit is not None: self.unit = unit if valid_min is not None: self.valid_min = valid_min if valid_max is not None: self.valid_max = valid_max if description is not None: self.description = description if enum is not None: self.enum = enum def _get_c_data_type(self): """Return the C data type code corresponding to this variable, based on the its data and optionally valid_min/max attributes. """ c_data_type = _get_c_data_type(self.data) dimension = getattr(self, "dimension", []) # TODO can we avoid getattr? if len(dimension) == 0: # Allow valid_min/valid_max to influence data type as well try: c_data_type = max(c_data_type, _get_c_data_type(self.valid_min)) except Exception: pass try: c_data_type = max(c_data_type, _get_c_data_type(self.valid_max)) except Exception: pass return c_data_type def _format_data_type(self): """Return the string representation of the C data type that would be used to store the variable, or "" if its data attribute is of an unsupported type. """ try: return _get_c_data_type_name(self._get_c_data_type()) except UnsupportedTypeError: return "" def __repr__(self): if not self.dimension: return "" % self._format_data_type() return "" % (self._format_data_type(), _format_dimensions(self.dimension, self.data)) def __str__(self): stream = StringIO() print("type =", self._format_data_type(), file=stream) if self.dimension: print("dimension =", _format_dimensions(self.dimension, self.data), file=stream) try: unit = self.unit except AttributeError: pass else: if unit is not None: print("unit = %r" % unit, file=stream) try: valid_min = self.valid_min except AttributeError: pass else: if valid_min is not None: print("valid_min = %r" % valid_min, file=stream) try: valid_max = self.valid_max except AttributeError: pass else: if valid_max is not None: print("valid_max = %r" % valid_max, file=stream) try: description = self.description except AttributeError: pass else: if description: print("description = %r" % description, file=stream) try: enum = self.enum except AttributeError: pass else: if enum: print("enum = %r" % enum, file=stream) if self.data is not None: if not isinstance(self.data, numpy.ndarray) and not numpy.isscalar(self.data): print("data = ", file=stream) elif numpy.isscalar(self.data): print("data = %r" % self.data, file=stream) elif not self.dimension and self.data.size == 1: print("data = %r" % self.data.flat[0], file=stream) elif self.data.size == 0: print("data = ", file=stream) else: print("data =", file=stream) print(str(self.data), file=stream) return stream.getvalue() class Product(object): """Python representation of a HARP product. A product consists of product attributes and variables. Any attribute of a Product instance of which the name does not start with an underscore is either a variable or a product attribute. Product attribute names are reserved and cannot be used for variables. The list of names reserved for product attributes is: source_product -- Name of the original product this product is derived from. history -- New-line separated list of invocations of HARP command line tools that have been performed on the product. """ # Product attribute names. All attribute names of this class that do not start with an underscore are assumed to be # HARP variable names, except for the reserved names listed below. _reserved_names = set(("source_product", "history")) def __init__(self, source_product="", history=""): if source_product: self.source_product = source_product if history: self.history = history self._variable_dict = OrderedDict() def _is_reserved_name(self, name): return name.startswith("_") or name in Product._reserved_names def _verify_key(self, key): if not isinstance(key, str): # The statement obj.__class__.__name__ works for both new-style and old-style classes. raise TypeError("key must be str, not %r" % key.__class__.__name__) if self._is_reserved_name(key): raise KeyError(key) def __setattr__(self, name, value): super(Product, self).__setattr__(name, value) if not self._is_reserved_name(name): self._variable_dict[str(name)] = value def __delattr__(self, name): super(Product, self).__delattr__(name) if not self._is_reserved_name(name): del self._variable_dict[str(name)] def __getitem__(self, key): self._verify_key(key) try: return getattr(self, key) except AttributeError: raise KeyError(key) def __setitem__(self, key, value): self._verify_key(key) setattr(self, key, value) def __delitem__(self, key): self._verify_key(key) try: delattr(self, key) except AttributeError: raise KeyError(key) def __len__(self): return len(self._variable_dict) def __iter__(self): return iter(self._variable_dict) def __reversed__(self): return reversed(self._variable_dict) def __contains__(self, name): return name in self._variable_dict def __repr__(self): return "" % list(self._variable_dict.keys()) def __str__(self): stream = StringIO() # Attributes. has_attributes = False try: source_product = self.source_product except AttributeError: pass else: if source_product: print("source product = %r" % source_product, file=stream) has_attributes = True try: history = self.history except AttributeError: pass else: if history: print("history = %r" % history, file=stream) has_attributes = True # Variables. if not self._variable_dict: return stream.getvalue() if has_attributes: stream.write("\n") for name, variable in _dict_iteritems(self._variable_dict): if not isinstance(variable, Variable): print("" % name, file=stream) continue if not isinstance(variable.data, numpy.ndarray) and not numpy.isscalar(variable.data): print("" % name, file=stream) continue if isinstance(variable.data, numpy.ndarray) and variable.data.size == 0: print("" % name, file=stream) continue # Data type and variable name. stream.write(variable._format_data_type() + " " + name) # Dimensions. if variable.dimension: stream.write(" " + _format_dimensions(variable.dimension, variable.data)) # Unit. try: unit = variable.unit except AttributeError: pass else: if unit is not None: stream.write(" [%s]" % unit) stream.write("\n") return stream.getvalue() def to_dict(self): """Convert product to a dictionary (OrderedDict).""" return to_dict(self) def to_xarray(self): """Convert product to an xarray dataset.""" import xarray # TODO no _lib.harp_type_int64? # TODO coord var heuristics, parameter overrides # TODO pass through more attrs? data = {} vardims = {} for varname, var in self._variable_dict.items(): # scalar/array data if len(var.dimension) == 0: conv_dims = () else: conv_dims = [] for dim, dimsize in zip(var.dimension, var.data.shape): conv_dim = dim or 'independent_%d' % dimsize conv_dims.append(conv_dim) vardims[conv_dim] = dimsize # attrs attrs = None unit = getattr(var, 'unit', None) if unit is not None: attrs = { 'unit': unit, } # xarray datavar data[varname] = (conv_dims, var.data, attrs) coord_vars = { 'time': 'datetime', 'longitude': 'longitude', 'latitude': 'latitude', } # promote coordinate vars dims = {} for dim, varname in coord_vars.items(): if varname in data: convdata = self[varname].data # if dim == 'time': # convdata = [datetime.datetime(2000, 1, 1) + datetime.timedelta(seconds=s) for s in convdata] dims[dim] = convdata del data[varname] return xarray.Dataset(data, dims) @staticmethod def from_xarray(dataset): """Convert dataset to HARP product.""" product = Product() for varname, var in (list(dataset.data_vars.items()) + list(dataset.coords.items())): dims = [None if d.startswith('independent_') else d for d in var.dims] variable = Variable(var.values, dims, var.attrs.get('unit')) setattr(product, varname, variable) return product def _get_c_library_filename(): """Return the filename of the HARP shared library depending on the current platform. """ from platform import system as _system if _system() == "Windows": return "harp.dll" if _system() == "Darwin": library_name = "libharp.dylib" else: library_name = "libharp.so" import os.path # check for library file in the parent directory (for pyinstaller bundles) library_path = os.path.normpath(os.path.join(os.path.dirname(__file__), "..", library_name)) if os.path.exists(library_path): return library_path # Debian/Ubuntu st,reply = subprocess.getstatusoutput('dpkg-architecture --query DEB_HOST_MULTIARCH') if st != 0: raise Exception("Multiarch check failed for libharp: %s" % reply) library_path = os.path.join(os.path.normpath("/usr/lib"), reply, "libharp.so.10") if os.path.exists(library_path): return library_path # assume the library to be in the parent library directory library_path = os.path.normpath(os.path.join(os.path.dirname(__file__), "../../..", library_name)) if not os.path.exists(library_path): # on RHEL the python path uses lib64, but the library might have gotten installed into lib alt_library_path = os.path.normpath(os.path.join(os.path.dirname(__file__), "../../../../lib", library_name)) if os.path.exists(alt_library_path): return alt_library_path # on Ubuntu the python prefix can be 'local/lib', but libraries get installed in 'lib' alt_library_path = os.path.normpath(os.path.join(os.path.dirname(__file__), "../../../../../lib", library_name)) if os.path.exists(alt_library_path): return alt_library_path return library_path def _get_filesystem_encoding(): """Return the encoding used by the filesystem.""" from sys import getdefaultencoding as _getdefaultencoding, getfilesystemencoding as _getfilesystemencoding encoding = _getfilesystemencoding() if encoding is None: encoding = _getdefaultencoding() return encoding def _init(): """Initialize the HARP Python interface.""" global _lib, _encoding, _py_dimension_type, _c_dimension_type, _py_data_type, _c_data_type_name from platform import system as _system # Initialize the HARP C library clib = _get_c_library_filename() _lib = _ffi.dlopen(clib) if os.getenv('CODA_DEFINITION') is None: # Set coda definition path relative to C library relpath = "../share/coda/definitions" if _system() == "Windows": _lib.harp_set_coda_definition_path_conditional(_encode_path(os.path.basename(clib)), _ffi.NULL, _encode_path(relpath)) else: _lib.harp_set_coda_definition_path_conditional(_encode_path(os.path.basename(clib)), _encode_path(os.path.dirname(clib)), _encode_path(relpath)) if os.getenv('UDUNITS2_XML_PATH') is None: # Set udunits2 xml path relpath = "../share/xml/udunits/udunits2.xml" if _system() == "Windows": _lib.harp_set_udunits2_xml_path_conditional(_encode_path(os.path.basename(clib)), _ffi.NULL, _encode_path(relpath)) else: _lib.harp_set_udunits2_xml_path_conditional(_encode_path(os.path.basename(clib)), _encode_path(os.path.dirname(clib)), _encode_path(relpath)) if _lib.harp_init() != 0: raise CLibraryError() # Set default encoding. _encoding = "ascii" # Initialize various look-up tables used throughout the HARP Python interface (i.e. this module). _py_dimension_type = \ { _lib.harp_dimension_independent: None, _lib.harp_dimension_time: "time", _lib.harp_dimension_latitude: "latitude", _lib.harp_dimension_longitude: "longitude", _lib.harp_dimension_vertical: "vertical", _lib.harp_dimension_spectral: "spectral" } _c_dimension_type = \ { None: _lib.harp_dimension_independent, "time": _lib.harp_dimension_time, "latitude": _lib.harp_dimension_latitude, "longitude": _lib.harp_dimension_longitude, "vertical": _lib.harp_dimension_vertical, "spectral": _lib.harp_dimension_spectral } _py_data_type = \ { _lib.harp_type_int8: numpy.int8, _lib.harp_type_int16: numpy.int16, _lib.harp_type_int32: numpy.int32, _lib.harp_type_float: numpy.float32, _lib.harp_type_double: numpy.float64, _lib.harp_type_string: numpy.object_ } _c_data_type_name = \ { _lib.harp_type_int8: "byte", _lib.harp_type_int16: "int", _lib.harp_type_int32: "long", _lib.harp_type_float: "float", _lib.harp_type_double: "double", _lib.harp_type_string: "string" } def _dict_iteritems(dictionary): """Get an iterator or view on the items of the specified dictionary. This method is Python 2 and Python 3 compatible. """ try: return iter(list(dictionary.items())) except AttributeError: return list(dictionary.items()) def _get_py_dimension_type(dimension_type): """Return the dimension name corresponding to the specified C dimension type code. """ try: return _py_dimension_type[dimension_type] except KeyError: raise UnsupportedDimensionError("unsupported C dimension type code '%d'" % dimension_type) def _get_c_dimension_type(dimension_type): """Return the C dimension type code corresponding to the specified dimension name. """ try: return _c_dimension_type[dimension_type] except KeyError: raise UnsupportedDimensionError("unsupported dimension %r" % dimension_type) def _get_py_data_type(data_type): """Return the Python type corresponding to the specified C data type code.""" try: return _py_data_type[data_type] except KeyError: raise UnsupportedTypeError("unsupported C data type code '%d'" % data_type) def _get_c_data_type(value): """Return the C data type code corresponding to the specified variable data value. """ if isinstance(value, (numpy.ndarray, numpy.generic)): # For NumPy ndarrays and scalars, determine the smallest HARP C data type that can safely contain elements of # the ndarray or scalar dtype. if numpy.issubdtype(value.dtype, numpy.object_): # NumPy object arrays are only used to contain variable length strings or byte strings. if all(isinstance(x, str) for x in value.flat): return _lib.harp_type_string elif all(isinstance(x, bytes) for x in value.flat): return _lib.harp_type_string else: raise UnsupportedTypeError("elements of a NumPy object array must be all str or all bytes") elif numpy.issubdtype(value.dtype, numpy.str_) or numpy.issubdtype(value.dtype, numpy.bytes_): return _lib.harp_type_string elif numpy.can_cast(value.dtype, numpy.int8): return _lib.harp_type_int8 elif numpy.can_cast(value.dtype, numpy.int16): return _lib.harp_type_int16 elif numpy.can_cast(value.dtype, numpy.int32): return _lib.harp_type_int32 elif numpy.can_cast(value.dtype, numpy.float32): return _lib.harp_type_float elif numpy.can_cast(value.dtype, numpy.float64): return _lib.harp_type_double else: raise UnsupportedTypeError("unsupported NumPy dtype '%s'" % value.dtype) elif numpy.isscalar(value): if isinstance(value, (str, bytes)): return _lib.harp_type_string try: numpy.array(value, dtype=numpy.int8) return _lib.harp_type_int8 except Exception: pass try: numpy.array(value, dtype=numpy.int16) return _lib.harp_type_int16 except Exception: pass try: numpy.array(value, dtype=numpy.int32) return _lib.harp_type_int32 except Exception: pass try: numpy.array(value, dtype=numpy.float32) return _lib.harp_type_float except Exception: pass try: numpy.array(value, dtype=numpy.float64) return _lib.harp_type_double except Exception: pass raise UnsupportedTypeError("unsupported type %r" % value.__class__.__name__) else: raise UnsupportedTypeError("unsupported type %r" % value.__class__.__name__) def _get_c_data_type_name(data_type): """Return the canonical name for the specified C data type code.""" try: return _c_data_type_name[data_type] except KeyError: raise UnsupportedTypeError("unsupported C data type code '%d'" % data_type) def _c_can_cast(c_data_type_src, c_data_type_dst): """Returns True if the source C data type can be cast to the destination C data type while preserving values. """ if c_data_type_dst == _lib.harp_type_int8: return (c_data_type_src == _lib.harp_type_int8) elif c_data_type_dst == _lib.harp_type_int16: return (c_data_type_src == _lib.harp_type_int8 or c_data_type_src == _lib.harp_type_int16) elif c_data_type_dst == _lib.harp_type_int32: return (c_data_type_src == _lib.harp_type_int8 or c_data_type_src == _lib.harp_type_int16 or c_data_type_src == _lib.harp_type_int32) elif c_data_type_dst == _lib.harp_type_float: return (c_data_type_src == _lib.harp_type_int8 or c_data_type_src == _lib.harp_type_int16 or c_data_type_src == _lib.harp_type_float or c_data_type_src == _lib.harp_type_double) elif c_data_type_dst == _lib.harp_type_double: return (c_data_type_src == _lib.harp_type_int8 or c_data_type_src == _lib.harp_type_int16 or c_data_type_src == _lib.harp_type_int32 or c_data_type_src == _lib.harp_type_float or c_data_type_src == _lib.harp_type_double) elif c_data_type_dst == _lib.harp_type_string: return (c_data_type_src == _lib.harp_type_string) else: return False def _encode_string_with_encoding(string, encoding="utf-8"): """Encode a unicode string using the specified encoding. By default, use the "surrogateescape" error handler to deal with encoding errors. This error handler ensures that invalid bytes encountered during decoding are converted to the same bytes during encoding, by decoding them to a special range of unicode code points. The "surrogateescape" error handler is available since Python 3.1. For earlier versions of Python 3, the "strict" error handler is used instead. """ try: try: return string.encode(encoding, "surrogateescape") except LookupError: # Either the encoding or the error handler is not supported; fall-through to the next try-block. pass try: return string.encode(encoding) except LookupError: # Here it is certain that the encoding is not supported. raise Error("unknown encoding '%s'" % encoding) except UnicodeEncodeError: raise Error("cannot encode '%s' using encoding '%s'" % (string, encoding)) def _decode_string_with_encoding(string, encoding="utf-8"): """Decode a byte string using the specified encoding. By default, use the "surrogateescape" error handler to deal with encoding errors. This error handler ensures that invalid bytes encountered during decoding are converted to the same bytes during encoding, by decoding them to a special range of unicode code points. The "surrogateescape" error handler is available since Python 3.1. For earlier versions of Python 3, the "strict" error handler is used instead. This may cause decoding errors if the input byte string contains bytes that cannot be decoded using the specified encoding. Since most HARP products use ASCII strings exclusively, it is unlikely this will occur often in practice. """ try: try: return string.decode(encoding, "surrogateescape") except LookupError: # Either the encoding or the error handler is not supported; fall-through to the next try-block. pass try: return string.decode(encoding) except LookupError: # Here it is certain that the encoding is not supported. raise Error("unknown encoding '%s'" % encoding) except UnicodeEncodeError: raise Error("cannot decode '%s' using encoding '%s'" % (string, encoding)) def _encode_path(path): """Encode the input unicode path using the filesystem encoding. On Python 2, this method returns the specified path unmodified. """ if isinstance(path, bytes): # This branch will be taken for instances of class str on Python 2 (since this is an alias for class bytes), and # on Python 3 for instances of class bytes. return path elif isinstance(path, str): # This branch will only be taken for instances of class str on Python 3. On Python 2 such instances will take # the branch above. return _encode_string_with_encoding(path, _get_filesystem_encoding()) else: raise TypeError("path must be bytes or str, not %r" % path.__class__.__name__) def _encode_string(string): """Encode the input unicode string using the package default encoding. On Python 2, this method returns the specified string unmodified. """ if isinstance(string, bytes): # This branch will be taken for instances of class str on Python 2 (since this is an alias for class bytes), and # on Python 3 for instances of class bytes. return string elif isinstance(string, str): # This branch will only be taken for instances of class str on Python 3. On Python 2 such instances will take # the branch above. return _encode_string_with_encoding(string, get_encoding()) else: raise TypeError("string must be bytes or str, not %r" % string.__class__.__name__) def _decode_string(string): """Decode the input byte string using the package default encoding. On Python 2, this method returns the specified byte string unmodified. """ if isinstance(string, str): # This branch will be taken for instances of class str on Python 2 and Python 3. return string elif isinstance(string, bytes): # This branch will only be taken for instances of class bytes on Python 3. On Python 2 such instances will take # the branch above. return _decode_string_with_encoding(string, get_encoding()) else: raise TypeError("string must be bytes or str, not %r" % string.__class__.__name__) def _format_dimensions(dimension, data): """Construct a formatted string from the specified dimensions and data that provides information about dimension types and lengths, or "" if this information cannot be determined. """ if not isinstance(data, numpy.ndarray) or data.ndim != len(dimension): return "{}" stream = StringIO() stream.write("{") for i in range(data.ndim): if dimension[i]: stream.write(dimension[i] + "=") stream.write(str(data.shape[i])) if (i + 1) < data.ndim: stream.write(", ") stream.write("}") return stream.getvalue() def _import_scalar(c_data_type, c_data): if c_data_type == _lib.harp_type_int8: return c_data.int8_data elif c_data_type == _lib.harp_type_int16: return c_data.int16_data elif c_data_type == _lib.harp_type_int32: return c_data.int32_data elif c_data_type == _lib.harp_type_float: return c_data.float_data elif c_data_type == _lib.harp_type_double: return c_data.double_data raise UnsupportedTypeError("unsupported C data type code '%d'" % c_data_type) def _import_array(c_data_type, c_num_elements, c_data): if c_data_type == _lib.harp_type_string: data = numpy.empty((c_num_elements,), dtype=numpy.object_) for i in range(c_num_elements): # NB. The _ffi.string() method returns a copy of the C string. data[i] = _decode_string(_ffi.string(c_data.string_data[i])) return data # NB. The _ffi.buffer() method, as well as the numpy.frombuffer() method, provide a view on the C array; neither # method incurs a copy. c_data_buffer = _ffi.buffer(c_data.ptr, c_num_elements * _lib.harp_get_size_for_type(c_data_type)) return numpy.copy(numpy.frombuffer(c_data_buffer, dtype=_get_py_data_type(c_data_type))) def _import_variable(c_variable): # Import variable data. data = _import_array(c_variable.data_type, c_variable.num_elements, c_variable.data) num_dimensions = c_variable.num_dimensions if num_dimensions == 0: variable = Variable(data.item()) else: data = data.reshape([c_variable.dimension[i] for i in range(num_dimensions)]) dimension = [_get_py_dimension_type(c_variable.dimension_type[i]) for i in range(num_dimensions)] variable = Variable(data, dimension) # Import variable attributes. if c_variable.unit != _ffi.NULL: variable.unit = _decode_string(_ffi.string(c_variable.unit)) if c_variable.data_type != _lib.harp_type_string: variable.valid_min = _import_scalar(c_variable.data_type, c_variable.valid_min) variable.valid_max = _import_scalar(c_variable.data_type, c_variable.valid_max) if c_variable.description: variable.description = _decode_string(_ffi.string(c_variable.description)) num_enum_values = c_variable.num_enum_values if num_enum_values > 0 and c_variable.enum_name != _ffi.NULL: variable.enum = [_decode_string(_ffi.string(c_variable.enum_name[i])) for i in range(num_enum_values)] return variable def _import_product(c_product): product = Product() # Import product attributes. if c_product.source_product: product.source_product = _decode_string(_ffi.string(c_product.source_product)) if c_product.history: product.history = _decode_string(_ffi.string(c_product.history)) # Import variables. for i in range(c_product.num_variables): c_variable_ptr = c_product.variable[i] variable = _import_variable(c_variable_ptr[0]) setattr(product, _decode_string(_ffi.string(c_variable_ptr[0].name)), variable) return product def _import_product_metadata(c_metadata): metadata = OrderedDict() metadata['filename'] = _decode_string(_ffi.string(c_metadata.filename)) if numpy.isinf(c_metadata.datetime_start): metadata['datetime_start'] = datetime.datetime.min else: metadata['datetime_start'] = datetime.datetime(2000, 1, 1) + datetime.timedelta(days=c_metadata.datetime_start) if numpy.isinf(c_metadata.datetime_stop): metadata['datetime_stop'] = datetime.datetime.max else: metadata['datetime_stop'] = datetime.datetime(2000, 1, 1) + datetime.timedelta(days=c_metadata.datetime_stop) metadata['time'] = c_metadata.dimension[0] metadata['latitude'] = c_metadata.dimension[1] metadata['longitude'] = c_metadata.dimension[2] metadata['vertical'] = c_metadata.dimension[3] metadata['spectral'] = c_metadata.dimension[4] metadata['format'] = _decode_string(_ffi.string(c_metadata.format)) metadata['source_product'] = _decode_string(_ffi.string(c_metadata.source_product)) if c_metadata.history: metadata['history'] = _decode_string(_ffi.string(c_metadata.history)) return metadata def _export_scalar(data, c_data_type, c_data): if c_data_type == _lib.harp_type_int8: c_data.int8_data = data elif c_data_type == _lib.harp_type_int16: c_data.int16_data = data elif c_data_type == _lib.harp_type_int32: c_data.int32_data = data elif c_data_type == _lib.harp_type_float: c_data.float_data = data elif c_data_type == _lib.harp_type_double: c_data.double_data = data else: raise UnsupportedTypeError("unsupported C data type code '%d'" % c_data_type) def _export_array(data, c_variable): if c_variable.data_type != _lib.harp_type_string: # NB. The _ffi.buffer() method as well as the numpy.frombuffer() method provide a view on the C array; neither # method incurs a copy. The numpy.copyto() method also works if the source array is a scalar, i.e. not an # instance of numpy.ndarray. size = c_variable.num_elements * _lib.harp_get_size_for_type(c_variable.data_type) shape = data.shape if isinstance(data, numpy.ndarray) else () c_data_buffer = _ffi.buffer(c_variable.data.ptr, size) c_data = numpy.reshape(numpy.frombuffer(c_data_buffer, dtype=_get_py_data_type(c_variable.data_type)), shape) numpy.copyto(c_data, data, casting="safe") elif isinstance(data, numpy.ndarray): for index, element in enumerate(data.flat): if _lib.harp_variable_set_string_data_element(c_variable, index, _encode_string(element)) != 0: raise CLibraryError() else: assert(c_variable.num_elements == 1) if _lib.harp_variable_set_string_data_element(c_variable, 0, _encode_string(data)) != 0: raise CLibraryError() def _export_variable(name, variable, c_product): data = getattr(variable, "data", None) if data is None: raise Error("no data or data is None") # Check dimensions dimension = getattr(variable, "dimension", []) if not dimension and isinstance(data, numpy.ndarray) and data.size != 1: raise Error("dimensions missing or incomplete") if dimension and (not isinstance(data, numpy.ndarray) or data.ndim != len(dimension)): raise Error("dimensions incorrect") # Determine C data type. c_data_type = variable._get_c_data_type() # Encode variable name. c_name = _encode_string(name) # Determine C dimension types and lengths. c_num_dimensions = len(dimension) c_dimension_type = [_get_c_dimension_type(dimension_name) for dimension_name in dimension] c_dimension = _ffi.NULL if not dimension else data.shape # Create C variable of the proper size. c_variable_ptr = _ffi.new("harp_variable **") if _lib.harp_variable_new(c_name, c_data_type, c_num_dimensions, c_dimension_type, c_dimension, c_variable_ptr) != 0: raise CLibraryError() # Add C variable to C product. if _lib.harp_product_add_variable(c_product, c_variable_ptr[0]) != 0: _lib.harp_variable_delete(c_variable_ptr[0]) raise CLibraryError() # The C variable has been successfully added to the C product. Therefore, the memory management of the C variable # is tied to the life time of the C product. If an error occurs, the memory occupied by the C variable will be # freed along with the C product. c_variable = c_variable_ptr[0] # Copy data into the C variable. _export_array(data, c_variable) # Variable attributes. if c_data_type != _lib.harp_type_string: try: valid_min = variable.valid_min except AttributeError: pass else: if isinstance(valid_min, numpy.ndarray): if valid_min.size == 1: valid_min = valid_min.flat[0] else: raise Error("valid_min attribute should be scalar") c_data_type_valid_min = _get_c_data_type(valid_min) if _c_can_cast(c_data_type_valid_min, c_data_type): _export_scalar(valid_min, c_data_type, c_variable.valid_min) else: raise Error("type '%s' of valid_min attribute incompatible with type '%s' of data" % (_get_c_data_type_name(c_data_type_valid_min), _get_c_data_type_name(c_data_type))) try: valid_max = variable.valid_max except AttributeError: pass else: if isinstance(valid_max, numpy.ndarray): if valid_max.size == 1: valid_max = valid_max.flat[0] else: raise Error("valid_max attribute should be scalar") c_data_type_valid_max = _get_c_data_type(valid_max) if _c_can_cast(c_data_type_valid_max, c_data_type): _export_scalar(valid_max, c_data_type, c_variable.valid_max) else: raise Error("type '%s' of valid_max attribute incompatible with type '%s' of data" % (_get_c_data_type_name(c_data_type_valid_max), _get_c_data_type_name(c_data_type))) try: unit = variable.unit except AttributeError: pass else: if unit is not None and _lib.harp_variable_set_unit(c_variable, _encode_string(unit)) != 0: raise CLibraryError() try: description = variable.description except AttributeError: pass else: if description and _lib.harp_variable_set_description(c_variable, _encode_string(description)) != 0: raise CLibraryError() try: enum = variable.enum except AttributeError: pass else: if enum and _lib.harp_variable_set_enumeration_values(c_variable, len(enum), [_ffi.new("char[]", _encode_string(name)) for name in enum]) != 0: raise CLibraryError() def _export_product(product, c_product): # Export product attributes. try: source_product = product.source_product except AttributeError: pass else: if source_product and _lib.harp_product_set_source_product(c_product, _encode_string(source_product)) != 0: raise CLibraryError() try: history = product.history except AttributeError: pass else: if history and _lib.harp_product_set_history(c_product, _encode_string(history)) != 0: raise CLibraryError() # Export variables. for name in product: try: _export_variable(name, product[name], c_product) except Error as _error: raise Error("variable '%r' could not be exported (%s)" % (name, str(_error))) def _update_history(product, command): line = datetime.datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%SZ") line += " [harp-%s] " % (version()) line += command try: product.history += "\n" + line except AttributeError: product.history = line def get_encoding(): """Return the encoding used to convert between unicode strings and C strings (only relevant when using Python 3). """ return _encoding def set_encoding(encoding): """Set the encoding used to convert between unicode strings and C strings (only relevant when using Python 3). """ global _encoding _encoding = encoding def version(): """Return the version of the HARP C library.""" return _decode_string(_ffi.string(_lib.libharp_version)) def to_dict(product): """Convert a Product to an OrderedDict. The OrderedDict representation provides direct access to the data associated with each variable. All product attributes and all variable attributes except the unit attribute are discarded as part of the conversion. The unit attribute of a variable is represented by adding a scalar variable of type string with the name of the corresponding variable suffixed with '_unit' as name and the unit as value. Arguments: product -- Product to convert. """ if not isinstance(product, Product): raise TypeError("product must be Product, not %r" % product.__class__.__name__) dictionary = OrderedDict() for name in product: variable = product[name] try: dictionary[name] = variable.data if variable.unit is not None: dictionary[name + "_unit"] = variable.unit except AttributeError: pass return dictionary def import_product(filename, operations="", options="", reduce_operations="", post_operations=""): """Import a product from a file. This will first try to import the file as an HDF4, HDF5, or netCDF file that complies to the HARP Data Format. If the file is not stored using the HARP format then it will try to import it using one of the available ingestion modules. If the filename argument is a list of filenames, a globbing (glob.glob()) pattern, or a list of globbing patterns then the harp.import_product() function will be called on each individual matching file. All imported products will then be appended into a single merged product and that merged product will be returned. If the filename argument is a .pth file, then the products referenced in the .pth file will be treated as a HARP Dataset and its merged content will be returned. Note that the `source_product` attribute of products in a HARP Dataset needs to be unique; if a dataset contains multiple products with the same `source_product` value then only the last product in the list will be kept. Arguments: filename -- Filename, file pattern, .pth file, or list of filenames/patterns/.pths of the product(s) to import operations -- Actions to apply as part of the import; should be specified as a semi-colon separated string of operations; in case a list of products is ingested these operations will be performed on each product individually before the data is merged. options -- Ingestion module specific options; should be specified as a semi- colon separated string of key=value pairs; only used if a file is not in HARP format. reduce_operations -- Actions to apply after each append; should be specified as a semi-colon separated string of operations; these operations will only be applied if the filename argument is a file pattern or a list of filenames/patterns; this advanced option allows for memory efficient application of time reduction operations (such as bin()) that would normally be provided as part of post_operations. post_operations -- Actions to apply after the list of products is merged; should be specified as a semi-colon separated string of operations; these operations will only be applied if the filename argument is a file pattern or a list of filenames/patterns. """ filenames = None if not (isinstance(filename, bytes) or isinstance(filename, str)): # Assume this is a list of filenames or patterns filenames = [] for file in filename: if '*' in file or '?' in file: # This is a globbing pattern filenames.extend(sorted(glob.glob(file))) else: filenames.append(file) elif '*' in filename or '?' in filename: # This is a globbing pattern filenames = sorted(glob.glob(filename)) elif filename.endswith('.pth'): filenames = [filename] if filenames is not None: expanded_filenames = [] for filename in filenames: if filename.endswith('.pth'): with open(filename, newline='') as csvfile: pathreader = csv.reader(csvfile) next(pathreader) # skip header for row in pathreader: expanded_filenames.append(row[0]) else: expanded_filenames.append(filename) filenames = expanded_filenames if len(filenames) == 0: raise Error("no files matching '%s'" % (filename)) # Return the merged concatenation of all products merged_product_ptr = None try: for file in filenames: c_product_ptr = _ffi.new("harp_product **") # Import the product as a C product. if _lib.harp_import(_encode_path(file), _encode_string(operations), _encode_string(options), c_product_ptr) != 0: raise CLibraryError() if _lib.harp_product_is_empty(c_product_ptr[0]) == 1: _lib.harp_product_delete(c_product_ptr[0]) else: if merged_product_ptr is None: merged_product_ptr = c_product_ptr # if this remains the only product, make sure it still looks like it was the result of a merge if _lib.harp_product_append(merged_product_ptr[0], _ffi.NULL) != 0: raise CLibraryError() else: try: if _lib.harp_product_append(merged_product_ptr[0], c_product_ptr[0]) != 0: raise CLibraryError() finally: _lib.harp_product_delete(c_product_ptr[0]) if reduce_operations: # perform reduction operations on the partially merged product after each append if _lib.harp_product_execute_operations(merged_product_ptr[0], _encode_string(reduce_operations)) != 0: raise CLibraryError() except Exception: if merged_product_ptr is not None: _lib.harp_product_delete(merged_product_ptr[0]) raise if merged_product_ptr is None: raise NoDataError() try: if post_operations: if _lib.harp_product_execute_operations(merged_product_ptr[0], _encode_string(post_operations)) != 0: raise CLibraryError() # Convert the merged C product into its Python representation. product = _import_product(merged_product_ptr[0]) finally: _lib.harp_product_delete(merged_product_ptr[0]) if operations or options or post_operations: # Update history command = "harp.import_product('{0}'".format(filename) if operations: command += ",operations='{0}'".format(operations) if options: command += ",options='{0}'".format(options) if reduce_operations: command += ",reduce_operations='{0}'".format(reduce_operations) if post_operations: command += ",post_operations='{0}'".format(post_operations) command += ")" _update_history(product, command) return product c_product_ptr = _ffi.new("harp_product **") # Import the product as a C product. if _lib.harp_import(_encode_path(filename), _encode_string(operations), _encode_string(options), c_product_ptr) != 0: raise CLibraryError() try: # Raise an exception if the imported C product contains no variables, or variables without data. if _lib.harp_product_is_empty(c_product_ptr[0]) == 1: raise NoDataError() # Convert the C product into its Python representation. product = _import_product(c_product_ptr[0]) if operations or options: # Update history command = "harp.import_product('{0}'".format(filename) if operations: command += ",operations='{0}'".format(operations) if options: command += ",options='{0}'".format(options) command += ")" _update_history(product, command) return product finally: _lib.harp_product_delete(c_product_ptr[0]) def import_product_metadata(filename, options=""): """Import specific metadata from a single file. This will try to extract the following information from a file. - datetime_start - datetime_stop - dimension lengths for time, latitude, longitude, vertical, and spectral - source_product If the file is not stored using the HARP format then it will try to import the metadata using one of the available ingestion modules. Arguments: filename -- Filename of the product from which to extract the metadata options -- Ingestion module specific options; should be specified as a semi- colon separated string of key=value pairs; only used if a file is not in HARP format. """ c_metadata_ptr = _ffi.new("harp_product_metadata **") # Import the product as a C product. if _lib.harp_import_product_metadata(_encode_path(filename), _encode_string(options), c_metadata_ptr) != 0: raise CLibraryError() try: # Convert the C metadata into its Python representation. return _import_product_metadata(c_metadata_ptr[0]) finally: _lib.harp_product_metadata_delete(c_metadata_ptr[0]) def export_product(product, filename, file_format="netcdf", operations="", hdf5_compression=0): """Export a HARP compliant product. Arguments: product -- Product to export. filename -- Filename of the exported product. file_format -- File format to use; one of 'netcdf', 'hdf4', or 'hdf5'. operations -- Actions to apply as part of the export; should be specified as a semi-colon separated string of operations. hdf5_compression -- Compression level when exporting to hdf5 (0=disabled, 1=low, ..., 9=high). """ if not isinstance(product, Product): raise TypeError("product must be Product, not %r" % product.__class__.__name__) if operations: # Update history command = "harp.export_product('{0}', operations='{1}')".format(filename, operations) _update_history(product, command) # Create C product. c_product_ptr = _ffi.new("harp_product **") if _lib.harp_product_new(c_product_ptr) != 0: raise CLibraryError() try: # Convert the Python product to its C representation. _export_product(product, c_product_ptr[0]) if operations: # Apply operations to the product before export if _lib.harp_product_execute_operations(c_product_ptr[0], _encode_string(operations)) != 0: raise CLibraryError() # Don't allow export of empty products if _lib.harp_product_is_empty(c_product_ptr[0]) == 1: raise NoDataError() # Export the C product to a file. if file_format == 'hdf5': _lib.harp_set_option_hdf5_compression(int(hdf5_compression)) if _lib.harp_export(_encode_path(filename), _encode_string(file_format), c_product_ptr[0]) != 0: raise CLibraryError() finally: _lib.harp_product_delete(c_product_ptr[0]) def _get_time_length(product): time_length = None for name in product: variable = product[name] dimension = getattr(variable, "dimension", []) data = numpy.asarray(variable.data) if dimension and len(data.shape) != len(dimension): raise Error("dimensions incorrect") for i in range(len(dimension)): if dimension[i] == "time": if time_length is None: time_length = data.shape[i] elif time_length != data.shape[i]: raise Error("inconsistent dimension lengths for 'time'") return time_length def _extend_variable_for_dim(variable, dim_index, new_length): shape = list(variable.data.shape) shape[dim_index] = new_length - shape[dim_index] if variable.data.dtype.kind == 'f': filler = numpy.full(shape, numpy.nan, dtype=variable.data.dtype) else: filler = numpy.zeros(shape, dtype=variable.data.dtype) variable.data = numpy.concatenate([variable.data, filler], axis=dim_index) def make_time_dependent(product): time_length = _get_time_length(product) if time_length is None: raise Error("product has no time dimension") for name in product: variable = product[name] dimension = getattr(variable, "dimension", []) if len(dimension) == 0 or dimension[0] != 'time': # add time dimension data = numpy.asarray(variable.data) data = data.reshape([1] + list(data.shape)) variable.data = numpy.repeat(data, time_length, axis=0) variable.dimension = ['time'] + dimension return product def concatenate(products): if len(products) == 0: raise Error("product list is empty") variable_names = [] for product in products: for name in product: if name not in variable_names and name != "index": variable_names.append(name) for name in variable_names: for product in products: if name not in product: raise Error("not all products contain variable '%s'" % (name,)) for product in products: make_time_dependent(product) target_product = Product() for name in variable_names: for product in products: source_variable = product[name] if name not in target_product: target_variable = Variable(source_variable.data, source_variable.dimension) if hasattr(source_variable, 'unit'): target_variable.unit = source_variable.unit if hasattr(source_variable, 'valid_min'): target_variable.valid_min = source_variable.valid_min if hasattr(source_variable, 'valid_max'): target_variable.valid_max = source_variable.valid_max if hasattr(source_variable, 'description'): target_variable.description = source_variable.description if hasattr(source_variable, 'enum'): target_variable.enum = source_variable.enum target_product[name] = target_variable else: target_variable = target_product[name] if hasattr(target_variable, 'unit'): if not hasattr(source_variable, 'unit') or target_variable.unit != source_variable.unit: raise Error("inconsistent units in appending variable '%s'" % (name,)) if len(target_variable.data.shape) != len(source_variable.data.shape): raise Error("inconsistent number of dimensions for appending variable '%s'" % (name,)) for i in range(len(target_variable.data.shape))[1:]: if target_variable.data.shape[i] < source_variable.data.shape[i]: _extend_variable_for_dim(target_variable, i, source_variable.data.shape[i]) if source_variable.data.shape[i] < target_variable.data.shape[i]: _extend_variable_for_dim(source_variable, i, target_variable.data.shape[i]) target_variable.data = numpy.append(target_variable.data, source_variable.data, axis=0) return target_product def execute_operations(products, operations="", post_operations=""): if isinstance(products, Product): if not operations: return products # Create C product. c_product_ptr = _ffi.new("harp_product **") if _lib.harp_product_new(c_product_ptr) != 0: raise CLibraryError() try: _export_product(products, c_product_ptr[0]) if _lib.harp_product_execute_operations(c_product_ptr[0], _encode_string(operations)) != 0: raise CLibraryError() product = _import_product(c_product_ptr[0]) return product finally: _lib.harp_product_delete(c_product_ptr[0]) else: # Return the merged concatenation of all products merged_product_ptr = None try: for product in products: c_product_ptr = _ffi.new("harp_product **") if _lib.harp_product_new(c_product_ptr) != 0: raise CLibraryError() try: _export_product(product, c_product_ptr[0]) if _lib.harp_product_execute_operations(c_product_ptr[0], _encode_string(operations)) != 0: raise CLibraryError() except Exception: _lib.harp_product_delete(c_product_ptr[0]) raise if _lib.harp_product_is_empty(c_product_ptr[0]) == 1: _lib.harp_product_delete(c_product_ptr[0]) elif merged_product_ptr is None: merged_product_ptr = c_product_ptr # if this remains the only product then make sure it still looks like it was the result of a merge if _lib.harp_product_append(merged_product_ptr[0], _ffi.NULL) != 0: raise CLibraryError() else: try: if _lib.harp_product_append(merged_product_ptr[0], c_product_ptr[0]) != 0: raise CLibraryError() finally: _lib.harp_product_delete(c_product_ptr[0]) except Exception: if merged_product_ptr is not None: _lib.harp_product_delete(merged_product_ptr[0]) raise if merged_product_ptr is None: raise NoDataError() try: if post_operations: if _lib.harp_product_execute_operations(merged_product_ptr[0], _encode_string(post_operations)) != 0: raise CLibraryError() # Convert the merged C product into its Python representation. product = _import_product(merged_product_ptr[0]) finally: _lib.harp_product_delete(merged_product_ptr[0]) return product def convert_unit(from_unit, to_unit, values): values = numpy.array(values, dtype=numpy.double) c_data = _ffi.cast("double *", values.ctypes.data) if _lib.harp_convert_unit_double(_encode_string(from_unit), _encode_string(to_unit), numpy.size(values), c_data) != 0: raise CLibraryError() return values # # Initialize the HARP Python interface. # _init()