From 212e7780f3b9e9d839d675167130b9fa3ac77196 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C4=81nis=20Gailis?= Date: Thu, 8 Mar 2018 14:52:20 +0100 Subject: [PATCH] Subset long_polygons Original extents information is preserved in subset_spatial implementation if a tuple or a string giving four values is used for PolygonLike. Plotting extents are fixed accordingly. Closes #541 --- cate/core/opimpl.py | 65 ++++++++++++++++++++++++++++++---------- cate/core/types.py | 3 +- cate/ds/esa_cci_odp.py | 1 - cate/ops/plot.py | 13 ++++---- cate/ops/plot_helpers.py | 22 ++++++++++++++ cate/ops/subset.py | 1 - 6 files changed, 78 insertions(+), 27 deletions(-) diff --git a/cate/core/opimpl.py b/cate/core/opimpl.py index 1c0b2b25e..bee523802 100644 --- a/cate/core/opimpl.py +++ b/cate/core/opimpl.py @@ -29,6 +29,9 @@ from jdcal import jd2gcal from shapely.geometry import Point, box, LineString, Polygon +from .types import PolygonLike +from ..util.misc import to_list + def normalize_impl(ds: xr.Dataset) -> xr.Dataset: """ @@ -395,8 +398,37 @@ def _lat_inverted(lat: xr.DataArray) -> bool: return False +def get_extents(region: PolygonLike.TYPE): + """ + Get extents of a PolygonLike, handles explicitly given + coordinates. + + :param region: PolygonLike to introspect + :return: ([lon_min, lat_min, lon_max, lat_max], boolean_explicit_coords) + """ + explicit_coords = False + try: + maybe_rectangle = to_list(region, dtype=float) + if maybe_rectangle is not None: + if len(maybe_rectangle) == 4: + lon_min, lat_min, lon_max, lat_max = maybe_rectangle + explicit_coords = True + except: + # The polygon must be convertible, but it's complex + polygon = PolygonLike.convert(region) + if not polygon.is_valid: + # Heal polygon, see #506 and Shapely User Manual + region = polygon.buffer(0) + else: + region = polygon + # Get the bounding box + lon_min, lat_min, lon_max, lat_max = region.bounds + + return([lon_min, lat_min, lon_max, lat_max], explicit_coords) + + def subset_spatial_impl(ds: xr.Dataset, - region: Polygon, + region: PolygonLike.TYPE, mask: bool = True) -> xr.Dataset: """ Do a spatial subset of the dataset @@ -407,30 +439,30 @@ def subset_spatial_impl(ds: xr.Dataset, not the polygon itself be masked with NaN. :return: Subset dataset """ + # Validate input + try: + polygon = PolygonLike.convert(region) + except BaseException as e: + raise e - if not region.is_valid: - # Heal polygon, see #506 and Shapely User Manual - region = region.buffer(0) + extents, explicit_coords = get_extents(region) - # Get the bounding box - lon_min, lat_min, lon_max, lat_max = region.bounds + lon_min, lat_min, lon_max, lat_max = extents # Validate the bounding box if (not (-90 <= lat_min <= 90)) or \ (not (-90 <= lat_max <= 90)) or \ (not (-180 <= lon_min <= 180)) or \ (not (-180 <= lon_max <= 180)): - raise ValueError('Provided polygon extent outside of geospatial' + raise ValueError('Provided polygon extends outside of geospatial' ' bounds: latitude [-90;90], longitude [-180;180]') - simple_polygon = False - if region.equals(box(lon_min, lat_min, lon_max, lat_max)): - # Don't do the computationally intensive masking if the provided - # region is a simple box-polygon, for which there will be nothing to - # mask. - simple_polygon = True + # Don't do the computationally intensive masking if the provided + # region is a simple box-polygon, for which there will be nothing to + # mask. + simple_polygon = polygon.equals(box(lon_min, lat_min, lon_max, lat_max)) - crosses_antimeridian = _crosses_antimeridian(region) + crosses_antimeridian = (lon_min > lon_max) if explicit_coords else _crosses_antimeridian(polygon) lat_inverted = _lat_inverted(ds.lat) if lat_inverted: lat_index = slice(lat_max, lat_min) @@ -445,7 +477,8 @@ def subset_spatial_impl(ds: xr.Dataset, if crosses_antimeridian: # Shapely messes up longitudes if the polygon crosses the antimeridian - lon_min, lon_max = lon_max, lon_min + if not explicit_coords: + lon_min, lon_max = lon_max, lon_min # Can't perform a simple selection with slice, hence we have to # construct an appropriate longitude indexer for selection @@ -466,7 +499,7 @@ def subset_spatial_impl(ds: xr.Dataset, # dimension return retset - if not mask or simple_polygon: + if not mask or simple_polygon or explicit_coords: # The polygon doesn't cross the IDL, it is a simple box -> Use a simple slice lon_slice = slice(lon_min, lon_max) indexers = {'lat': lat_index, 'lon': lon_slice} diff --git a/cate/core/types.py b/cate/core/types.py index 6aef58c50..b9fb0945d 100644 --- a/cate/core/types.py +++ b/cate/core/types.py @@ -50,7 +50,6 @@ def some_op(file: PathLike.TYPE) -> bool: import xarray from shapely.errors import ShapelyError -from .opimpl import adjust_temporal_attrs_impl from ..util.misc import to_list, to_datetime_range, to_datetime from ..util.safe import safe_eval @@ -624,6 +623,8 @@ class DatasetLike(Like[xarray.Dataset]): @classmethod def convert(cls, value: Any) -> Optional[xarray.Dataset]: # Can be optional + from cate.core.opimpl import adjust_temporal_attrs_impl + if value is None: return None if isinstance(value, xarray.Dataset): diff --git a/cate/ds/esa_cci_odp.py b/cate/ds/esa_cci_odp.py index cc6fe44c9..55d64686a 100644 --- a/cate/ds/esa_cci_odp.py +++ b/cate/ds/esa_cci_odp.py @@ -619,7 +619,6 @@ def open_dataset(self, var_names: VarNamesLike.TYPE = None, protocol: str = None) -> Any: time_range = TimeRangeLike.convert(time_range) if time_range else None - region = PolygonLike.convert(region) if region else None var_names = VarNamesLike.convert(var_names) if var_names else None selected_file_list = self._find_files(time_range) diff --git a/cate/ops/plot.py b/cate/ops/plot.py index 71e598696..3085172d8 100644 --- a/cate/ops/plot.py +++ b/cate/ops/plot.py @@ -79,7 +79,7 @@ from cate.ops.plot_helpers import get_var_data from cate.ops.plot_helpers import in_notebook -from cate.ops.plot_helpers import check_bounding_box +from cate.ops.plot_helpers import handle_plot_polygon PLOT_FILE_EXTENSIONS = ['eps', 'jpeg', 'jpg', 'pdf', 'pgf', 'png', 'ps', 'raw', 'rgba', 'svg', @@ -159,12 +159,9 @@ def plot_map(ds: xr.Dataset, properties = DictLike.convert(properties) or {} extents = None - region = PolygonLike.convert(region) - if region: - lon_min, lat_min, lon_max, lat_max = region.bounds - if not check_bounding_box(lat_min, lat_max, lon_min, lon_max): - raise ValueError('Provided plot extents do not form a valid bounding box ' - 'within [-180.0,+180.0,-90.0,+90.0]') + bounds = handle_plot_polygon(region) + if bounds: + lon_min, lat_min, lon_max, lat_max = bounds extents = [lon_min, lon_max, lat_min, lat_max] # See http://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html# @@ -194,7 +191,7 @@ def plot_map(ds: xr.Dataset, figure = plt.figure(figsize=(8, 4)) ax = plt.axes(projection=proj) if extents: - ax.set_extent(extents) + ax.set_extent(extents, ccrs.PlateCarree()) else: ax.set_global() diff --git a/cate/ops/plot_helpers.py b/cate/ops/plot_helpers.py index 84e04452f..ba69a9dda 100644 --- a/cate/ops/plot_helpers.py +++ b/cate/ops/plot_helpers.py @@ -29,6 +29,28 @@ ========== """ +from cate.core.types import PolygonLike +from cate.core.opimpl import get_extents + + +def handle_plot_polygon(region: PolygonLike.TYPE = None): + """ + Return extents of the given PolygonLike. + + :param region: PolygonLike to introspect + :return: extents + """ + if region is None: + return None + + extents, explicit_coords = get_extents(region) + + lon_min, lat_min, lon_max, lat_max = extents + + if not check_bounding_box(lat_min, lat_max, lon_min, lon_max): + raise ValueError('Provided plot extents do not form a valid bounding box ' + 'within [-180.0,+180.0,-90.0,+90.0]') + return extents def check_bounding_box(lat_min: float, diff --git a/cate/ops/subset.py b/cate/ops/subset.py index 2b9e4da06..46b4dcd92 100644 --- a/cate/ops/subset.py +++ b/cate/ops/subset.py @@ -51,7 +51,6 @@ def subset_spatial(ds: xr.Dataset, :param mask: Should values falling in the bounding box of the polygon but not the polygon itself be masked with NaN. :return: Subset dataset """ - region = PolygonLike.convert(region) return adjust_spatial_attrs(subset_spatial_impl(ds, region, mask))