diff --git a/README.md b/README.md index 37e8be205..7059403a4 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![Build Status](https://github.com/pytroll/satpy/workflows/CI/badge.svg?branch=master)](https://github.com/pytroll/satpy/actions?query=workflow%3A%22CI%22) +[![Build Status](https://github.com/pytroll/pyresample/workflows/CI/badge.svg?branch=master)](https://github.com/pytroll/pyresample/actions?query=workflow%3A%22CI%22) [![Coverage Status](https://coveralls.io/repos/github/pytroll/pyresample/badge.svg?branch=master)](https://coveralls.io/github/pytroll/pyresample?branch=master) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3372769.svg)](https://doi.org/10.5281/zenodo.3372769) diff --git a/RELEASING.md b/RELEASING.md index f620a8486..acc72a605 100644 --- a/RELEASING.md +++ b/RELEASING.md @@ -6,7 +6,7 @@ 4. run `loghub` and update the `CHANGELOG.md` file: ``` -loghub pytroll/pyresample --token $LOGHUB_GITHUB_TOKEN -st v0.8.0 -plg bug "Bugs fixed" -plg enhancement "Features added" -plg documentation "Documentation changes" +loghub pytroll/pyresample --token $LOGHUB_GITHUB_TOKEN -st v0.8.0 -plg bug "Bugs fixed" -plg enhancement "Features added" -plg documentation "Documentation changes -plg backwards-incompatibility "Backward incompatible changes" -plg refactor "Refactoring" ``` This uses a `LOGHUB_GITHUB_TOKEN` environment variable. This must be created diff --git a/pyresample/_spatial_mp.py b/pyresample/_spatial_mp.py index 4a154562b..6ebbf5a64 100644 --- a/pyresample/_spatial_mp.py +++ b/pyresample/_spatial_mp.py @@ -22,9 +22,6 @@ import numpy as np import pyproj import multiprocessing as mp -import warnings - -from pyresample.utils import proj4_str_to_dict, is_pyproj2 try: @@ -113,17 +110,11 @@ class BaseProj(pyproj.Proj): """Helper class for easier backwards compatibility.""" def __init__(self, projparams=None, preserve_units=True, **kwargs): - if is_pyproj2(): - # have to have this because pyproj uses __new__ - # subclasses would fail when calling __init__ otherwise - super(BaseProj, self).__init__(projparams=projparams, - preserve_units=preserve_units, - **kwargs) - - def is_latlong(self): - if is_pyproj2(): - return self.crs.is_geographic - return super(BaseProj, self).is_latlong() + # have to have this because pyproj uses __new__ + # subclasses would fail when calling __init__ otherwise + super(BaseProj, self).__init__(projparams=projparams, + preserve_units=preserve_units, + **kwargs) class Proj(BaseProj): @@ -131,7 +122,7 @@ class Proj(BaseProj): def __call__(self, data1, data2, inverse=False, radians=False, errcheck=False, nprocs=1): - if self.is_latlong(): + if self.crs.is_geographic: return data1, data2 return super(Proj, self).__call__(data1, data2, inverse=inverse, radians=radians, errcheck=errcheck) @@ -146,7 +137,7 @@ def __init__(self, *args, **kwargs): def __call__(self, data1, data2, inverse=False, radians=False, errcheck=False, nprocs=2, chunk=None, schedule='guided'): - if self.is_latlong(): + if self.crs.is_geographic: return data1, data2 grid_shape = data1.shape diff --git a/pyresample/area_config.py b/pyresample/area_config.py index 23d10fbd7..eaeb46458 100644 --- a/pyresample/area_config.py +++ b/pyresample/area_config.py @@ -21,10 +21,13 @@ import io import warnings import pathlib +from typing import Any, Union import numpy as np import yaml from pyresample.utils import proj4_str_to_dict +from pyproj.crs import CRS, CRSError +from pyproj import Proj, Transformer try: @@ -457,7 +460,6 @@ def create_area_def(area_id, projection, width=None, height=None, area_extent=No they represent [projection x distance from center[0] to center[0]+dx, projection y distance from center[1] to center[1]+dy] """ - from pyresample._spatial_mp import Proj description = kwargs.pop('description', area_id) proj_id = kwargs.pop('proj_id', None) @@ -466,16 +468,16 @@ def create_area_def(area_id, projection, width=None, height=None, area_extent=No if isinstance(projection, dict) and 'EPSG' in projection: projection = "EPSG:{}".format(projection['EPSG']) - # Get a proj4_dict from either a proj4_dict or a proj4_string. - proj_dict = _get_proj_data(projection) try: - p = Proj(projection, preserve_units=True) - except RuntimeError: + crs = _get_proj_data(projection) + p = Proj(crs, preserve_units=True) + except (RuntimeError, CRSError): + # Assume that an invalid projection will be "fixed" by a dynamic area definition later return _make_area(area_id, description, proj_id, projection, shape, area_extent, **kwargs) # If no units are provided, try to get units used in proj_dict. If still none are provided, use meters. if units is None: - units = proj_dict.get('units', 'm' if not p.is_latlong() else 'degrees') + units = _get_proj_units(crs) # Allow height and width to be provided for more consistency across functions in pyresample. if height is not None or width is not None: @@ -490,14 +492,14 @@ def create_area_def(area_id, projection, width=None, height=None, area_extent=No area_extent = _verify_list('area_extent', area_extent, 4) # Converts from lat/lon to projection coordinates (x,y) if not in projection coordinates. Returns tuples. - center = _convert_units(center, 'center', units, p, proj_dict) - upper_left_extent = _convert_units(upper_left_extent, 'upper_left_extent', units, p, proj_dict) + center = _convert_units(center, 'center', units, p, crs) + upper_left_extent = _convert_units(upper_left_extent, 'upper_left_extent', units, p, crs) if area_extent is not None: # convert area extent, pass as (X, Y) area_extent_ll = area_extent[:2] area_extent_ur = area_extent[2:] - area_extent_ll = _convert_units(area_extent_ll, 'area_extent', units, p, proj_dict) - area_extent_ur = _convert_units(area_extent_ur, 'area_extent', units, p, proj_dict) + area_extent_ll = _convert_units(area_extent_ll, 'area_extent', units, p, crs) + area_extent_ur = _convert_units(area_extent_ur, 'area_extent', units, p, crs) area_extent = area_extent_ll + area_extent_ur # Fills in missing information to attempt to create an area definition. @@ -505,12 +507,19 @@ def create_area_def(area_id, projection, width=None, height=None, area_extent=No area_extent, shape, resolution = \ _extrapolate_information(area_extent, shape, center, radius, resolution, upper_left_extent, units, - p, proj_dict) + p, crs) return _make_area(area_id, description, proj_id, projection, shape, area_extent, resolution=resolution, **kwargs) -def _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **kwargs): +def _make_area( + area_id: str, + description: str, + proj_id: str, + projection: Union[dict, CRS], + shape: tuple, + area_extent: tuple, + **kwargs): """Handles the creation of an area definition for create_area_def.""" from pyresample.geometry import AreaDefinition from pyresample.geometry import DynamicAreaDefinition @@ -523,14 +532,14 @@ def _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **k if shape is not None: height, width = shape if None not in (area_extent, shape): - return AreaDefinition(area_id, description, proj_id, proj_dict, width, height, area_extent, **kwargs) + return AreaDefinition(area_id, description, proj_id, projection, width, height, area_extent, **kwargs) - return DynamicAreaDefinition(area_id=area_id, description=description, projection=proj_dict, width=width, + return DynamicAreaDefinition(area_id=area_id, description=description, projection=projection, width=width, height=height, area_extent=area_extent, rotation=kwargs.get('rotation'), resolution=resolution, optimize_projection=optimize_projection) -def _get_proj_data(projection): +def _get_proj_data(projection: Any) -> CRS: """Takes a proj4_dict or proj4_string and returns a proj4_dict and a Proj function. There is special handling for the "EPSG:XXXX" case where "XXXX" is an @@ -547,14 +556,20 @@ def _get_proj_data(projection): """ if isinstance(projection, dict) and 'EPSG' in projection: projection = "EPSG:{}".format(projection['EPSG']) + return CRS.from_user_input(projection) - if isinstance(projection, str): - proj_dict = proj4_str_to_dict(projection) - elif isinstance(projection, dict): - proj_dict = projection + +def _get_proj_units(crs): + if crs.is_geographic: + unit_name = 'degrees' else: - raise TypeError('Wrong type for projection: {0}. Expected dict or string.'.format(type(projection))) - return proj_dict + unit_name = crs.axis_info[0].unit_name + return { + 'metre': 'm', + 'meter': 'm', + 'kilometre': 'km', + 'kilometer': 'km', + }.get(unit_name, unit_name) def _sign(num): @@ -586,7 +601,10 @@ def _round_poles(center, units, p): return center -def _distance_from_center_forward(var, center, p): +def _distance_from_center_forward( + var: tuple, + center: tuple, + p: Proj): """Convert distances in degrees to projection units.""" # Interprets radius and resolution as distances between latitudes/longitudes. # Since the distance between longitudes and latitudes is not constant in @@ -597,7 +615,7 @@ def _distance_from_center_forward(var, center, p): center_as_angle = p(*center, inverse=True, errcheck=True) pole = 90 # If on a pole, use northern/southern latitude for both height and width. - if abs(abs(center_as_angle[1]) - pole) < 1e-8: + if abs(abs(center_as_angle[1]) - pole) < 1e-3: direction_of_poles = _sign(center_as_angle[1]) var = (center[1] - p(0, center_as_angle[1] - direction_of_poles * abs(var[0]), errcheck=True)[1], @@ -611,20 +629,25 @@ def _distance_from_center_forward(var, center, p): return var -def _convert_units(var, name, units, p, proj_dict, inverse=False, center=None): +def _convert_units( + var, + name: str, + units: str, + p: Proj, + crs: CRS, + inverse: bool = False, + center=None): """Converts units from lon/lat to projection coordinates (meters). If `inverse` it True then the inverse calculation is done. """ - from pyproj import transform - from pyresample._spatial_mp import Proj if var is None: return None if isinstance(var, DataArray): units = var.units var = tuple(var.data.tolist()) - if p.is_latlong() and not ('deg' == units or 'degrees' == units): + if crs.is_geographic and not ('deg' == units or 'degrees' == units): raise ValueError('latlon/latlong projection cannot take {0} as units: {1}'.format(units, name)) # Check if units are an angle. is_angle = ('deg' == units or 'degrees' == units) @@ -634,10 +657,11 @@ def _convert_units(var, name, units, p, proj_dict, inverse=False, center=None): if not is_angle: if units == 'meters' or units == 'metres': units = 'm' - if proj_dict.get('units', 'm') != units: - tmp_proj_dict = proj_dict.copy() + if _get_proj_units(crs) != units: + tmp_proj_dict = crs.to_dict() tmp_proj_dict['units'] = units - var = transform(Proj(tmp_proj_dict, preserve_units=True), p, *var) + transformer = Transformer.from_crs(tmp_proj_dict, p.crs) + var = transformer.transform(*var) if name == 'center': var = _round_poles(var, units, p) # Return either degrees or meters depending on if the inverse is true or not. @@ -646,7 +670,10 @@ def _convert_units(var, name, units, p, proj_dict, inverse=False, center=None): if is_angle and not inverse: if name in ('radius', 'resolution'): var = _distance_from_center_forward(var, center, p) - else: + elif not crs.is_geographic: + # only convert to meters + # this allows geographic projections to use coordinates outside + # normal lon/lat ranges (ex. -90/90) var = p(*var, errcheck=True) # Don't convert if inverse is False: Want meters. elif not is_angle and inverse: @@ -706,13 +733,13 @@ def _validate_variable(var, new_var, var_name, input_list): return new_var -def _extrapolate_information(area_extent, shape, center, radius, resolution, upper_left_extent, units, p, proj_dict): +def _extrapolate_information(area_extent, shape, center, radius, resolution, upper_left_extent, units, p, crs): """Attempts to find shape and area_extent based on data provided. Parameters are used in a specific order to determine area_extent and shape. The area_extent and shape are later used to create an `AreaDefinition`. Providing some parameters may have no effect if other parameters could be - to used determine area_extent and shape. The order of the parameters used + used to determine area_extent and shape. The order of the parameters used is: 1. area_extent @@ -731,7 +758,7 @@ def _extrapolate_information(area_extent, shape, center, radius, resolution, upp new_center = ((area_extent[2] + area_extent[0]) / 2, (area_extent[3] + area_extent[1]) / 2) center = _validate_variable(center, new_center, 'center', ['area_extent']) # If radius is given in an angle without center it will raise an exception, and to verify, it must be in meters. - radius = _convert_units(radius, 'radius', units, p, proj_dict, center=center) + radius = _convert_units(radius, 'radius', units, p, crs, center=center) new_radius = ((area_extent[2] - area_extent[0]) / 2, (area_extent[3] - area_extent[1]) / 2) radius = _validate_variable(radius, new_radius, 'radius', ['area_extent']) new_upper_left_extent = (area_extent[0], area_extent[3]) @@ -740,13 +767,13 @@ def _extrapolate_information(area_extent, shape, center, radius, resolution, upp # Output used below, but nowhere else is upper_left_extent made. Thus it should go as early as possible. elif None not in (upper_left_extent, center): # Function 1-B - radius = _convert_units(radius, 'radius', units, p, proj_dict, center=center) + radius = _convert_units(radius, 'radius', units, p, crs, center=center) new_radius = (center[0] - upper_left_extent[0], upper_left_extent[1] - center[1]) radius = _validate_variable(radius, new_radius, 'radius', ['upper_left_extent', 'center']) else: - radius = _convert_units(radius, 'radius', units, p, proj_dict, center=center) + radius = _convert_units(radius, 'radius', units, p, crs, center=center) # Convert resolution to meters if given as an angle. If center is not found, an exception is raised. - resolution = _convert_units(resolution, 'resolution', units, p, proj_dict, center=center) + resolution = _convert_units(resolution, 'resolution', units, p, crs, center=center) # Inputs unaffected by data below: area_extent is not an input. However, output is used below. if radius is not None and resolution is not None: # Function 2-A diff --git a/pyresample/ewa/_ll2cr.pyx b/pyresample/ewa/_ll2cr.pyx index 4ebdbb3bc..919e06e58 100644 --- a/pyresample/ewa/_ll2cr.pyx +++ b/pyresample/ewa/_ll2cr.pyx @@ -28,8 +28,6 @@ import numpy cimport cython cimport numpy -from pyresample._spatial_mp import BaseProj - # column and rows can only be doubles for now until the PROJ.4 is linked directly so float->double casting can be done # inside the loop ctypedef fused cr_dtype: diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 0c1babb02..b1f6357d2 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -37,8 +37,9 @@ from pyresample import CHUNK_SIZE from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj, Proj_MP from pyresample.boundary import AreaDefBoundary, Boundary, SimpleBoundary -from pyresample.utils import (proj4_str_to_dict, proj4_dict_to_str, - convert_proj_floats, proj4_radius_parameters, +from pyresample.utils import (proj4_dict_to_str, + proj4_radius_parameters, + get_geostationary_height, check_slice_orientation, load_cf_area) from pyresample.area_config import create_area_def @@ -47,10 +48,7 @@ except ImportError: DataArray = np.ndarray -try: - from pyproj import CRS -except ImportError: - CRS = None +from pyproj import CRS logger = getLogger(__name__) @@ -841,9 +839,9 @@ class DynamicAreaDefinition(object): description: The description of the area. projection: - The dictionary or string of projection parameters. Doesn't have to - be complete. If not complete, ``proj_info`` must be provided to - ``freeze`` to "fill in" any missing parameters. + The dictionary or string or CRS object of projection parameters. + Doesn't have to be complete. If not complete, ``proj_info`` must + be provided to ``freeze`` to "fill in" any missing parameters. width: x dimension in number of pixels, aka number of grid columns height: @@ -885,31 +883,16 @@ def __init__(self, area_id=None, description=None, projection=None, # check if non-dict projections are valid # dicts may be updated later if not isinstance(self._projection, dict): - Proj(projection) + CRS(projection) def _get_proj_dict(self): projection = self._projection - if CRS is not None: - try: - crs = CRS(projection) - except RuntimeError: - # could be incomplete dictionary - return projection - if hasattr(crs, 'to_dict'): - # pyproj 2.2+ - proj_dict = crs.to_dict() - else: - proj_dict = proj4_str_to_dict(crs.to_proj4()) - else: - if isinstance(projection, str): - proj_dict = proj4_str_to_dict(projection) - elif isinstance(projection, dict): - proj_dict = projection.copy() - else: - raise TypeError('Wrong type for projection: {0}. Expected ' - 'dict or string.'.format(type(projection))) - - return proj_dict + try: + crs = CRS(projection) + except RuntimeError: + # could be incomplete dictionary + return projection + return crs.to_dict() @property def pixel_size_x(self): @@ -1115,20 +1098,8 @@ def __init__(self, area_id, description, proj_id, projection, width, height, self.pixel_size_x = (area_extent[2] - area_extent[0]) / float(width) self.pixel_size_y = (area_extent[3] - area_extent[1]) / float(height) self._area_extent = tuple(area_extent) - if CRS is not None: - self.crs_wkt = CRS(projection).to_wkt() - self._proj_dict = None - self.crs = self._crs # see _crs property for details - else: - if isinstance(projection, str): - proj_dict = proj4_str_to_dict(projection) - elif isinstance(projection, dict): - # use the float-converted dict to pass to Proj - projection = convert_proj_floats(projection.items()) - proj_dict = projection - else: - raise TypeError('Wrong type for projection: {0}. Expected dict or string.'.format(type(projection))) - self._proj_dict = proj_dict + self.crs_wkt = CRS(projection).to_wkt() + self._proj_dict = None # Calculate area_extent in lon lat proj = Proj(projection) @@ -1154,7 +1125,7 @@ def __init__(self, area_id, description, proj_id, projection, width, height, self.dtype = dtype @property - def _crs(self): + def crs(self): """Wrap the `crs` property in a helper property. The :class:`pyproj.crs.CRS` object is not thread-safe. To avoid @@ -1162,32 +1133,16 @@ def _crs(self): is requested (the `self.crs` property). The alternative of storing it as a normal instance attribute could cause issues between threads. - For backwards compatibility, we only create the `.crs` property if - pyproj 2.0+ is installed. Users can then check - `hasattr(area_def, 'crs')` to easily support older versions of - pyresample and pyproj. - """ return CRS.from_wkt(self.crs_wkt) - @property - def _crs_or_proj_dict(self): - return self.crs if hasattr(self, 'crs') else self.proj_dict - - @property - def _wkt_or_proj_str(self): - return self.crs_wkt if hasattr(self, 'crs_wkt') else self.proj_str - @property def is_geostationary(self): """Whether this area is in a geostationary satellite projection or not.""" - if hasattr(self, 'crs'): - coord_operation = self.crs.coordinate_operation - if coord_operation is None: - return False - return 'geostationary' in coord_operation.method_name.lower() - else: - return self.proj_dict.get('proj') == 'geos' + coord_operation = self.crs.coordinate_operation + if coord_operation is None: + return False + return 'geostationary' in coord_operation.method_name.lower() @property def proj_dict(self): @@ -1196,12 +1151,8 @@ def proj_dict(self): This is no longer the preferred way of describing CRS information. Switch to the `crs` or `crs_wkt` properties for the most flexibility. """ - if self._proj_dict is None and hasattr(self, 'crs'): - if hasattr(self.crs, 'to_dict'): - # pyproj 2.2+ - self._proj_dict = self.crs.to_dict() - else: - self._proj_dict = proj4_str_to_dict(self.crs.to_proj4()) + if self._proj_dict is None: + self._proj_dict = self.crs.to_dict() return self._proj_dict @property @@ -1217,7 +1168,7 @@ def copy(self, **override_kwargs): kwargs = {'area_id': self.area_id, 'description': self.description, 'proj_id': self.proj_id, - 'projection': self._wkt_or_proj_str, + 'projection': self.crs_wkt, 'width': self.width, 'height': self.height, 'area_extent': self.area_extent, @@ -1560,11 +1511,11 @@ def to_cartopy_crs(self): self.area_extent[2], self.area_extent[1], self.area_extent[3]) - if hasattr(self, 'crs') and self.crs.to_epsg() is not None: + if self.crs.to_epsg() is not None: proj_params = "EPSG:{}".format(self.crs.to_epsg()) else: - proj_params = self.proj_str - if Proj(proj_params).is_latlong(): + proj_params = self.crs.to_proj4() + if self.crs.is_geographic: # Convert area extent from degrees to radians bounds = np.deg2rad(bounds) crs = from_proj(proj_params, bounds=bounds) @@ -1586,12 +1537,10 @@ def dump(self, filename=None): Returns: If file is None returns yaml str """ - if hasattr(self, 'crs') and self.crs.to_epsg() is not None: + if self.crs.to_epsg() is not None: proj_dict = {'EPSG': self.crs.to_epsg()} else: - proj_dict = self.proj_dict - # pyproj 2.0+ adds a '+type=crs' parameter - proj_dict.pop('type', None) + proj_dict = self.crs.to_dict() res = OrderedDict(description=self.description, projection=OrderedDict(proj_dict), @@ -1616,6 +1565,8 @@ def dump(self, filename=None): def create_areas_def_legacy(self): """Create area definition in legacy format.""" + warnings.warn("Pyresample's legacy areas file format is deprecated. " + "Use the 'YAML' format instead.") proj_dict = self.proj_dict proj_str = ','.join(["%s=%s" % (str(k), str(proj_dict[k])) for k in sorted(proj_dict.keys())]) @@ -1651,8 +1602,7 @@ def update_hash(self, the_hash=None): """Update a hash, or return a new one if needed.""" if the_hash is None: the_hash = hashlib.sha1() - crs_wkt = self._wkt_or_proj_str - the_hash.update(crs_wkt.encode('utf-8')) + the_hash.update(self.crs_wkt.encode('utf-8')) the_hash.update(np.array(self.shape)) the_hash.update(np.array(self.area_extent)) return the_hash @@ -1932,9 +1882,7 @@ def projection_y_coords(self): def outer_boundary_corners(self): """Return the lon,lat of the outer edges of the corner points.""" from pyresample.spherical_geometry import Coordinate - proj_def = self._wkt_or_proj_str - proj = Proj(proj_def) - + proj = Proj(self.crs) corner_lons, corner_lats = proj((self.area_extent[0], self.area_extent[2], self.area_extent[2], self.area_extent[0]), (self.area_extent[3], self.area_extent[3], @@ -2000,20 +1948,18 @@ def get_lonlats(self, nprocs=None, data_slice=None, cache=False, dtype=None, chu # but if the user provided nprocs then this doesn't make sense raise ValueError("Can't specify 'nprocs' and 'chunks' at the same time") - # Proj.4 definition of target area projection - proj_def = self._wkt_or_proj_str if hasattr(target_x, 'chunks'): # we are using dask arrays, map blocks to th from dask.array import map_blocks res = map_blocks(invproj, target_x, target_y, chunks=(target_x.chunks[0], target_x.chunks[1], 2), - new_axis=[2], proj_dict=proj_def).astype(dtype) + new_axis=[2], proj_dict=self.crs_wkt).astype(dtype) return res[:, :, 0], res[:, :, 1] if nprocs > 1: - target_proj = Proj_MP(proj_def) + target_proj = Proj_MP(self.crs) else: - target_proj = Proj(proj_def) + target_proj = Proj(self.crs) # Get corresponding longitude and latitude values lons, lats = target_proj(target_x, target_y, inverse=True, nprocs=nprocs) @@ -2060,8 +2006,8 @@ def get_area_slices(self, area_to_cover, shape_divisible_by=None): raise NotImplementedError('Only AreaDefinitions can be used') # Intersection only required for two different projections - proj_def_to_cover = area_to_cover._crs_or_proj_dict - proj_def = self._crs_or_proj_dict + proj_def_to_cover = area_to_cover.crs + proj_def = self.crs if proj_def_to_cover == proj_def: logger.debug('Projections for data and slice areas are' ' identical: %s', @@ -2127,9 +2073,8 @@ def __getitem__(self, key): (self.pixel_upper_left[0] + (xslice.stop - 0.5) * self.pixel_size_x), (self.pixel_upper_left[1] - (yslice.start - 0.5) * self.pixel_size_y)) - proj_def = self._wkt_or_proj_str new_area = AreaDefinition(self.area_id, self.description, - self.proj_id, proj_def, + self.proj_id, self.crs, total_cols, total_rows, new_area_extent) @@ -2166,8 +2111,7 @@ def geocentric_resolution(self, ellps='WGS84', radius=None): x, y = self.get_proj_vectors() mid_col_x = np.repeat(x[mid_col], y.size) mid_row_y = np.repeat(y[mid_row], x.size) - proj_def = self._crs_or_proj_dict - src = Proj(proj_def) + src = Proj(self.crs) if radius: dst = Proj("+proj=cart +a={} +b={}".format(radius, radius)) else: @@ -2222,11 +2166,11 @@ def _make_slice_divisible(sli, max_size, factor=2): def get_geostationary_angle_extent(geos_area): """Get the max earth (vs space) viewing angles in x and y.""" # get some projection parameters - proj_def = geos_area._crs_or_proj_dict - a, b = proj4_radius_parameters(proj_def) + a, b = proj4_radius_parameters(geos_area.crs) + h = get_geostationary_height(geos_area.crs) req = a / 1000.0 rp = b / 1000.0 - h = geos_area.proj_dict['h'] / 1000.0 + req + h = h / 1000.0 + req # compute some constants aeq = 1 - req ** 2 / (h ** 2) @@ -2247,6 +2191,7 @@ def get_geostationary_bounding_box(geos_area, nb_points=50): """ xmax, ymax = get_geostationary_angle_extent(geos_area) + h = get_geostationary_height(geos_area.crs) # generate points around the north hemisphere in satellite projection # make it a bit smaller so that we stay inside the valid area @@ -2255,13 +2200,13 @@ def get_geostationary_bounding_box(geos_area, nb_points=50): ll_x, ll_y, ur_x, ur_y = geos_area.area_extent - x *= geos_area.proj_dict['h'] - y *= geos_area.proj_dict['h'] + x *= h + y *= h x = np.clip(np.concatenate([x, x[::-1]]), min(ll_x, ur_x), max(ll_x, ur_x)) y = np.clip(np.concatenate([y, -y]), min(ll_y, ur_y), max(ll_y, ur_y)) - return Proj(**geos_area.proj_dict)(x, y, inverse=True) + return Proj(geos_area.crs)(x, y, inverse=True) def combine_area_extents_vertical(area1, area2): @@ -2287,13 +2232,12 @@ def combine_area_extents_vertical(area1, area2): def concatenate_area_defs(area1, area2, axis=0): """Append *area2* to *area1* and return the results.""" - different_items = (set(area1.proj_dict.items()) ^ - set(area2.proj_dict.items())) + crs_is_equal = area1.crs == area2.crs if axis == 0: same_size = area1.width == area2.width else: raise NotImplementedError('Only vertical contatenation is supported.') - if different_items or not same_size: + if not crs_is_equal or not same_size: raise IncompatibleAreas("Can't concatenate area definitions with " "different projections: " "{0} and {1}".format(area1, area2)) @@ -2305,7 +2249,7 @@ def concatenate_area_defs(area1, area2, axis=0): else: raise NotImplementedError('Only vertical contatenation is supported.') return AreaDefinition(area1.area_id, area1.description, area1.proj_id, - area1.proj_dict, x_size, y_size, + area1.crs, x_size, y_size, area_extent) @@ -2322,10 +2266,18 @@ def __init__(self, *definitions, **kwargs): super(StackedAreaDefinition, self).__init__(nprocs=nprocs) self.dtype = kwargs.get('dtype', np.float64) self.defs = [] - self.proj_dict = {} + self.crs_wkt = None for definition in definitions: self.append(definition) + @property + def crs(self): + return CRS.from_wkt(self.crs_wkt) + + @property + def proj_dict(self): + return self.crs.to_dict() + @property def width(self): """Return width of the area definition.""" @@ -2367,10 +2319,10 @@ def append(self, definition): if definition.height == 0: return if not self.defs: - self.proj_dict = definition.proj_dict - elif self.proj_dict != definition.proj_dict: + self.crs_wkt = definition.crs_wkt + elif self.crs != definition.crs: raise NotImplementedError('Cannot append areas:' - ' Proj.4 dict mismatch') + ' CRS mismatch') try: self.defs[-1] = concatenate_area_defs(self.defs[-1], definition) except (IncompatibleAreas, IndexError): diff --git a/pyresample/test/test_dask_ewa.py b/pyresample/test/test_dask_ewa.py index 54a3eaa59..9bd8774f1 100644 --- a/pyresample/test/test_dask_ewa.py +++ b/pyresample/test/test_dask_ewa.py @@ -25,10 +25,7 @@ import pytest import numpy as np -try: - from pyproj import CRS -except ImportError: - CRS = None +from pyproj import CRS da = pytest.importorskip("dask.array") xr = pytest.importorskip("xarray") @@ -114,9 +111,8 @@ def get_test_data(input_shape=(100, 50), output_shape=(200, 100), output_proj=No geo_dims = ('y', 'x') if input_dims else None swath_def = _get_test_swath_def(input_area_shape, chunk_size, geo_dims) ds1.attrs['area'] = swath_def - if CRS is not None: - crs = CRS.from_string('+proj=latlong +datum=WGS84 +ellps=WGS84') - ds1 = ds1.assign_coords(crs=crs) + crs = CRS.from_string('+proj=latlong +datum=WGS84 +ellps=WGS84') + ds1 = ds1.assign_coords(crs=crs) target_area = _get_test_target_area(output_shape, output_proj) return ds1, swath_def, target_area @@ -145,17 +141,15 @@ def _coord_and_crs_checks(new_data, target_area, has_bands=False): assert 'x' in new_data.coords if has_bands: assert 'bands' in new_data.coords - if CRS is not None: - assert 'crs' in new_data.coords - assert isinstance(new_data.coords['crs'].item(), CRS) - assert 'lcc' in new_data.coords['crs'].item().to_proj4() - assert new_data.coords['y'].attrs['units'] == 'meter' - assert new_data.coords['x'].attrs['units'] == 'meter' - if hasattr(target_area, 'crs'): - assert target_area.crs is new_data.coords['crs'].item() - if has_bands: - np.testing.assert_equal(new_data.coords['bands'].values, - ['R', 'G', 'B']) + assert 'crs' in new_data.coords + assert isinstance(new_data.coords['crs'].item(), CRS) + assert 'lcc' in new_data.coords['crs'].item().to_proj4() + assert new_data.coords['y'].attrs['units'] == 'meter' + assert new_data.coords['x'].attrs['units'] == 'meter' + assert target_area.crs == new_data.coords['crs'].item() + if has_bands: + np.testing.assert_equal(new_data.coords['bands'].values, + ['R', 'G', 'B']) def _get_num_chunks(source_swath, resampler_class): diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py index 58e5da322..018dae868 100644 --- a/pyresample/test/test_geometry.py +++ b/pyresample/test/test_geometry.py @@ -35,6 +35,8 @@ from unittest.mock import MagicMock, patch import unittest +import pyproj +import pytest try: from pyproj import CRS @@ -133,10 +135,7 @@ def test_cartopy_crs(self): self.assertEqual(crs.threshold, thresh_exp) # EPSG projection - projections = ['+init=EPSG:6932'] - if utils.is_pyproj2(): - projections.append('EPSG:6932') - + projections = ['+init=EPSG:6932', 'EPSG:6932'] for projection in projections: area = geometry.AreaDefinition( area_id='ease-sh-2.0', @@ -205,9 +204,10 @@ def test_dump(self): expected['projection'][proj_key]) # EPSG - projections = {'+init=epsg:3006': 'init: epsg:3006'} - if utils.is_pyproj2(): - projections['EPSG:3006'] = 'EPSG: 3006' + projections = { + '+init=epsg:3006': 'init: epsg:3006', + 'EPSG:3006': 'EPSG: 3006', + } for projection, epsg_yaml in projections.items(): area_def = geometry.AreaDefinition('baws300_sweref99tm', 'BAWS, 300m resolution, sweref99tm', @@ -272,9 +272,10 @@ def test_parse_area_file(self): self.assertEqual(area_def, expected) # EPSG - projections = {'+init=epsg:3006': 'init: epsg:3006'} - if utils.is_pyproj2(): - projections['EPSG:3006'] = 'EPSG: 3006' + projections = { + '+init=epsg:3006': 'init: epsg:3006', + 'EPSG:3006': 'EPSG: 3006', + } for projection, epsg_yaml in projections.items(): expected = geometry.AreaDefinition('baws300_sweref99tm', 'BAWS, 300m resolution, sweref99tm', 'sweref99tm', @@ -1216,9 +1217,7 @@ def test_get_area_slices(self): self.assertEqual(slice(0, y_size, None), slice_y) # totally different area - projections = [{"init": 'EPSG:4326'}] - if utils.is_pyproj2(): - projections.append('EPSG:4326') + projections = [{"init": 'EPSG:4326'}, 'EPSG:4326'] for projection in projections: area_to_cover = geometry.AreaDefinition( 'epsg4326', 'Global equal latitude/longitude grid for global sphere', @@ -1306,18 +1305,13 @@ def test_proj_str(self): ) # EPSG - if utils.is_pyproj2(): - # With pyproj 2.0+ we expand EPSG to full parameter list - full_proj = ('+datum=WGS84 +lat_0=-90 +lon_0=0 +no_defs ' - '+proj=laea +type=crs +units=m +x_0=0 +y_0=0') - projections = [ - ('+init=EPSG:6932', full_proj), - ('EPSG:6932', full_proj) - ] - else: - projections = [ - ('+init=EPSG:6932', '+init=EPSG:6932'), - ] + # With pyproj 2.0+ we expand EPSG to full parameter list + full_proj = ('+datum=WGS84 +lat_0=-90 +lon_0=0 +no_defs ' + '+proj=laea +type=crs +units=m +x_0=0 +y_0=0') + projections = [ + ('+init=EPSG:6932', full_proj), + ('EPSG:6932', full_proj) + ] for projection, expected_proj in projections: area = geometry.AreaDefinition( area_id='ease-sh-2.0', @@ -1328,37 +1322,36 @@ def test_proj_str(self): area_extent=[-40000., -40000., 40000., 40000.]) self.assertEqual(area.proj_str, expected_proj) - if utils.is_pyproj2(): - # CRS with towgs84 in it - # we remove towgs84 if they are all 0s - projection = {'proj': 'laea', 'lat_0': 52, 'lon_0': 10, 'x_0': 4321000, 'y_0': 3210000, - 'ellps': 'GRS80', 'towgs84': '0,0,0,0,0,0,0', 'units': 'm', 'no_defs': True} - area = geometry.AreaDefinition( - area_id='test_towgs84', - description='', - proj_id='', - projection=projection, - width=123, height=123, - area_extent=[-40000., -40000., 40000., 40000.]) - self.assertEqual(area.proj_str, - '+ellps=GRS80 +lat_0=52 +lon_0=10 +no_defs +proj=laea ' - # '+towgs84=0.0,0.0,0.0,0.0,0.0,0.0,0.0 ' - '+type=crs +units=m ' - '+x_0=4321000 +y_0=3210000') - projection = {'proj': 'laea', 'lat_0': 52, 'lon_0': 10, 'x_0': 4321000, 'y_0': 3210000, - 'ellps': 'GRS80', 'towgs84': '0,5,0,0,0,0,0', 'units': 'm', 'no_defs': True} - area = geometry.AreaDefinition( - area_id='test_towgs84', - description='', - proj_id='', - projection=projection, - width=123, height=123, - area_extent=[-40000., -40000., 40000., 40000.]) - self.assertEqual(area.proj_str, - '+ellps=GRS80 +lat_0=52 +lon_0=10 +no_defs +proj=laea ' - '+towgs84=0.0,5.0,0.0,0.0,0.0,0.0,0.0 ' - '+type=crs +units=m ' - '+x_0=4321000 +y_0=3210000') + # CRS with towgs84 in it + # we remove towgs84 if they are all 0s + projection = {'proj': 'laea', 'lat_0': 52, 'lon_0': 10, 'x_0': 4321000, 'y_0': 3210000, + 'ellps': 'GRS80', 'towgs84': '0,0,0,0,0,0,0', 'units': 'm', 'no_defs': True} + area = geometry.AreaDefinition( + area_id='test_towgs84', + description='', + proj_id='', + projection=projection, + width=123, height=123, + area_extent=[-40000., -40000., 40000., 40000.]) + self.assertEqual(area.proj_str, + '+ellps=GRS80 +lat_0=52 +lon_0=10 +no_defs +proj=laea ' + # '+towgs84=0.0,0.0,0.0,0.0,0.0,0.0,0.0 ' + '+type=crs +units=m ' + '+x_0=4321000 +y_0=3210000') + projection = {'proj': 'laea', 'lat_0': 52, 'lon_0': 10, 'x_0': 4321000, 'y_0': 3210000, + 'ellps': 'GRS80', 'towgs84': '0,5,0,0,0,0,0', 'units': 'm', 'no_defs': True} + area = geometry.AreaDefinition( + area_id='test_towgs84', + description='', + proj_id='', + projection=projection, + width=123, height=123, + area_extent=[-40000., -40000., 40000., 40000.]) + self.assertEqual(area.proj_str, + '+ellps=GRS80 +lat_0=52 +lon_0=10 +no_defs +proj=laea ' + '+towgs84=0.0,5.0,0.0,0.0,0.0,0.0,0.0 ' + '+type=crs +units=m ' + '+x_0=4321000 +y_0=3210000') def test_striding(self): """Test striding AreaDefinitions.""" @@ -1539,7 +1532,6 @@ def test_area_def_init_projection(self): # self.assertEqual(crs, area_def.crs) def test_areadef_immutable(self): - import pytest area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', {'a': '6378144.0', 'b': '6356759.0', @@ -1830,7 +1822,6 @@ def test_get_edge_lonlats(self): def test_compute_optimal_bb(self): """Test computing the bb area.""" - from pyresample.utils import is_pyproj2 import xarray as xr nplats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469], [80.84000396728516, 60.74200439453125, 34.08500289916992], @@ -1850,11 +1841,10 @@ def test_compute_optimal_bb(self): proj_dict = {'gamma': 0.0, 'lonc': -11.391744043133668, 'ellps': 'WGS84', 'proj': 'omerc', 'alpha': 9.185764390923012, 'lat_0': -0.2821013754097188} - if is_pyproj2(): - # pyproj2 adds some extra defaults - proj_dict.update({'x_0': 0, 'y_0': 0, 'units': 'm', - 'k': 1, 'gamma': 0, - 'no_defs': None, 'type': 'crs'}) + # pyproj2 adds some extra defaults + proj_dict.update({'x_0': 0, 'y_0': 0, 'units': 'm', + 'k': 1, 'gamma': 0, + 'no_defs': None, 'type': 'crs'}) assert_np_dict_allclose(res.proj_dict, proj_dict) self.assertEqual(res.shape, (6, 3)) @@ -1867,11 +1857,10 @@ def test_compute_optimal_bb(self): proj_dict = {'gamma': 0.0, 'lonc': -11.391744043133668, 'ellps': 'WGS84', 'proj': 'omerc', 'alpha': 9.185764390923012, 'lat_0': -0.2821013754097188} - if is_pyproj2(): - # pyproj2 adds some extra defaults - proj_dict.update({'x_0': 0, 'y_0': 0, 'units': 'm', - 'k': 1, 'gamma': 0, - 'no_defs': None, 'type': 'crs'}) + # pyproj2 adds some extra defaults + proj_dict.update({'x_0': 0, 'y_0': 0, 'units': 'm', + 'k': 1, 'gamma': 0, + 'no_defs': None, 'type': 'crs'}) assert_np_dict_allclose(res.proj_dict, proj_dict) self.assertEqual(res.shape, (6, 3)) @@ -1948,7 +1937,7 @@ def test_swath_def_geocentric_resolution(self): self.assertRaises(RuntimeError, sd.geocentric_resolution) -class TestStackedAreaDefinition(unittest.TestCase): +class TestStackedAreaDefinition: """Test the StackedAreaDefition.""" def test_append(self): @@ -1972,13 +1961,13 @@ def test_append(self): ) adef = geometry.StackedAreaDefinition(area1, area2) - self.assertEqual(len(adef.defs), 1) - self.assertTupleEqual(adef.defs[0].area_extent, - (3738502.0095458371, 4179561.259167064, - -1830246.0673044831, 3251436.5796920112)) + assert len(adef.defs) == 1 + assert adef.defs[0].area_extent == (3738502.0095458371, + 4179561.259167064, + -1830246.0673044831, + 3251436.5796920112) # same - area3 = geometry.AreaDefinition("area3", 'area3', "geosmsg", {'a': '6378169.0', 'b': '6356583.8', 'h': '35785831.0', 'lon_0': '0.0', @@ -1987,12 +1976,12 @@ def test_append(self): (3738502.0095458371, 3251436.5796920112, -1830246.0673044831, 2787374.2399544837)) adef.append(area3) - self.assertEqual(len(adef.defs), 1) - self.assertTupleEqual(adef.defs[0].area_extent, - (3738502.0095458371, 4179561.259167064, - -1830246.0673044831, 2787374.2399544837)) - - self.assertIsInstance(adef.squeeze(), geometry.AreaDefinition) + assert len(adef.defs) == 1 + assert adef.defs[0].area_extent == (3738502.0095458371, + 4179561.259167064, + -1830246.0673044831, + 2787374.2399544837) + assert isinstance(adef.squeeze(), geometry.AreaDefinition) # transition area4 = geometry.AreaDefinition("area4", 'area4', "geosmsg", @@ -2004,24 +1993,26 @@ def test_append(self): -1000.3358822065015, 2323311.9002169576)) adef.append(area4) - self.assertEqual(len(adef.defs), 2) - self.assertTupleEqual(adef.defs[-1].area_extent, - (5567747.7409681147, 2787374.2399544837, - -1000.3358822065015, 2323311.9002169576)) + assert len(adef.defs) == 2 + assert adef.defs[-1].area_extent == (5567747.7409681147, + 2787374.2399544837, + -1000.3358822065015, + 2323311.9002169576) - self.assertEqual(adef.height, 4 * 464) - self.assertIsInstance(adef.squeeze(), geometry.StackedAreaDefinition) + assert adef.height == 4 * 464 + assert isinstance(adef.squeeze(), geometry.StackedAreaDefinition) adef2 = geometry.StackedAreaDefinition() - self.assertEqual(len(adef2.defs), 0) + assert len(adef2.defs) == 0 adef2.append(adef) - self.assertEqual(len(adef2.defs), 2) - self.assertTupleEqual(adef2.defs[-1].area_extent, - (5567747.7409681147, 2787374.2399544837, - -1000.3358822065015, 2323311.9002169576)) + assert len(adef2.defs) == 2 + assert adef2.defs[-1].area_extent == (5567747.7409681147, + 2787374.2399544837, + -1000.3358822065015, + 2323311.9002169576) - self.assertEqual(adef2.height, 4 * 464) + assert adef2.height == 4 * 464 def test_get_lonlats(self): """Test get_lonlats on StackedAreaDefinition.""" @@ -2043,7 +2034,7 @@ def test_get_lonlats(self): -1000.3358822065015, 2323311.9002169576)) final_area = geometry.StackedAreaDefinition(area3, area4) - self.assertEqual(len(final_area.defs), 2) + assert len(final_area.defs) == 2 lons, lats = final_area.get_lonlats() lons0, lats0 = final_area.defs[0].get_lonlats() lons1, lats1 = final_area.defs[1].get_lonlats() @@ -2073,22 +2064,22 @@ def test_combine_area_extents(self): area2 = MagicMock() area2.area_extent = (1, 6, 3, 2) res = combine_area_extents_vertical(area1, area2) - self.assertListEqual(res, [1, 6, 3, 4]) + assert res == [1, 6, 3, 4] area1 = MagicMock() area1.area_extent = (1, 2, 3, 4) area2 = MagicMock() area2.area_extent = (1, 4, 3, 6) res = combine_area_extents_vertical(area1, area2) - self.assertListEqual(res, [1, 2, 3, 6]) + assert res == [1, 2, 3, 6] # Non contiguous area extends shouldn't be combinable area1 = MagicMock() area1.area_extent = (1, 2, 3, 4) area2 = MagicMock() area2.area_extent = (1, 5, 3, 7) - self.assertRaises(IncompatibleAreas, - combine_area_extents_vertical, area1, area2) + pytest.raises(IncompatibleAreas, combine_area_extents_vertical, + area1, area2) def test_append_area_defs_fail(self): """Fail appending areas.""" @@ -2101,8 +2092,7 @@ def test_append_area_defs_fail(self): area2.width = 4 area2.height = 6 # res = combine_area_extents_vertical(area1, area2) - self.assertRaises(IncompatibleAreas, - concatenate_area_defs, area1, area2) + pytest.raises(IncompatibleAreas, concatenate_area_defs, area1, area2) @patch('pyresample.geometry.AreaDefinition') def test_append_area_defs(self, adef): @@ -2110,13 +2100,13 @@ def test_append_area_defs(self, adef): x_size = random.randrange(6425) area1 = MagicMock() area1.area_extent = (1, 2, 3, 4) - area1.proj_dict = {"proj": 'A'} + area1.crs = 'some_crs' area1.height = random.randrange(6425) area1.width = x_size area2 = MagicMock() area2.area_extent = (1, 4, 3, 6) - area2.proj_dict = {"proj": 'A'} + area2.crs = 'some_crs' area2.height = random.randrange(6425) area2.width = x_size @@ -2124,111 +2114,179 @@ def test_append_area_defs(self, adef): area_extent = [1, 2, 3, 6] y_size = area1.height + area2.height adef.assert_called_once_with(area1.area_id, area1.description, area1.proj_id, - area1.proj_dict, area1.width, y_size, area_extent) - - def test_create_area_def(self): + area1.crs, area1.width, y_size, area_extent) + + @staticmethod + def _compare_area_defs(actual, expected): + actual_dict = actual.proj_dict + if 'EPSG' in actual_dict or 'init' in actual_dict: + # Use formal definition of EPSG projections to make them comparable to the base definition + proj_def = pyproj.Proj(actual_dict).definition_string().strip() + actual = actual.copy(projection=proj_def) + + # Remove extra attributes from the formal definition + # for key in ['x_0', 'y_0', 'no_defs', 'b', 'init']: + # actual.proj_dict.pop(key, None) + + assert actual == expected + + @pytest.mark.parametrize( + 'projection', + [ + {'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'}, + '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m', + '+init=EPSG:3409', + 'EPSG:3409', + ]) + @pytest.mark.parametrize( + 'center', + [ + [0, 0], + 'a', + (1, 2, 3), + ]) + @pytest.mark.parametrize('units', ['meters', 'degrees']) + def test_create_area_def_base_combinations(self, projection, center, units): """Test create_area_def and the four sub-methods that call it in AreaDefinition.""" from pyresample.geometry import AreaDefinition - from pyresample.geometry import DynamicAreaDefinition - from pyresample.area_config import DataArray from pyresample.area_config import create_area_def as cad - from pyresample import utils - import pyproj area_id = 'ease_sh' description = 'Antarctic EASE grid' - projection_list = [{'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'}, - '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m', - '+init=EPSG:3409'] - if utils.is_pyproj2(): - projection_list.append('EPSG:3409') proj_id = 'ease_sh' shape = (425, 850) - upper_left_extent = (-5326849.0625, 5326849.0625) - center_list = [[0, 0], 'a', (1, 2, 3)] area_extent = (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625) - resolution = (12533.7625, 25067.525) - radius = [5326849.0625, 5326849.0625] - units_list = ['meters', 'degrees'] - base_def = AreaDefinition(area_id, description, '', projection_list[0], shape[1], shape[0], area_extent) + base_def = AreaDefinition( + area_id, description, '', + {'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'}, + shape[1], shape[0], area_extent) # Tests that incorrect lists do not create an area definition, that both projection strings and # dicts are accepted, and that degrees and meters both create the same area definition. # area_list used to check that areas are all correct at the end. - area_list = [] - from itertools import product - for projection, units, center in product(projection_list, units_list, center_list): - # essentials = center, radius, upper_left_extent, resolution, shape. - if 'm' in units: - # Meters. - essentials = [[0, 0], [5326849.0625, 5326849.0625], (-5326849.0625, 5326849.0625), - (12533.7625, 25067.525), (425, 850)] - else: - # Degrees. - essentials = [(0.0, -90.0), 49.4217406986, (-45.0, -17.516001139327766), - (0.11271481862984278, 0.22542974631297721), (425, 850)] - # If center is valid, use it. - if len(center) == 2: - center = essentials[0] - try: - area_list.append(cad(area_id, projection, proj_id=proj_id, upper_left_extent=essentials[2], - center=center, shape=essentials[4], resolution=essentials[3], - radius=essentials[1], description=description, units=units, rotation=45)) - except ValueError: - pass - self.assertEqual(len(area_list), 8 if utils.is_pyproj2() else 6) + # essentials = center, radius, upper_left_extent, resolution, shape. + if 'm' in units: + # Meters. + essentials = [[0, 0], [5326849.0625, 5326849.0625], (-5326849.0625, 5326849.0625), + (12533.7625, 25067.525), (425, 850)] + else: + # Degrees. + essentials = [(0.0, -90.0), 49.4217406986, (-45.0, -17.516001139327766), + (0.11271481862984278, 0.22542974631297721), (425, 850)] + # If center is valid, use it. + if len(center) == 2: + center = essentials[0] + + args = (area_id, projection) + kwargs = dict( + proj_id=proj_id, + upper_left_extent=essentials[2], + center=center, + shape=essentials[4], + resolution=essentials[3], + radius=essentials[1], + description=description, + units=units, + rotation=45, + ) + + should_fail = isinstance(center, str) or len(center) != 2 + if should_fail: + pytest.raises(ValueError, cad, *args, **kwargs) + else: + area_def = cad(*args, **kwargs) + self._compare_area_defs(area_def, base_def) + + def test_create_area_def_extra_combinations(self): + """Test extra combinations of create_area_def parameters.""" + from xarray import DataArray + from pyresample import create_area_def as cad + from pyresample.geometry import AreaDefinition + + projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m' + area_id = 'ease_sh' + description = 'Antarctic EASE grid' + shape = (425, 850) + upper_left_extent = (-5326849.0625, 5326849.0625) + area_extent = (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625) + resolution = (12533.7625, 25067.525) + radius = [5326849.0625, 5326849.0625] + base_def = AreaDefinition( + area_id, description, '', + {'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'}, + shape[1], shape[0], area_extent) # Tests that specifying units through xarrays works. - area_list.append(cad(area_id, projection_list[1], shape=shape, - area_extent=DataArray((-135.0, -17.516001139327766, - 45.0, -17.516001139327766), - attrs={'units': 'degrees'}))) + area_def = cad(area_id, projection, shape=shape, + area_extent=DataArray( + (-135.0, -17.516001139327766, 45.0, -17.516001139327766), + attrs={'units': 'degrees'})) + self._compare_area_defs(area_def, base_def) + # Tests area functions 1-A and 2-A. - area_list.append(cad(area_id, projection_list[1], resolution=resolution, area_extent=area_extent)) + area_def = cad(area_id, projection, resolution=resolution, area_extent=area_extent) + self._compare_area_defs(area_def, base_def) + # Tests area function 1-B. Also test that DynamicAreaDefinition arguments don't crash AreaDefinition. - area_list.append(cad(area_id, projection_list[1], shape=shape, center=center_list[0], - upper_left_extent=upper_left_extent, optimize_projection=None)) + area_def = cad(area_id, projection, shape=shape, center=[0, 0], + upper_left_extent=upper_left_extent, optimize_projection=None) + self._compare_area_defs(area_def, base_def) + # Tests area function 1-C. - area_list.append(cad(area_id, projection_list[1], shape=shape, center=center_list[0], radius=radius)) + area_def = cad(area_id, projection, shape=shape, center=[0, 0], + radius=radius) + self._compare_area_defs(area_def, base_def) + # Tests area function 1-D. - area_list.append(cad(area_id, projection_list[1], shape=shape, - radius=radius, upper_left_extent=upper_left_extent)) + area_def = cad(area_id, projection, shape=shape, + radius=radius, upper_left_extent=upper_left_extent) + self._compare_area_defs(area_def, base_def) + # Tests all 4 user cases. - area_list.append(AreaDefinition.from_extent(area_id, projection_list[1], shape, area_extent)) - area_list.append(AreaDefinition.from_circle(area_id, projection_list[1], center_list[0], radius, - resolution=resolution)) - area_list.append(AreaDefinition.from_area_of_interest(area_id, projection_list[1], shape, center_list[0], - resolution)) - area_list.append(AreaDefinition.from_ul_corner(area_id, projection_list[1], shape, upper_left_extent, - resolution)) - # Tests non-poles using degrees and mercator. - area_def = cad(area_id, '+a=6371228.0 +units=m +lon_0=0 +proj=merc +lat_0=0', - center=(0, 0), radius=45, resolution=(1, 0.9999291722135637), units='degrees') - self.assertTrue(isinstance(area_def, AreaDefinition)) - self.assertTrue(np.allclose(area_def.area_extent, (-5003950.7698, -5615432.0761, 5003950.7698, 5615432.0761))) - self.assertEqual(area_def.shape, (101, 90)) - # Checks every area definition made - for area_def in area_list: - if 'EPSG' in area_def.proj_dict or 'init' in area_def.proj_dict: - # Use formal definition of EPSG projections to make them comparable to the base definition - proj_def = pyproj.Proj(area_def.proj_str).definition_string().strip() - area_def = area_def.copy(projection=proj_def) - - # Remove extra attributes from the formal definition - if 'R' in area_def.proj_dict: - # pyproj < 2 - area_def.proj_dict['a'] = area_def.proj_dict.pop('R') - for key in ['x_0', 'y_0', 'no_defs', 'b', 'init']: - area_def.proj_dict.pop(key, None) - - self.assertEqual(area_def, base_def) - - # Makes sure if shape or area_extent is found/given, a DynamicAreaDefinition is made. - self.assertTrue(isinstance(cad(area_id, projection_list[1], shape=shape), DynamicAreaDefinition)) - self.assertTrue(isinstance(cad(area_id, projection_list[1], area_extent=area_extent), DynamicAreaDefinition)) + area_def = AreaDefinition.from_extent(area_id, projection, shape, area_extent) + self._compare_area_defs(area_def, base_def) + + area_def = AreaDefinition.from_circle(area_id, projection, [0, 0], radius, + resolution=resolution) + self._compare_area_defs(area_def, base_def) + area_def = AreaDefinition.from_area_of_interest(area_id, projection, + shape, [0, 0], + resolution) + self._compare_area_defs(area_def, base_def) + area_def = AreaDefinition.from_ul_corner(area_id, projection, shape, + upper_left_extent, + resolution) + self._compare_area_defs(area_def, base_def) + + def test_create_area_def_nonpole_center(self): + """Test that a non-pole center can be used.""" + from pyresample import create_area_def as cad + from pyresample.geometry import AreaDefinition + area_def = cad('ease_sh', '+a=6371228.0 +units=m +lon_0=0 +proj=merc +lat_0=0', + center=(0, 0), radius=45, + resolution=(1, 0.9999291722135637), + units='degrees') + assert isinstance(area_def, AreaDefinition) + np.testing.assert_allclose(area_def.area_extent, (-5003950.7698, -5615432.0761, 5003950.7698, 5615432.0761)) + assert area_def.shape == (101, 90) + + def test_create_area_def_dynamic_areas(self): + """Test certain parameter combinations produce a DynamicAreaDefinition.""" + from pyresample import create_area_def as cad + from pyresample.geometry import DynamicAreaDefinition + projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m' + shape = (425, 850) + area_extent = (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625) + assert isinstance(cad('ease_sh', projection, shape=shape), DynamicAreaDefinition) + assert isinstance(cad('ease_sh', projection, area_extent=area_extent), DynamicAreaDefinition) + + def test_create_area_def_dynamic_omerc(self): + """Test 'omerc' projections work in 'create_area_def'.""" + from pyresample import create_area_def as cad + from pyresample.geometry import DynamicAreaDefinition area_def = cad('omerc_bb', {'ellps': 'WGS84', 'proj': 'omerc'}) - self.assertTrue(isinstance(area_def, DynamicAreaDefinition)) + assert isinstance(area_def, DynamicAreaDefinition) def _check_final_area_lon_lat_with_chunks(final_area, lons, lats, chunks): @@ -2331,12 +2389,12 @@ def test_get_geostationary_bbox(self): """Get the geostationary bbox.""" geos_area = MagicMock() lon_0 = 0 - geos_area.proj_dict = {'a': 6378169.00, - 'b': 6356583.80, - 'h': 35785831.00, - 'lon_0': lon_0, - 'proj': 'geos'} - geos_area._crs_or_proj_dict = geos_area.proj_dict + proj_dict = {'a': 6378169.00, + 'b': 6356583.80, + 'h': 35785831.00, + 'lon_0': lon_0, + 'proj': 'geos'} + geos_area.crs = CRS(proj_dict) geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.] lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20) @@ -2356,15 +2414,14 @@ def test_get_geostationary_bbox(self): np.testing.assert_allclose(lat, elat) geos_area = MagicMock() - del geos_area.crs_wkt lon_0 = 10 - geos_area.proj_dict = {'a': 6378169.00, - 'b': 6356583.80, - 'h': 35785831.00, - 'lon_0': lon_0, - 'proj': 'geos'} + proj_dict = {'a': 6378169.00, + 'b': 6356583.80, + 'h': 35785831.00, + 'lon_0': lon_0, + 'proj': 'geos'} + geos_area.crs = CRS(proj_dict) geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.] - geos_area._crs_or_proj_dict = geos_area.proj_dict lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20) np.testing.assert_allclose(lon, elon + lon_0) @@ -2372,7 +2429,7 @@ def test_get_geostationary_bbox(self): def test_get_geostationary_angle_extent(self): """Get max geostationary angles.""" geos_area = MagicMock() - geos_area.proj_dict = { + proj_dict = { 'proj': 'geos', 'sweep': 'x', 'lon_0': -89.5, @@ -2380,28 +2437,29 @@ def test_get_geostationary_angle_extent(self): 'b': 6356583.80, 'h': 35785831.00, 'units': 'm'} - geos_area._crs_or_proj_dict = geos_area.proj_dict + geos_area.crs = CRS(proj_dict) expected = (0.15185342867090912, 0.15133555510297725) np.testing.assert_allclose(expected, geometry.get_geostationary_angle_extent(geos_area)) - geos_area.proj_dict['a'] = 1000.0 - geos_area.proj_dict['b'] = 1000.0 - geos_area.proj_dict['h'] = np.sqrt(2) * 1000.0 - 1000.0 + proj_dict['a'] = 1000.0 + proj_dict['b'] = 1000.0 + proj_dict['h'] = np.sqrt(2) * 1000.0 - 1000.0 + geos_area.crs = CRS(proj_dict) expected = (np.deg2rad(45), np.deg2rad(45)) np.testing.assert_allclose(expected, geometry.get_geostationary_angle_extent(geos_area)) - geos_area.proj_dict = { + proj_dict = { 'proj': 'geos', 'sweep': 'x', 'lon_0': -89.5, 'ellps': 'GRS80', 'h': 35785831.00, 'units': 'm'} - geos_area._crs_or_proj_dict = geos_area.proj_dict + geos_area.crs = CRS(proj_dict) expected = (0.15185277703584374, 0.15133971368991794) np.testing.assert_allclose(expected, geometry.get_geostationary_angle_extent(geos_area)) diff --git a/pyresample/test/test_grid.py b/pyresample/test/test_grid.py index 4c75e4409..01ea8ab93 100644 --- a/pyresample/test/test_grid.py +++ b/pyresample/test/test_grid.py @@ -201,12 +201,9 @@ def test_single_lonlat(self): def test_proj4_string(self): """Test 'proj_str' property of AreaDefinition.""" - from pyresample.utils import is_pyproj2 proj4_string = self.area_def.proj_str - expected_string = '+a=6378144.0 +b=6356759.0 +lat_ts=50.0 +lon_0=8.0 +proj=stere +lat_0=50.0' - if is_pyproj2(): - expected_string = '+a=6378144 +k=1 +lat_0=50 +lon_0=8 ' \ - '+no_defs +proj=stere +rf=298.253168108487 ' \ - '+type=crs +units=m +x_0=0 +y_0=0' + expected_string = '+a=6378144 +k=1 +lat_0=50 +lon_0=8 ' \ + '+no_defs +proj=stere +rf=298.253168108487 ' \ + '+type=crs +units=m +x_0=0 +y_0=0' self.assertEqual( frozenset(proj4_string.split()), frozenset(expected_string.split())) diff --git a/pyresample/test/test_utils.py b/pyresample/test/test_utils.py index b0b372fd1..01baa93cd 100644 --- a/pyresample/test/test_utils.py +++ b/pyresample/test/test_utils.py @@ -19,12 +19,12 @@ import os import unittest -from unittest import mock import io import pathlib from tempfile import NamedTemporaryFile import numpy as np +from pyproj import CRS import uuid from pyresample.test.utils import create_test_longitude, create_test_latitude @@ -44,18 +44,13 @@ class TestLegacyAreaParser(unittest.TestCase): def test_area_parser_legacy(self): """Test legacy area parser.""" from pyresample import parse_area_file - from pyresample.utils import is_pyproj2 ease_nh, ease_sh = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'ease_nh', 'ease_sh') - if is_pyproj2(): - # pyproj 2.0+ adds some extra parameters - projection = ("{'R': '6371228', 'lat_0': '90', 'lon_0': '0', " - "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " - "'units': 'm', 'x_0': '0', 'y_0': '0'}") - else: - projection = ("{'a': '6371228.0', 'lat_0': '90.0', " - "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}") + # pyproj 2.0+ adds some extra parameters + projection = ("{'R': '6371228', 'lat_0': '90', 'lon_0': '0', " + "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " + "'units': 'm', 'x_0': '0', 'y_0': '0'}") nh_str = """Area ID: ease_nh Description: Arctic EASE grid Projection ID: ease_nh @@ -66,13 +61,9 @@ def test_area_parser_legacy(self): self.assertEqual(ease_nh.__str__(), nh_str) self.assertIsInstance(ease_nh.proj_dict['lat_0'], (int, float)) - if is_pyproj2(): - projection = ("{'R': '6371228', 'lat_0': '-90', 'lon_0': '0', " - "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " - "'units': 'm', 'x_0': '0', 'y_0': '0'}") - else: - projection = ("{'a': '6371228.0', 'lat_0': '-90.0', " - "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}") + projection = ("{'R': '6371228', 'lat_0': '-90', 'lon_0': '0', " + "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " + "'units': 'm', 'x_0': '0', 'y_0': '0'}") sh_str = """Area ID: ease_sh Description: Antarctic EASE grid Projection ID: ease_sh @@ -85,16 +76,11 @@ def test_area_parser_legacy(self): def test_load_area(self): from pyresample import load_area - from pyresample.utils import is_pyproj2 ease_nh = load_area(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'ease_nh') - if is_pyproj2(): - # pyproj 2.0+ adds some extra parameters - projection = ("{'R': '6371228', 'lat_0': '90', 'lon_0': '0', " - "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " - "'units': 'm', 'x_0': '0', 'y_0': '0'}") - else: - projection = ("{'a': '6371228.0', 'lat_0': '90.0', " - "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}") + # pyproj 2.0+ adds some extra parameters + projection = ("{'R': '6371228', 'lat_0': '90', 'lon_0': '0', " + "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " + "'units': 'm', 'x_0': '0', 'y_0': '0'}") nh_str = """Area ID: ease_nh Description: Arctic EASE grid Projection ID: ease_nh @@ -131,15 +117,10 @@ def test_area_parser_yaml(self): 'test_latlong') ease_nh, ease_sh, test_m, test_deg, test_latlong = test_areas - from pyresample.utils import is_pyproj2 - if is_pyproj2(): - # pyproj 2.0+ adds some extra parameters - projection = ("{'R': '6371228', 'lat_0': '-90', 'lon_0': '0', " - "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " - "'units': 'm', 'x_0': '0', 'y_0': '0'}") - else: - projection = ("{'a': '6371228.0', 'lat_0': '-90.0', " - "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}") + # pyproj 2.0+ adds some extra parameters + projection = ("{'R': '6371228', 'lat_0': '-90', 'lon_0': '0', " + "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " + "'units': 'm', 'x_0': '0', 'y_0': '0'}") nh_str = """Area ID: ease_nh Description: Arctic EASE grid Projection: {} @@ -172,14 +153,10 @@ def test_area_parser_yaml(self): Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""".format(projection) self.assertEqual(test_deg.__str__(), deg_str) - if is_pyproj2(): - # pyproj 2.0+ adds some extra parameters - projection = ("{'ellps': 'WGS84', 'no_defs': 'None', " - "'pm': '-81.36', 'proj': 'longlat', " - "'type': 'crs'}") - else: - projection = ("{'ellps': 'WGS84', 'lat_0': '27.12', " - "'lon_0': '-81.36', 'proj': 'longlat'}") + # pyproj 2.0+ adds some extra parameters + projection = ("{'ellps': 'WGS84', 'no_defs': 'None', " + "'pm': '-81.36', 'proj': 'longlat', " + "'type': 'crs'}") latlong_str = """Area ID: test_latlong Description: Basic latlong grid Projection: {} @@ -207,7 +184,7 @@ def test_dynamic_area_parser_yaml(self): self.assertIsInstance(test_area, DynamicAreaDefinition) self.assertTrue(hasattr(test_area, 'resolution')) - self.assertEqual(test_area.resolution, (1.0, 1.0)) + np.testing.assert_allclose(test_area.resolution, (1.0, 1.0)) def test_multiple_file_content(self): from pyresample import parse_area_file @@ -353,14 +330,6 @@ def test_proj4_radius_parameters_provided(self): np.testing.assert_almost_equal(a, 6378273) np.testing.assert_almost_equal(b, 6356889.44891) - # test again but force pyproj <2 behavior - with mock.patch.object(utils.proj4, 'CRS', None): - a, b = utils.proj4.proj4_radius_parameters( - '+proj=stere +a=6378273 +b=6356889.44891', - ) - np.testing.assert_almost_equal(a, 6378273) - np.testing.assert_almost_equal(b, 6356889.44891) - def test_proj4_radius_parameters_ellps(self): """Test proj4_radius_parameters with ellps.""" from pyresample import utils @@ -370,14 +339,6 @@ def test_proj4_radius_parameters_ellps(self): np.testing.assert_almost_equal(a, 6378137.) np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) - # test again but force pyproj <2 behavior - with mock.patch.object(utils.proj4, 'CRS', None): - a, b = utils.proj4.proj4_radius_parameters( - '+proj=stere +ellps=WGS84', - ) - np.testing.assert_almost_equal(a, 6378137.) - np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) - def test_proj4_radius_parameters_default(self): """Test proj4_radius_parameters with default parameters.""" from pyresample import utils @@ -388,15 +349,6 @@ def test_proj4_radius_parameters_default(self): np.testing.assert_almost_equal(a, 6378137.) np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) - # test again but force pyproj <2 behavior - with mock.patch.object(utils.proj4, 'CRS', None): - a, b = utils.proj4.proj4_radius_parameters( - '+proj=lcc +lat_0=10 +lat_1=10', - ) - # WGS84 - np.testing.assert_almost_equal(a, 6378137.) - np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) - def test_proj4_radius_parameters_spherical(self): """Test proj4_radius_parameters in case of a spherical earth.""" from pyresample import utils @@ -406,14 +358,6 @@ def test_proj4_radius_parameters_spherical(self): np.testing.assert_almost_equal(a, 6378273.) np.testing.assert_almost_equal(b, 6378273.) - # test again but force pyproj <2 behavior - with mock.patch.object(utils.proj4, 'CRS', None): - a, b = utils.proj4.proj4_radius_parameters( - '+proj=stere +R=6378273', - ) - np.testing.assert_almost_equal(a, 6378273.) - np.testing.assert_almost_equal(b, 6378273.) - def test_convert_proj_floats(self): from collections import OrderedDict import pyresample.utils as utils @@ -431,38 +375,34 @@ def test_convert_proj_floats(self): def test_proj4_str_dict_conversion(self): from pyresample import utils - proj_str = "+proj=lcc +ellps=WGS84 +lon_0=-95 +no_defs" + proj_str = "+proj=lcc +ellps=WGS84 +lon_0=-95 +lat_1=25.5 +no_defs" proj_dict = utils.proj4.proj4_str_to_dict(proj_str) proj_str2 = utils.proj4.proj4_dict_to_str(proj_dict) proj_dict2 = utils.proj4.proj4_str_to_dict(proj_str2) self.assertDictEqual(proj_dict, proj_dict2) - self.assertIsInstance(proj_dict['lon_0'], float) - self.assertIsInstance(proj_dict2['lon_0'], float) + self.assertIsInstance(proj_dict['lon_0'], (float, int)) + self.assertIsInstance(proj_dict2['lon_0'], (float, int)) + self.assertIsInstance(proj_dict['lat_1'], float) + self.assertIsInstance(proj_dict2['lat_1'], float) # EPSG proj_str = '+init=EPSG:4326' - proj_dict_exp = {'init': 'EPSG:4326'} proj_dict = utils.proj4.proj4_str_to_dict(proj_str) - self.assertEqual(proj_dict, proj_dict_exp) - self.assertEqual(utils.proj4.proj4_dict_to_str(proj_dict), proj_str) # round-trip + proj_str2 = utils.proj4.proj4_dict_to_str(proj_dict) + proj_dict2 = utils.proj4.proj4_str_to_dict(proj_str2) + # pyproj usually expands EPSG definitions so we can only round trip + self.assertEqual(proj_dict, proj_dict2) proj_str = 'EPSG:4326' - proj_dict_exp = {'init': 'EPSG:4326'} proj_dict_exp2 = {'proj': 'longlat', 'datum': 'WGS84', 'no_defs': None, 'type': 'crs'} proj_dict = utils.proj4.proj4_str_to_dict(proj_str) - if 'init' in proj_dict: - # pyproj <2.0 - self.assertEqual(proj_dict, proj_dict_exp) - else: - # pyproj 2.0+ - self.assertEqual(proj_dict, proj_dict_exp2) + self.assertEqual(proj_dict, proj_dict_exp2) # input != output for this style of EPSG code # EPSG to PROJ.4 can be lossy # self.assertEqual(utils._proj4.proj4_dict_to_str(proj_dict), proj_str) # round-trip def test_def2yaml_converter(self): from pyresample import parse_area_file, convert_def_to_yaml - from pyresample.utils import is_pyproj2 import tempfile def_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg') filehandle, yaml_file = tempfile.mkstemp() @@ -471,14 +411,6 @@ def test_def2yaml_converter(self): convert_def_to_yaml(def_file, yaml_file) areas_new = set(parse_area_file(yaml_file)) areas = parse_area_file(def_file) - for area in areas: - if is_pyproj2(): - # pyproj 2.0 adds units back in - # pyproj <2 doesn't - continue - # initialize _proj_dict - area.proj_dict # noqa - area._proj_dict.pop('units', None) areas_old = set(areas) areas_new = {area.area_id: area for area in areas_new} areas_old = {area.area_id: area for area in areas_old} @@ -488,19 +420,13 @@ def test_def2yaml_converter(self): def test_get_area_def_from_raster(self): from pyresample import utils - from rasterio.crs import CRS + from rasterio.crs import CRS as RCRS from affine import Affine x_size = 791 y_size = 718 transform = Affine(300.0379266750948, 0.0, 101985.0, 0.0, -300.041782729805, 2826915.0) - crs = CRS(init='epsg:3857') - if utils.is_pyproj2(): - # pyproj 2.0+ expands CRS parameters - from pyproj import CRS - proj_dict = CRS(3857).to_dict() - else: - proj_dict = crs.to_dict() + crs = RCRS(init='epsg:3857') source = tmptiff(x_size, y_size, transform, crs=crs) area_id = 'area_id' proj_id = 'proj_id' @@ -512,14 +438,14 @@ def test_get_area_def_from_raster(self): self.assertEqual(area_def.description, description) self.assertEqual(area_def.width, x_size) self.assertEqual(area_def.height, y_size) - self.assertDictEqual(proj_dict, area_def.proj_dict) + self.assertEqual(crs, area_def.crs) self.assertTupleEqual(area_def.area_extent, (transform.c, transform.f + transform.e * y_size, transform.c + transform.a * x_size, transform.f)) def test_get_area_def_from_raster_extracts_proj_id(self): - from rasterio.crs import CRS + from rasterio.crs import CRS as RCRS from pyresample import utils - crs = CRS(init='epsg:3857') + crs = RCRS(init='epsg:3857') source = tmptiff(crs=crs) area_def = utils.rasterio.get_area_def_from_raster(source) epsg3857_names = ( @@ -552,10 +478,7 @@ def test_get_area_def_from_raster_non_georef_respects_proj_dict(self): source = tmptiff(transform=transform) proj_dict = {'init': 'epsg:3857'} area_def = utils.rasterio.get_area_def_from_raster(source, proj_dict=proj_dict) - if utils.is_pyproj2(): - from pyproj import CRS - proj_dict = CRS(3857).to_dict() - self.assertDictEqual(area_def.proj_dict, proj_dict) + self.assertEqual(area_def.crs, CRS(3857)) class TestProjRotation(unittest.TestCase): @@ -729,7 +652,7 @@ def test_load_cf_from_wrong_filepath(self): # wrong case #2: the path exists, but is not a netCDF file cf_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml') - self.assertRaises(OSError, load_cf_area, cf_file) + self.assertRaises((ValueError, OSError), load_cf_area, cf_file) def test_load_cf_parameters_errors(self): from pyresample.utils import load_cf_area diff --git a/pyresample/utils/__init__.py b/pyresample/utils/__init__.py index aff43e277..98eb76e80 100644 --- a/pyresample/utils/__init__.py +++ b/pyresample/utils/__init__.py @@ -24,7 +24,8 @@ import pyproj import warnings -from .proj4 import (proj4_dict_to_str, proj4_str_to_dict, convert_proj_floats, proj4_radius_parameters) # noqa +from .proj4 import (proj4_dict_to_str, proj4_str_to_dict, convert_proj_floats, # noqa + proj4_radius_parameters, get_geostationary_height) # noqa from .rasterio import get_area_def_from_raster # noqa from .cf import load_cf_area # noqa @@ -241,11 +242,6 @@ def recursive_dict_update(d, u): return d -def is_pyproj2(): - """Determine whether the current pyproj version is >= 2.0""" - return pyproj.__version__ >= '2' - - def check_slice_orientation(sli): """Check that the slice is slicing the right way.""" if sli.start > sli.stop: diff --git a/pyresample/utils/cf.py b/pyresample/utils/cf.py index 9ea709fea..4e09e4554 100644 --- a/pyresample/utils/cf.py +++ b/pyresample/utils/cf.py @@ -315,11 +315,9 @@ def _load_cf_area_one_variable_areadef(axis_info, crs, unit, grid_mapping_variab # get area extent from the x and y info extent = _get_area_extent_from_cf_axis(axis_info['x'], axis_info['y']) - # transform the crs objecto a proj_dict (might not be needed in future versions of pyresample) - proj_dict = crs.to_dict() - # finally prepare the AreaDefinition object - return geometry.AreaDefinition.from_extent(grid_mapping_variable, proj_dict, shape, extent, units=unit) + return geometry.AreaDefinition.from_extent(grid_mapping_variable, crs, + shape, extent, units=unit) def _load_cf_area_one_variable(nc_handle, variable, y=None, x=None): diff --git a/pyresample/utils/proj4.py b/pyresample/utils/proj4.py index 8e8be8456..09b52da9b 100644 --- a/pyresample/utils/proj4.py +++ b/pyresample/utils/proj4.py @@ -16,12 +16,10 @@ # You should have received a copy of the GNU Lesser General Public License along # with this program. If not, see . +import math from collections import OrderedDict -try: - from pyproj.crs import CRS -except ImportError: - CRS = None +from pyproj import CRS def convert_proj_floats(proj_pairs): @@ -51,18 +49,8 @@ def proj4_str_to_dict(proj4_str): Note: Key only parameters will be assigned a value of `True`. """ # convert EPSG codes to equivalent PROJ4 string definition - if proj4_str.startswith('EPSG:') and CRS is not None: - crs = CRS(proj4_str) - if hasattr(crs, 'to_dict'): - # pyproj 2.2+ - return crs.to_dict() - proj4_str = crs.to_proj4() - elif proj4_str.startswith('EPSG:'): - # legacy +init= PROJ4 string and no pyproj 2.0+ to help convert - proj4_str = "+init={}".format(proj4_str) - - pairs = (x.split('=', 1) for x in proj4_str.replace('+', '').split(" ")) - return convert_proj_floats(pairs) + crs = CRS(proj4_str) + return crs.to_dict() def proj4_dict_to_str(proj4_dict, sort=False): @@ -90,44 +78,16 @@ def proj4_radius_parameters(proj4_dict): Returns: a (float), b (float): equatorial and polar radius """ - if CRS is not None: - import math - crs = CRS(proj4_dict) - a = crs.ellipsoid.semi_major_metre - b = crs.ellipsoid.semi_minor_metre - if not math.isnan(b): - return a, b - # older versions of pyproj didn't always have a valid minor radius - proj4_dict = crs.to_dict() - - if isinstance(proj4_dict, str): - new_info = proj4_str_to_dict(proj4_dict) - else: - new_info = proj4_dict.copy() - - # load information from PROJ.4 about the ellipsis if possible - - from pyproj import Geod - - if 'ellps' in new_info: - geod = Geod(**new_info) - new_info['a'] = geod.a - new_info['b'] = geod.b - elif 'a' not in new_info or 'b' not in new_info: - - if 'rf' in new_info and 'f' not in new_info: - new_info['f'] = 1. / float(new_info['rf']) - - if 'a' in new_info and 'f' in new_info: - new_info['b'] = float(new_info['a']) * (1 - float(new_info['f'])) - elif 'b' in new_info and 'f' in new_info: - new_info['a'] = float(new_info['b']) / (1 - float(new_info['f'])) - elif 'R' in new_info: - new_info['a'] = new_info['R'] - new_info['b'] = new_info['R'] - else: - geod = Geod(**{'ellps': 'WGS84'}) - new_info['a'] = geod.a - new_info['b'] = geod.b - - return float(new_info['a']), float(new_info['b']) + crs = proj4_dict + if not isinstance(crs, CRS): + crs = CRS(crs) + a = crs.ellipsoid.semi_major_metre + b = crs.ellipsoid.semi_minor_metre + if not math.isnan(b): + return a, b + + +def get_geostationary_height(geos_area_crs): + params = geos_area_crs.coordinate_operation.params + h_param = [p for p in params if 'satellite height' in p.name.lower()][0] + return h_param.value diff --git a/setup.py b/setup.py index aa477d30c..f5aca91be 100644 --- a/setup.py +++ b/setup.py @@ -26,7 +26,7 @@ from setuptools import Extension, find_packages, setup from Cython.Build import cythonize -requirements = ['setuptools>=3.2', 'pyproj>=1.9.5.1', 'configobj', +requirements = ['setuptools>=3.2', 'pyproj>=2.2', 'configobj', 'pykdtree>=1.3.1', 'pyyaml', 'numpy>=1.10.0', ] extras_require = {'numexpr': ['numexpr'],