From 3ddedde7422465c7acf667bfee37e61017556e8e Mon Sep 17 00:00:00 2001 From: bekozi Date: Mon, 3 Aug 2020 08:11:40 -0500 Subject: [PATCH] ENH: Improved handling of self-intersecting polygons (#514) - Added ESMF mesh conversion/validation script - Minor documentation improvements --- doc/sphinx_examples/chunked_rwg_help.sh | 2 +- examples/convert_validate_esmf_mesh.py | 135 +++ src/ocgis/exc.py | 965 +++++++++--------- .../test_spatial/test_grid_chunker.py | 104 +- .../test_ocgis/test_variable/test_geom.py | 48 +- src/ocgis/util/logging_ocgis.py | 5 +- src/ocgis/variable/geom.py | 92 +- 7 files changed, 808 insertions(+), 543 deletions(-) create mode 100644 examples/convert_validate_esmf_mesh.py diff --git a/doc/sphinx_examples/chunked_rwg_help.sh b/doc/sphinx_examples/chunked_rwg_help.sh index 20301fed6..3ba18ad6a 100644 --- a/doc/sphinx_examples/chunked_rwg_help.sh +++ b/doc/sphinx_examples/chunked_rwg_help.sh @@ -1,6 +1,6 @@ $ ocli chunked-rwg --help -Usage: ocli.py chunked-rwg [OPTIONS] +Usage: ocli chunked-rwg [OPTIONS] Generate regridding weights using a spatial decomposition. diff --git a/examples/convert_validate_esmf_mesh.py b/examples/convert_validate_esmf_mesh.py new file mode 100644 index 000000000..d09c64227 --- /dev/null +++ b/examples/convert_validate_esmf_mesh.py @@ -0,0 +1,135 @@ +""" +Test conversion of GeoPackage GIS data to ESMF Unstructured validating each mesh element along the way. +""" +import os +import random +import shutil +import tempfile + +import ESMF +import numpy as np + +import ocgis + +ESMF.Manager(debug=True) + +# Path to the source GIS data assuming a GeoPackage +GEOPKG = 'hdma_global_catch_0.gpkg' +# This is the template for catchment-specific directories. If the test operation is successful, this will deleted. If it +# is not successful, the directory will be left behind. +OUTDIR_TEMPLATE = os.path.join(tempfile.gettempdir(), 'individual-element-nc', 'hruid-tmp-{}') +# Debugging to simulate errors +DEBUG = False +# Number of nodes for each virtual polygon +NODE_THRESHOLD = 5000 +# De-duplicate nodes (serial only). Leave this off for initial testing since there is only one node coordinates array +PACK = False +# Whether to split holes/interiors. Start with False just to use exteriors +SPLIT_INTERIORS = False + + +def do_destroy(to_destroy): + for t in to_destroy: + try: + t.destroy() + except: + pass + + +def do_esmf(ncpath, exedir, dst_logdir): + # Will return True if the operation is successful + was_successful = False + # Destroy the ESMF objects so deleting them if there is a failure is smoother + mesh, src, dst, regrid = [None]*4 + try: + # Fake debug error simulation + if DEBUG and random.random() > 0.5: + raise ValueError('a fake debug error for testing') + # Create the ESMF mesh + mesh = ESMF.Mesh(filename=ncpath, filetype=ESMF.constants.FileFormat.ESMFMESH) + # Create the source + src = ESMF.Field(mesh, ndbounds=np.array([1, 1]), meshloc=ESMF.constants.MeshLoc.ELEMENT) + # Create the destination + dst = ESMF.Field(mesh, ndbounds=np.array([1, 1]), meshloc=ESMF.constants.MeshLoc.ELEMENT) + # This will create the route handle and return some weights + regrid = ESMF.Regrid(srcfield=src, dstfield=dst, regrid_method=ESMF.constants.RegridMethod.CONSERVE, factors=True) + factors = regrid.get_weights_dict(deep_copy=True) + assert factors is not None + was_successful = True + except Exception as e: + print('ERROR: {}: ESMF read/regridding problem with file path {}'.format(str(e), ncpath)) + # Copy over the ESMF log file + srcpath = os.path.join(exedir, 'PET0.ESMF_LogFile') + dstpath = os.path.join(dst_logdir, 'PET0.ESMF_LogFile') + if os.path.exists(srcpath): + shutil.copy2(srcpath, dstpath) + raise + finally: + do_destroy([mesh, src, dst, regrid]) + return was_successful + + +def do_record_test(exedir, record): + # The record's unique HRU identifier + hruid = record['properties']['hruid'] + # The current output directory + curr_outdir = OUTDIR_TEMPLATE.format(hruid) + # Create the directory + os.makedirs(curr_outdir, exist_ok=True) + # Make that the current working directory + os.chdir(curr_outdir) + # We need to transform the coordinate system from WGS84 to Spherical for ESMF + crs = ocgis.crs.CoordinateReferenceSystem(value=record['meta']['crs']) + field = ocgis.Field.from_records([record], crs=crs) + field.update_crs(ocgis.crs.Spherical()) + # Convert the field geometry to an unstructured grid format based on the UGRID spec. + gc = field.geom.convert_to( + pack=PACK, + node_threshold=NODE_THRESHOLD, + split_interiors=SPLIT_INTERIORS, + remove_self_intersects=False, + allow_splitting_excs=False + ) + # Path to the output netCDF file for the current element + out_element_nc = os.path.join(curr_outdir, "esmf-element_hruid-{}.nc".format(hruid)) + # Add the center coordinate to make ESMF happy (even though we are not using it) + centerCoords = np.array([field.geom.v()[0].centroid.x, field.geom.v()[0].centroid.y]).reshape(1, 2) + ocgis.Variable(name='centerCoords', value=centerCoords, dimensions=['elementCount', 'coordDim'], attrs={'units': 'degrees'}, parent=gc.parent) + # When writing the data to file, convert to ESMF unstructured format. + gc.parent.write(out_element_nc, driver='netcdf-esmf-unstruct') + # Run the simple regridding test + success = do_esmf(out_element_nc, exedir, curr_outdir) + if success: + # If successful, remove the directory + assert 'hruid-tmp-' in curr_outdir + shutil.rmtree(curr_outdir) + else: + # If it's not successful, leave the directory. Write the shapefile so it's easy to look at. Also send the record + # string to a file. + field.geom.write_vector('02-problem-hruid-{}.shp'.format(hruid)) + record_out = '01-problem-record-hruid-{}.out'.format(hruid) + with open(record_out, 'w') as f: + f.write(str(record)) + # Change back to the execution directory + os.chdir(exedir) + # Try to remove the log file if it exists. + if os.path.exists('PET0.ESMF_LogFile'): + os.remove('PET0.ESMF_LogFile') + + +def main(): + # The execution directory + exedir = os.getcwd() + # An iterator over the geometries in the GeoPackage + gc = ocgis.GeomCabinetIterator(path=GEOPKG, driver_kwargs={'feature_class': 'Catchment'}) + # Number of records in the GeoPackage + len_gc = len(gc) + for ctr, record in enumerate(gc): + # Print an update every 100 iterations + if ctr % 100 == 0: + print('INFO: Index {} of {}'.format(ctr, len_gc)) + do_record_test(exedir, record) + + +if __name__ == "__main__": + main() diff --git a/src/ocgis/exc.py b/src/ocgis/exc.py index 2103af75e..3716199be 100644 --- a/src/ocgis/exc.py +++ b/src/ocgis/exc.py @@ -1,480 +1,485 @@ -from ocgis import messages - - -class OcgException(Exception): - """Base class for all OCGIS exceptions.""" - omessage = None - - def __init__(self, message=None): - if self.omessage is not None and message is None: - message = self.omessage - self.message = message - - def __str__(self): - return self.message - - -class OcgWarning(Warning): - """Base class for all OCGIS warnings.""" - - -######################################################################################################################## - - -class BoundsAlreadyAvailableError(OcgException): - """Raised when an attempt is made to extrapolate bounds and they are already present.""" - - def __str__(self): - msg = 'Bounds/corners already available.' - return msg - - -class CannotFormatTimeError(OcgException): - """ - Raised when datetime objects from numeric are blocked by "format_time". - """ - - def __init__(self, property_name): - self.property_name = property_name - - def __str__(self): - msg = 'Attempted to retrieve datetime values from "{0}" with "format_time" as "False". Set "format_time" to "True".'.format( - self.property_name) - return msg - - -class MaskedDataFound(OcgException): - """Raised when data is masked.""" - - def __str__(self): - msg = 'Data is masked.' - return msg - - -class MultipleElementsFound(OcgException): - """ - Raised when multiple elements are encountered in a :class:`ocgis.interface.base.dimension.spatial.SpatialDimension` - object. - - :param sdim: The incoming spatial dimension object. - :type sdim: :class:`ocgis.interface.base.dimension.spatial.SpatialDimension` - """ - - def __init__(self, sdim): - self.sdim = sdim - - def __str__(self): - msg = 'Shape of the spatial dimension object is: {0}'.format(self.sdim.shape) - return msg - - -class RequestableFeature(OcgException): - """Raised when a feature is not implemented but could be provided there is a use case available.""" - - def __init__(self, message=None): - from ocgis import constants - self.message = 'This feature is not implemented. Do you need it? Please post an issue: {}. {}' - if message is None: - second_format = '' - else: - second_format = "Custom message is: '{}'".format(message) - self.message = self.message.format(constants.GITHUB_ISSUES, second_format) - - -class ShapeError(OcgException): - """ - Raised when an array has an incompatible shape with an operation. - """ - - -class SingleElementError(ShapeError): - """ - Raised when an operation requires more than a single element. - """ - - -class CalculationException(OcgException): - def __init__(self, function_klass, message=None): - self.function_klass = function_klass - OcgException.__init__(self, message=message) - - def __str__(self): - msg = 'The function class "{0}" raised an exception with message: "{1}"'.format(self.function_klass.__name__, - self.message) - return (msg) - - -class VariableInCollectionError(OcgException): - def __init__(self, variable_or_name): - try: - name = variable_or_name.name - except AttributeError: - name = variable_or_name - self.name = name - - def __str__(self): - msg = 'Variable name already in collection: {0}'.format(self.name) - return msg - - -class VariableShapeMismatch(OcgException): - def __init__(self, variable, collection_shape): - self.variable = variable - self.collection_shape = collection_shape - - def __str__(self): - msg = 'Variable with alias "{0}" has shape {1}. Collection shape is {2}'.format(self.variable.alias, - self.variable.shape, - self.collection_shape) - return msg - - -class SampleSizeNotImplemented(CalculationException): - pass - - -class InterpreterException(OcgException): - pass - - -class InterpreterNotRecognized(InterpreterException): - pass - - -class EmptyIterationError(OcgException): - def __init__(self, obj): - self.message = 'Iteration on the object "{0}" requested, but the object is empty.'.format(obj) - - -class CFException(OcgException): - pass - - -class ProjectionCoordinateNotFound(CFException): - def __init__(self, target): - self.message = 'The projection coordinate "{0}" was not found in the dataset.'.format(target) - - -class ProjectionDoesNotMatch(CFException): - pass - - -class CRSDepthNotImplemented(CFException): - """Raised when a CRS depth check is not supported.""" - - def __init__(self, depth): - self.message = "Depth check not supported: {}".format(depth) - - -class DimensionNotFound(CFException): - def __init__(self, name): - self.message = "Dimension not found: '{}'".format(name) - - -class DefinitionValidationError(OcgException): - """Raised when validation fails on :class:`~ocgis.OcgOperations`. - - :param ocg_argument: The origin of the exception. - :type ocg_argument: :class:`ocgis.driver.definition.AbstractParameter`, str - :param msg: The message related to the exception to display in the exception's template. - :type msg: str - """ - - def __init__(self, ocg_argument, msg): - self.ocg_argument = ocg_argument - - fmt = ('OcgOperations validation raised an exception on the argument/operation ' - '"{0}" with the message: {1}') - try: - msg = fmt.format(ocg_argument.name, msg) - except AttributeError: - try: - msg = fmt.format(ocg_argument._name, msg) - except AttributeError: - msg = fmt.format(ocg_argument, msg) - - self.message = msg - - -class ParameterFormattingError(OcgException): - pass - - -class UniqueIdNotFound(OcgException): - def __init__(self): - self.message = 'No unique ids found.' - - -class DummyDimensionEncountered(OcgException): - pass - - -class ResolutionError(OcgException): - pass - - -class SubsetException(OcgException): - """Base class for all subset exceptions.""" - pass - - -class OcgisEnvironmentError(OcgException): - def __init__(self, env_parm, msg): - self.env_parm = env_parm - self.msg = msg - - def __str__(self): - new_msg = 'Error when setting the ocgis.env variable {0}. The message is: {1}'.format(self.env_parm.name, - self.msg) - return (new_msg) - - -class SpatialWrappingError(OcgException): - """Raised for errors related to wrapping or unwrapping a geographic coordinate system.""" - pass - - -class MaskedDataError(SubsetException): - def __init__(self): - self.message = 'Geometric intersection returned all masked values.' - - -class ExtentError(SubsetException): - def __init__(self, message=None): - self.message = message or 'Geometry intersection is empty.' - - -class TemporalExtentError(SubsetException): - def __init__(self): - self.message = 'Temporal subset returned empty.' - - -class EmptyDataNotAllowed(SubsetException): - """Raised when the empty set for a geometry is returned and ``allow_empty`` is ``False``.""" - - def __init__(self): - self.message = 'Intersection returned empty, but empty data not allowed.' - - -class EmptyData(SubsetException): - def __init__(self, message=None, origin=None): - self.message = message or 'Empty data returned.' - self.origin = origin - - -class EmptySubsetError(SubsetException): - def __init__(self, origin=None): - self.origin = origin - - def __str__(self): - msg = 'A subset operation on variable "{0}" returned empty.'.format(self.origin) - return msg - - -class AllElementsMaskedError(OcgException): - """Raised when all elements are masked.""" - - def __str__(self): - return "All elements are masked." - - -class PayloadProtectedError(OcgException): - def __init__(self, name): - self.name = name - super(PayloadProtectedError, self).__init__() - - def __str__(self): - return 'Payload with name "{}" may not be loaded from source until "protected" is False.'.format(self.name) - - -class NoUnitsError(OcgException): - """Raised when a :class:`cfunits.Units` object is constructed from ``None``.""" - - def __init__(self, variable, message=None): - self.variable = variable - super(NoUnitsError, self).__init__(message=message) - - def __str__(self): - if self.message is None: - msg = 'Variable "{0}" has not been assigned units in the source metadata. Set the "units" attribute to continue.' - msg = msg.format(self.variable) - else: - msg = self.message - return (msg) - - -class UnitsValidationError(OcgException): - """Raised when units validation fails.""" - - def __init__(self, variable, required_units, calculation_key): - self.variable = variable - self.required_units = required_units - self.calculation_key = calculation_key - - def __str__(self): - msg = ('There was an error in units validation for calculation with key "{3}". The units on variable "{0}" ' - '(units="{2}") do not match the required units "{1}". The units should be conformed or overloaded if ' - 'incorrectly attributed.').format(self.variable.alias, self.required_units, self.variable.units, - self.calculation_key) - return msg - - -class IncompleteSeasonError(OcgException): - def __init__(self, season, year=None, month=None): - self.season = season - self.year = year - self.month = month - - def __str__(self): - if self.year is not None: - msg = 'The season specification "{0}" is missing the year "{1}".'.format(self.season, self.year + 1) - if self.month is not None: - msg = 'The season specification "{0}" is missing the month "{1}".'.format(self.season, self.month) - return msg - - -class VariableNotFoundError(OcgException): - """Raised when a variable name is not found in the target dataset.""" - - def __init__(self, uri, variable): - self.uri = uri - self.variable = variable - - def __str__(self): - msg = 'The variable "{0}" was not found in the dataset with URI: {1}'.format(self.variable, self.uri) - return msg - - -class VariableNotInCollection(OcgException): - """Raised when a variable is not present in a collection.""" - - def __init__(self, variable_name): - self.message = "The variable '{}' was not found in the collection.".format(variable_name) - - -class RegriddingError(OcgException): - """Raised for exceptions related to ESMPy-enabled regridding.""" - pass - - -class CornersInconsistentError(RegriddingError): - """Raised when corners are not available on all sources and/or destination fields.""" - pass - - -class RequestValidationError(OcgException): - """Raised when validation fails on a parameter when creating a :class:`~ocgis.RequestDataset` object.""" - - def __init__(self, keyword, message): - self.keyword = keyword - self.message = message - - def __str__(self): - message = 'Validation failed on the parameter "{0}" with the message: {1}'.format(self.keyword, self.message) - return message - - -class NoDataVariablesFound(RequestValidationError): - """Raised when no data variables are found in the target dataset.""" - - def __init__(self): - super(NoDataVariablesFound, self).__init__('variable', messages.M1) - - -class NoInteriorsError(OcgException): - def __init__(self, message=None): - if message is None: - message = 'Polygon has no holes or interiors. Nothing to do.' - super(NoInteriorsError, self).__init__(message) - - -class GridDeficientError(OcgException): - """Raised when a grid is missing parameters necessary to create a geometry array.""" - - -class DimensionMismatchError(OcgException): - """Raised when a variable's dimensions do not match those in the existing collection.""" - - def __init__(self, dim_name, vc_name, message=None): - self.dim_name = dim_name - self.vc_name = vc_name - - super(DimensionMismatchError, self).__init__(message) - - def __str__(self): - msg = 'The dimension "{}" does not match the dimension in variable collection "{}".'.format(self.dim_name, - self.vc_name) - return msg - - -class DimensionsRequiredError(OcgException): - """Raised when a variable requires dimensions.""" - - def __init__(self, message=None): - if message is None: - message = "Variables with dimension count greater than 0 (ndim > 0) require dimensions. Initialize the " \ - "variable with dimensions or call 'create_dimensions' before size inquiries (i.e. ndim, shape)." - super(DimensionsRequiredError, self).__init__(message=message) - - -class OcgDistError(OcgException): - """Raised for MPI-related exceptions.""" - - -class EmptyObjectError(OcgException): - """Raised when an empty object is not allowed.""" - - -class SubcommNotFoundError(OcgDistError): - """Raised when a subcommunicator is not found.""" - - def __init__(self, name): - message = "Subcommunicator '{}' not found.".format(name) - super(SubcommNotFoundError, self).__init__(message=message) - - -class SubcommAlreadyCreatedError(OcgDistError): - """Raised when a subcommunicator name already exists.""" - - def __init__(self, name): - message = "Subcommunicator '{}' already created.".format(name) - super(SubcommAlreadyCreatedError, self).__init__(message=message) - - -class CRSNotEquivalenError(OcgException): - """Raised when coordinate systems are not equivalent (not compatible for transform).""" - - def __init__(self, lhs, rhs): - msg = '{} is not equivalent to {}.'.format(lhs, rhs) - super(CRSNotEquivalenError, self).__init__(message=msg) - - -class DimensionMapError(OcgException): - """Raised when there is an issue with a dimension map entry.""" - - def __init__(self, entry_key, message): - msg = "Error with entry key '{}': {}".format(entry_key, message) - super(DimensionMapError, self).__init__(message=msg) - - -class VariableMissingMetadataError(OcgException): - """Raised when variable metadata cannot be found.""" - - def __init__(self, variable_name): - msg = 'Variable is missing metadata: {}'.format(variable_name) - super(VariableMissingMetadataError, self).__init__(message=msg) - - -class WrappedStateEvalTargetMissing(OcgException): - """Raised when attempting to retrieve the wrapped state of a field and no evaluation target is available.""" - - def __init__(self): - super(WrappedStateEvalTargetMissing, self).__init__(message='Target has no spatial information to evaluate.') - - -class NoTouching(OcgException): - """ - Raised when a subset geometry only touches an element. Used to prevent duplicate cells when creating a spatial - decomposition. - """ - omessage = "No touching allowed by spatial operation." +from ocgis import messages + + +class OcgException(Exception): + """Base class for all OCGIS exceptions.""" + omessage = None + + def __init__(self, message=None): + if self.omessage is not None and message is None: + message = self.omessage + self.message = message + + def __str__(self): + return self.message + + +class OcgWarning(Warning): + """Base class for all OCGIS warnings.""" + + +######################################################################################################################## + + +class BoundsAlreadyAvailableError(OcgException): + """Raised when an attempt is made to extrapolate bounds and they are already present.""" + + def __str__(self): + msg = 'Bounds/corners already available.' + return msg + + +class CannotFormatTimeError(OcgException): + """ + Raised when datetime objects from numeric are blocked by "format_time". + """ + + def __init__(self, property_name): + self.property_name = property_name + + def __str__(self): + msg = 'Attempted to retrieve datetime values from "{0}" with "format_time" as "False". Set "format_time" to "True".'.format( + self.property_name) + return msg + + +class MaskedDataFound(OcgException): + """Raised when data is masked.""" + + def __str__(self): + msg = 'Data is masked.' + return msg + + +class MultipleElementsFound(OcgException): + """ + Raised when multiple elements are encountered in a :class:`ocgis.interface.base.dimension.spatial.SpatialDimension` + object. + + :param sdim: The incoming spatial dimension object. + :type sdim: :class:`ocgis.interface.base.dimension.spatial.SpatialDimension` + """ + + def __init__(self, sdim): + self.sdim = sdim + + def __str__(self): + msg = 'Shape of the spatial dimension object is: {0}'.format(self.sdim.shape) + return msg + + +class RequestableFeature(OcgException): + """Raised when a feature is not implemented but could be provided there is a use case available.""" + + def __init__(self, message=None): + from ocgis import constants + self.message = 'This feature is not implemented. Do you need it? Please post an issue: {}. {}' + if message is None: + second_format = '' + else: + second_format = "Custom message is: '{}'".format(message) + self.message = self.message.format(constants.GITHUB_ISSUES, second_format) + + +class ShapeError(OcgException): + """ + Raised when an array has an incompatible shape with an operation. + """ + + +class SingleElementError(ShapeError): + """ + Raised when an operation requires more than a single element. + """ + + +class CalculationException(OcgException): + def __init__(self, function_klass, message=None): + self.function_klass = function_klass + OcgException.__init__(self, message=message) + + def __str__(self): + msg = 'The function class "{0}" raised an exception with message: "{1}"'.format(self.function_klass.__name__, + self.message) + return (msg) + + +class VariableInCollectionError(OcgException): + def __init__(self, variable_or_name): + try: + name = variable_or_name.name + except AttributeError: + name = variable_or_name + self.name = name + + def __str__(self): + msg = 'Variable name already in collection: {0}'.format(self.name) + return msg + + +class VariableShapeMismatch(OcgException): + def __init__(self, variable, collection_shape): + self.variable = variable + self.collection_shape = collection_shape + + def __str__(self): + msg = 'Variable with alias "{0}" has shape {1}. Collection shape is {2}'.format(self.variable.alias, + self.variable.shape, + self.collection_shape) + return msg + + +class SampleSizeNotImplemented(CalculationException): + pass + + +class InterpreterException(OcgException): + pass + + +class InterpreterNotRecognized(InterpreterException): + pass + + +class EmptyIterationError(OcgException): + def __init__(self, obj): + self.message = 'Iteration on the object "{0}" requested, but the object is empty.'.format(obj) + + +class CFException(OcgException): + pass + + +class ProjectionCoordinateNotFound(CFException): + def __init__(self, target): + self.message = 'The projection coordinate "{0}" was not found in the dataset.'.format(target) + + +class ProjectionDoesNotMatch(CFException): + pass + + +class CRSDepthNotImplemented(CFException): + """Raised when a CRS depth check is not supported.""" + + def __init__(self, depth): + self.message = "Depth check not supported: {}".format(depth) + + +class DimensionNotFound(CFException): + def __init__(self, name): + self.message = "Dimension not found: '{}'".format(name) + + +class DefinitionValidationError(OcgException): + """Raised when validation fails on :class:`~ocgis.OcgOperations`. + + :param ocg_argument: The origin of the exception. + :type ocg_argument: :class:`ocgis.driver.definition.AbstractParameter`, str + :param msg: The message related to the exception to display in the exception's template. + :type msg: str + """ + + def __init__(self, ocg_argument, msg): + self.ocg_argument = ocg_argument + + fmt = ('OcgOperations validation raised an exception on the argument/operation ' + '"{0}" with the message: {1}') + try: + msg = fmt.format(ocg_argument.name, msg) + except AttributeError: + try: + msg = fmt.format(ocg_argument._name, msg) + except AttributeError: + msg = fmt.format(ocg_argument, msg) + + self.message = msg + + +class ParameterFormattingError(OcgException): + pass + + +class UniqueIdNotFound(OcgException): + def __init__(self): + self.message = 'No unique ids found.' + + +class DummyDimensionEncountered(OcgException): + pass + + +class ResolutionError(OcgException): + pass + + +class SubsetException(OcgException): + """Base class for all subset exceptions.""" + pass + + +class OcgisEnvironmentError(OcgException): + def __init__(self, env_parm, msg): + self.env_parm = env_parm + self.msg = msg + + def __str__(self): + new_msg = 'Error when setting the ocgis.env variable {0}. The message is: {1}'.format(self.env_parm.name, + self.msg) + return (new_msg) + + +class SpatialWrappingError(OcgException): + """Raised for errors related to wrapping or unwrapping a geographic coordinate system.""" + pass + + +class MaskedDataError(SubsetException): + def __init__(self): + self.message = 'Geometric intersection returned all masked values.' + + +class ExtentError(SubsetException): + def __init__(self, message=None): + self.message = message or 'Geometry intersection is empty.' + + +class TemporalExtentError(SubsetException): + def __init__(self): + self.message = 'Temporal subset returned empty.' + + +class EmptyDataNotAllowed(SubsetException): + """Raised when the empty set for a geometry is returned and ``allow_empty`` is ``False``.""" + + def __init__(self): + self.message = 'Intersection returned empty, but empty data not allowed.' + + +class EmptyData(SubsetException): + def __init__(self, message=None, origin=None): + self.message = message or 'Empty data returned.' + self.origin = origin + + +class EmptySubsetError(SubsetException): + def __init__(self, origin=None): + self.origin = origin + + def __str__(self): + msg = 'A subset operation on variable "{0}" returned empty.'.format(self.origin) + return msg + + +class AllElementsMaskedError(OcgException): + """Raised when all elements are masked.""" + + def __str__(self): + return "All elements are masked." + + +class PayloadProtectedError(OcgException): + def __init__(self, name): + self.name = name + super(PayloadProtectedError, self).__init__() + + def __str__(self): + return 'Payload with name "{}" may not be loaded from source until "protected" is False.'.format(self.name) + + +class NoUnitsError(OcgException): + """Raised when a :class:`cfunits.Units` object is constructed from ``None``.""" + + def __init__(self, variable, message=None): + self.variable = variable + super(NoUnitsError, self).__init__(message=message) + + def __str__(self): + if self.message is None: + msg = 'Variable "{0}" has not been assigned units in the source metadata. Set the "units" attribute to continue.' + msg = msg.format(self.variable) + else: + msg = self.message + return (msg) + + +class UnitsValidationError(OcgException): + """Raised when units validation fails.""" + + def __init__(self, variable, required_units, calculation_key): + self.variable = variable + self.required_units = required_units + self.calculation_key = calculation_key + + def __str__(self): + msg = ('There was an error in units validation for calculation with key "{3}". The units on variable "{0}" ' + '(units="{2}") do not match the required units "{1}". The units should be conformed or overloaded if ' + 'incorrectly attributed.').format(self.variable.alias, self.required_units, self.variable.units, + self.calculation_key) + return msg + + +class IncompleteSeasonError(OcgException): + def __init__(self, season, year=None, month=None): + self.season = season + self.year = year + self.month = month + + def __str__(self): + if self.year is not None: + msg = 'The season specification "{0}" is missing the year "{1}".'.format(self.season, self.year + 1) + if self.month is not None: + msg = 'The season specification "{0}" is missing the month "{1}".'.format(self.season, self.month) + return msg + + +class VariableNotFoundError(OcgException): + """Raised when a variable name is not found in the target dataset.""" + + def __init__(self, uri, variable): + self.uri = uri + self.variable = variable + + def __str__(self): + msg = 'The variable "{0}" was not found in the dataset with URI: {1}'.format(self.variable, self.uri) + return msg + + +class VariableNotInCollection(OcgException): + """Raised when a variable is not present in a collection.""" + + def __init__(self, variable_name): + self.message = "The variable '{}' was not found in the collection.".format(variable_name) + + +class RegriddingError(OcgException): + """Raised for exceptions related to ESMPy-enabled regridding.""" + pass + + +class CornersInconsistentError(RegriddingError): + """Raised when corners are not available on all sources and/or destination fields.""" + pass + + +class RequestValidationError(OcgException): + """Raised when validation fails on a parameter when creating a :class:`~ocgis.RequestDataset` object.""" + + def __init__(self, keyword, message): + self.keyword = keyword + self.message = message + + def __str__(self): + message = 'Validation failed on the parameter "{0}" with the message: {1}'.format(self.keyword, self.message) + return message + + +class NoDataVariablesFound(RequestValidationError): + """Raised when no data variables are found in the target dataset.""" + + def __init__(self): + super(NoDataVariablesFound, self).__init__('variable', messages.M1) + + +class NoInteriorsError(OcgException): + def __init__(self, message=None): + if message is None: + message = 'Polygon has no holes or interiors. Nothing to do.' + super(NoInteriorsError, self).__init__(message) + + +class GridDeficientError(OcgException): + """Raised when a grid is missing parameters necessary to create a geometry array.""" + + +class DimensionMismatchError(OcgException): + """Raised when a variable's dimensions do not match those in the existing collection.""" + + def __init__(self, dim_name, vc_name, message=None): + self.dim_name = dim_name + self.vc_name = vc_name + + super(DimensionMismatchError, self).__init__(message) + + def __str__(self): + msg = 'The dimension "{}" does not match the dimension in variable collection "{}".'.format(self.dim_name, + self.vc_name) + return msg + + +class DimensionsRequiredError(OcgException): + """Raised when a variable requires dimensions.""" + + def __init__(self, message=None): + if message is None: + message = "Variables with dimension count greater than 0 (ndim > 0) require dimensions. Initialize the " \ + "variable with dimensions or call 'create_dimensions' before size inquiries (i.e. ndim, shape)." + super(DimensionsRequiredError, self).__init__(message=message) + + +class OcgDistError(OcgException): + """Raised for MPI-related exceptions.""" + + +class EmptyObjectError(OcgException): + """Raised when an empty object is not allowed.""" + + +class SubcommNotFoundError(OcgDistError): + """Raised when a subcommunicator is not found.""" + + def __init__(self, name): + message = "Subcommunicator '{}' not found.".format(name) + super(SubcommNotFoundError, self).__init__(message=message) + + +class SubcommAlreadyCreatedError(OcgDistError): + """Raised when a subcommunicator name already exists.""" + + def __init__(self, name): + message = "Subcommunicator '{}' already created.".format(name) + super(SubcommAlreadyCreatedError, self).__init__(message=message) + + +class CRSNotEquivalenError(OcgException): + """Raised when coordinate systems are not equivalent (not compatible for transform).""" + + def __init__(self, lhs, rhs): + msg = '{} is not equivalent to {}.'.format(lhs, rhs) + super(CRSNotEquivalenError, self).__init__(message=msg) + + +class DimensionMapError(OcgException): + """Raised when there is an issue with a dimension map entry.""" + + def __init__(self, entry_key, message): + msg = "Error with entry key '{}': {}".format(entry_key, message) + super(DimensionMapError, self).__init__(message=msg) + + +class VariableMissingMetadataError(OcgException): + """Raised when variable metadata cannot be found.""" + + def __init__(self, variable_name): + msg = 'Variable is missing metadata: {}'.format(variable_name) + super(VariableMissingMetadataError, self).__init__(message=msg) + + +class WrappedStateEvalTargetMissing(OcgException): + """Raised when attempting to retrieve the wrapped state of a field and no evaluation target is available.""" + + def __init__(self): + super(WrappedStateEvalTargetMissing, self).__init__(message='Target has no spatial information to evaluate.') + + +class NoTouching(OcgException): + """ + Raised when a subset geometry only touches an element. Used to prevent duplicate cells when creating a spatial + decomposition. + """ + omessage = "No touching allowed by spatial operation." + + +class SelfIntersectsRemovalError(OcgException): + """Raised when node removal to fix self-intersections causes a geometry validation error.""" + omessage = 'Polygon not valid after removing self-intersecting points.' diff --git a/src/ocgis/test/test_ocgis/test_spatial/test_grid_chunker.py b/src/ocgis/test/test_ocgis/test_spatial/test_grid_chunker.py index 3c09453c6..c017f267a 100644 --- a/src/ocgis/test/test_ocgis/test_spatial/test_grid_chunker.py +++ b/src/ocgis/test/test_ocgis/test_spatial/test_grid_chunker.py @@ -246,56 +246,62 @@ def _create_mRequestDataset_(): @attr('esmf', 'slow') def test_system_negative_values_in_spherical_grid(self): original_dir = os.getcwd() - xcn = np.arange(-10, 350, step=10, dtype=float) - xc = np.arange(0, 360, step=10, dtype=float) - yc = np.arange(-90, 100, step=10, dtype=float) - - xvn = Variable("lon", xcn, dimensions=["lon"]) - xv = Variable("lon", xc, dimensions=["lon"]) - yv = Variable("lat", yc, dimensions=["lat"]) - - gridn = Grid(x=xvn.copy(), y=yv.copy(), crs=Spherical()) - gridu = Grid(x=xv.copy(), y=yv.copy(), crs=Spherical()) - gridw = create_gridxy_global(5, with_bounds=False, crs=Spherical()) - grids = [gridn, gridu, gridw] - for ctr, (src, dst) in enumerate(itertools.product(grids, grids)): - os.chdir(self.current_dir_output) - gdirname = "grid-ctr-{}".format(ctr) - self.dprint(gdirname) - griddir = os.path.join(self.current_dir_output, gdirname) - os.mkdir(gdirname) - os.chdir(gdirname) - - srcgridname = "gridn.nc" - src.parent.write(srcgridname) - dstgridname = "grid.nc" - dst.parent.write(dstgridname) - - nchunks_dst = [(4, 1), (3, 1), (2, 1), (1, 1)] - for ctr, n in enumerate(nchunks_dst): - os.chdir(griddir) - dirname = 'ctr-{}'.format(ctr) - os.mkdir(dirname) - os.chdir(dirname) - wd = os.getcwd() - self.dprint("current chunks", n) - g = GridChunker(src, dst, nchunks_dst=n, genweights=True, paths={'wd': wd}, - esmf_kwargs={'regrid_method': 'BILINEAR'}) - if not g.is_one_chunk: - g.write_chunks() - g.create_merged_weight_file(os.path.join(griddir, "ctr-{}".format(ctr), "merged-weights.nc")) - else: - g.write_esmf_weights(os.path.join(griddir, srcgridname), - os.path.join(griddir, dstgridname), - os.path.join(griddir, "global-weights.nc")) + try: + xcn = np.arange(-10, 350, step=10, dtype=float) + xc = np.arange(0, 360, step=10, dtype=float) + yc = np.arange(-90, 100, step=10, dtype=float) + + xvn = Variable("lon", xcn, dimensions=["lon"]) + xv = Variable("lon", xc, dimensions=["lon"]) + yv = Variable("lat", yc, dimensions=["lat"]) + + gridn = Grid(x=xvn.copy(), y=yv.copy(), crs=Spherical()) + gridu = Grid(x=xv.copy(), y=yv.copy(), crs=Spherical()) + gridw = create_gridxy_global(5, with_bounds=False, crs=Spherical()) + grids = [gridn, gridu, gridw] + for ctr, (src, dst) in enumerate(itertools.product(grids, grids)): + os.chdir(self.current_dir_output) + gdirname = "grid-ctr-{}".format(ctr) + self.dprint(gdirname) + griddir = os.path.join(self.current_dir_output, gdirname) + os.mkdir(gdirname) + os.chdir(gdirname) + + srcgridname = "gridn.nc" + src.parent.write(srcgridname) + dstgridname = "grid.nc" + dst.parent.write(dstgridname) + + nchunks_dst = [ + (4, 1), + (3, 1), + (2, 1), + (1, 1) + ] + for ctr, n in enumerate(nchunks_dst): + os.chdir(griddir) + dirname = 'ctr-{}'.format(ctr) + os.mkdir(dirname) + os.chdir(dirname) + wd = os.getcwd() + self.dprint("current chunks", n) + g = GridChunker(src, dst, nchunks_dst=n, genweights=True, paths={'wd': wd}, + esmf_kwargs={'regrid_method': 'BILINEAR'}) + if not g.is_one_chunk: + g.write_chunks() + g.create_merged_weight_file(os.path.join(griddir, "ctr-{}".format(ctr), "merged-weights.nc")) + else: + g.write_esmf_weights(os.path.join(griddir, srcgridname), + os.path.join(griddir, dstgridname), + os.path.join(griddir, "global-weights.nc")) - os.chdir(griddir) - for ctr in range(0, len(nchunks_dst)-1): - src_filename = os.path.join(griddir, "ctr-{}".format(ctr), "merged-weights.nc") - dst_filename = os.path.join(griddir, "global-weights.nc") - self.assertWeightFilesEquivalent(src_filename, dst_filename) - - os.chdir(original_dir) + os.chdir(griddir) + for ctr in range(0, len(nchunks_dst)-1): + src_filename = os.path.join(griddir, "ctr-{}".format(ctr), "merged-weights.nc") + dst_filename = os.path.join(griddir, "global-weights.nc") + self.assertWeightFilesEquivalent(src_filename, dst_filename) + finally: + os.chdir(original_dir) def test_system_scrip_destination_splitting(self): """Test splitting a SCRIP destination grid.""" diff --git a/src/ocgis/test/test_ocgis/test_variable/test_geom.py b/src/ocgis/test/test_ocgis/test_variable/test_geom.py index c83b78284..99ad8d98e 100644 --- a/src/ocgis/test/test_ocgis/test_variable/test_geom.py +++ b/src/ocgis/test/test_ocgis/test_variable/test_geom.py @@ -25,10 +25,45 @@ from ocgis.variable.crs import WGS84, Spherical, Cartesian from ocgis.variable.dimension import Dimension from ocgis.variable.geom import GeometryVariable, GeometryProcessor, get_split_polygon_by_node_threshold, \ - GeometrySplitter + GeometrySplitter, do_remove_self_intersects_multi from ocgis.vmachine.mpi import OcgDist, MPI_RANK, variable_scatter, MPI_SIZE, variable_gather, MPI_COMM +class FixtureSelfIntersectingPolygon(object): + + @property + def fixture_self_intersecting_polygon_coords(self): + coords = [(0, 0), (0, 3), (2, 3), (2, 2), (1, 2), (1, 1), (2, 1), (2, 2), (3, 2), (3, 0), (0, 0)] + return coords + + @property + def fixture_self_intersecting_polygon_coords_first(self): + # Fixture has the first coordinate as the repeating one. + coords = [(2, 2), (3, 2), (3, 0), (0, 0), (0, 3), (2, 3), (2, 2), (1, 2), (1, 1), (2, 1), (2, 2)] + return coords + + +class TestGeom(TestBase, FixtureSelfIntersectingPolygon): + + def run_do_remove_self_intersects(self, fixture, debug=False): + poly = Polygon(fixture) + if debug: + gvar = GeometryVariable.from_shapely(poly) + gvar.write_vector('/tmp/vector.shp') + + new_poly = do_remove_self_intersects_multi(poly) + self.assertEqual(np.array(poly.exterior.coords).shape[0] - 1, np.array(new_poly.exterior.coords).shape[0]) + self.assertTrue(new_poly.is_valid) + if debug: + GeometryVariable.from_shapely(new_poly).write_vector('/tmp/new_vector.shp') + + def test_do_remove_self_intersects(self): + self.run_do_remove_self_intersects(self.fixture_self_intersecting_polygon_coords, debug=False) + + def test_do_remove_self_intersects_first(self): + self.run_do_remove_self_intersects(self.fixture_self_intersecting_polygon_coords_first, debug=False) + + class TestGeometryProcessor(AbstractTestInterface): def test_iter_intersection(self): @@ -159,7 +194,7 @@ def test_iter_interiors(self): self.assertEqual(actual, [(2.5, 10.5, 3.5, 15.5)]) -class TestGeometryVariable(AbstractTestInterface, FixturePolygonWithHole): +class TestGeometryVariable(AbstractTestInterface, FixturePolygonWithHole, FixtureSelfIntersectingPolygon): @staticmethod def get_geometryvariable_with_parent(): @@ -405,6 +440,15 @@ def test_convert_to_geometry_coordinates_polygons(self): for actual_geom, desired_geom in zip(actual.get_geometry_iterable(), geom.get_value().flat): self.assertEqual(actual_geom[1], desired_geom) + def test_convert_to_self_intersecting(self): + poly = Polygon(self.fixture_self_intersecting_polygon_coords) + gvar = GeometryVariable.from_shapely(poly) + without_si = gvar.convert_to(remove_self_intersects=True) + gvar2 = without_si.convert_to() + desired = do_remove_self_intersects_multi(poly) + self.assertPolygonSimilar(gvar2.v()[0], desired) + # gvar2.write_vector('/tmp/without_si.shp') + @attr('mpi') def test_create_ugid_global(self): ompi = OcgDist() diff --git a/src/ocgis/util/logging_ocgis.py b/src/ocgis/util/logging_ocgis.py index a46313082..ddac815af 100644 --- a/src/ocgis/util/logging_ocgis.py +++ b/src/ocgis/util/logging_ocgis.py @@ -4,8 +4,9 @@ import sys from warnings import warn -import ocgis import six + +import ocgis from ocgis import env, vm from ocgis.exc import OcgWarning @@ -101,7 +102,7 @@ def __call__(self, msg=None, logger=None, level=logging.INFO, alias=None, ugid=N else: if level == logging.WARN: wmsg = '{0}: {1}'.format(exc.__class__.__name__, str(exc)) - dest_logger.warn(wmsg) + dest_logger.warning(wmsg) else: dest_logger.exception(msg) raise exc diff --git a/src/ocgis/variable/geom.py b/src/ocgis/variable/geom.py index 65c0005d2..ab6972b1d 100644 --- a/src/ocgis/variable/geom.py +++ b/src/ocgis/variable/geom.py @@ -20,7 +20,7 @@ from ocgis.constants import KeywordArgument, HeaderName, VariableName, DimensionName, ConversionTarget, DriverKey, \ WrappedState, AttributeName, WrapAction from ocgis.environment import ogr -from ocgis.exc import EmptySubsetError, RequestableFeature, NoInteriorsError +from ocgis.exc import EmptySubsetError, RequestableFeature, NoInteriorsError, SelfIntersectsRemovalError from ocgis.spatial.base import AbstractSpatialVariable, create_split_polygons from ocgis.util.addict import Dict from ocgis.util.helpers import iter_array, get_trimmed_array_by_mask, get_swap_chain, find_index, \ @@ -34,7 +34,19 @@ CreateGeometryFromWkb, Geometry, wkbGeometryCollection, wkbPoint = ogr.CreateGeometryFromWkb, ogr.Geometry, \ ogr.wkbGeometryCollection, ogr.wkbPoint GEOM_TYPE_MAPPING = {'Polygon': Polygon, 'Point': Point, 'MultiPoint': MultiPoint, 'MultiPolygon': MultiPolygon} -_LOCAL_LOG = "spatial.base" +_LOCAL_LOG = "geom" + + +def count_interiors(geometry): + if not isinstance(geometry, Polygon) and not isinstance(geometry, MultiPolygon): + return 0 + try: + ret = len(geometry.interiors) + except AttributeError: + ret = 0 + for g in geometry: + ret += len(g.interiors) + return ret class GeometrySplitter(AbstractOcgisObject): @@ -48,13 +60,7 @@ def __init__(self, geometry): @property def interior_count(self): - try: - ret = len(self.geometry.interiors) - except AttributeError: - ret = 0 - for g in self.geometry: - ret += len(g.interiors) - return ret + return count_interiors(self.geometry) def create_split_vector_dict(self, interior): minx, miny, maxx, maxy = self.geometry.buffer(self._buffer_split).bounds @@ -361,6 +367,8 @@ def convert_to(self, target=ConversionTarget.GEOMETRY_COORDS, **kwargs): :keyword bool allow_splitting_excs: ``(=False)`` If ``True``, log and emit a warning when geometry processing errors are encountered during splitting operations. These exceptions may occur when holes/interiors are encountered or a node threshold is exceeded. + :keyword bool remove_self_intersects: ``(=False)`` If ``True``, attempt to remove self-intersections from + polygon objects. """ # TODO: IMPLEMENT: Line conversion. @@ -394,6 +402,7 @@ def convert_to(self, target=ConversionTarget.GEOMETRY_COORDS, **kwargs): to_crs = kwargs.pop('to_crs', None) start_index = kwargs.pop('start_index', 0) allow_splitting_excs = kwargs.pop('allow_splitting_excs', False) + remove_self_intersects = kwargs.pop('remove_self_intersects', False) assert len(kwargs) == 0 if to_crs is not None and not use_geometry_iterator: raise ValueError("'to_crs' only applies when using a geometry iterator") @@ -409,6 +418,7 @@ def convert_to(self, target=ConversionTarget.GEOMETRY_COORDS, **kwargs): # Problematic indices in the conversion. These may be removed if allow_splitting_excs is true. removed_indices = [] + has_z = False if target == ConversionTarget.GEOMETRY_COORDS: if use_geometry_iterator: has_z_itr = self._request_dataset.driver.get_variable_value(self, as_geometry_iterator=True) @@ -464,6 +474,15 @@ def convert_to(self, target=ConversionTarget.GEOMETRY_COORDS, **kwargs): geom = to_transform.get_value()[0] if geom.geom_type in polygon_types: + # Identify and remove self-intersections if requested. These can create virtual holes/interiors + # in polygon objects leading to issues with splitting for node reduction and hole removal. + if remove_self_intersects: + try: + geom = do_remove_self_intersects_multi(geom) + except SelfIntersectsRemovalError as e: + if allow_splitting_excs: + removed_indices.append(idx) + continue try: gsplitter = GeometrySplitter(geom) except NoInteriorsError: @@ -1414,3 +1433,58 @@ def geometryvariable_get_mask_from_intersects(gvar, geometry, use_spatial_index= ref_fill_mask[global_index[idx]] = bool_value return fill + + +def do_remove_self_intersects(poly, try_again=True): + if not isinstance(poly, Polygon): + exc = ValueError("only Polygons supported") + try: + raise exc + finally: + vm.abort(exc=exc) + if count_interiors(poly) > 0: + raise SelfIntersectsRemovalError('not supported when polygon has interiors') + + coords = np.array(poly.exterior.coords[0:-1]) + coords_to_remove = deque() + indices_to_remove = deque() + for idx in range(coords.shape[0]): + ctr = 0 + for idx2 in range(coords.shape[0]): + if np.all(coords[idx, :] == coords[idx2, :]): + ctr += 1 + if ctr > 1 and not any([np.all(c == coords[idx, :]) for c in coords_to_remove]): + coords_to_remove.append(coords[idx, :]) + indices_to_remove.append(idx) + + if len(coords_to_remove) == 0: + new_poly = poly + else: + new_poly_points = np.zeros((coords.shape[0]-len(coords_to_remove)+1, 2), dtype=float) + offset = 0 + for idx in range(coords.shape[0]): + if idx not in indices_to_remove: + new_poly_points[offset] = coords[idx, :] + offset += 1 + else: + indices_to_remove.remove(idx) + new_poly_points[-1, :] = new_poly_points[0, :] + new_poly = Polygon(new_poly_points) + + if not new_poly.is_valid: + new_poly = new_poly.buffer(0) + if not new_poly.is_valid: + if try_again: + new_poly = do_remove_self_intersects(new_poly, try_again=False) + else: + raise SelfIntersectsRemovalError() + return new_poly + + +def do_remove_self_intersects_multi(poly): + if isinstance(poly, MultiPolygon): + newpoly = [do_remove_self_intersects(p) for p in poly] + newpoly = MultiPolygon(newpoly) + else: + newpoly = do_remove_self_intersects(poly) + return newpoly