From 89a53783d35402aba5fecd2fc9e2b61feeb412c9 Mon Sep 17 00:00:00 2001 From: Jon Thielen Date: Wed, 26 Aug 2020 17:41:54 -0500 Subject: [PATCH] Use PyProj instead of Cartopy for internal coordinate transforms and rename crs coordinate --- src/metpy/calc/cross_sections.py | 15 +++++---- src/metpy/calc/tools.py | 18 +++++----- src/metpy/interpolate/slices.py | 15 +++++---- src/metpy/plots/mapping.py | 6 ++++ src/metpy/xarray.py | 56 +++++++++++++++++++++++-------- tests/calc/test_calc_tools.py | 44 +++++++++++++----------- tests/calc/test_cross_sections.py | 6 ++-- tests/calc/test_kinematics.py | 6 ++-- tests/interpolate/test_slices.py | 12 +++---- tests/test_xarray.py | 32 +++++++++--------- tutorials/xarray_tutorial.py | 20 +++++++---- 11 files changed, 139 insertions(+), 91 deletions(-) diff --git a/src/metpy/calc/cross_sections.py b/src/metpy/calc/cross_sections.py index bbae1520762..590f263edf5 100644 --- a/src/metpy/calc/cross_sections.py +++ b/src/metpy/calc/cross_sections.py @@ -36,9 +36,7 @@ def distances_from_cross_section(cross): """ if check_axis(cross.metpy.x, 'longitude') and check_axis(cross.metpy.y, 'latitude'): # Use pyproj to obtain x and y distances - from pyproj import Geod - - g = Geod(cross.metpy.cartopy_crs.proj4_init) + g = cross.metpy.pyproj_crs.get_geod() lon = cross.metpy.x lat = cross.metpy.y @@ -83,10 +81,13 @@ def latitude_from_cross_section(cross): if check_axis(y, 'latitude'): return y else: - import cartopy.crs as ccrs - latitude = ccrs.Geodetic().transform_points(cross.metpy.cartopy_crs, - cross.metpy.x.values, - y.values)[..., 1] + from pyproj import Proj + latitude = Proj(cross.metpy.pyproj_crs)( + cross.metpy.x.values, + y.values, + inverse=True, + radians=False + )[1] latitude = xr.DataArray(latitude * units.degrees_north, coords=y.coords, dims=y.dims) return latitude diff --git a/src/metpy/calc/tools.py b/src/metpy/calc/tools.py index eadaa0b2dea..2a485c17c11 100644 --- a/src/metpy/calc/tools.py +++ b/src/metpy/calc/tools.py @@ -751,7 +751,7 @@ def take(indexer): @exporter.export @preprocess_xarray -def lat_lon_grid_deltas(longitude, latitude, y_dim=-2, x_dim=-1, **kwargs): +def lat_lon_grid_deltas(longitude, latitude, y_dim=-2, x_dim=-1, geod=None): r"""Calculate the actual delta between grid points that are in latitude/longitude format. Parameters @@ -766,8 +766,9 @@ def lat_lon_grid_deltas(longitude, latitude, y_dim=-2, x_dim=-1, **kwargs): axis number for the y dimesion, defaults to -2. x_dim: int axis number for the x dimension, defaults to -1. - kwargs - Other keyword arguments to pass to :class:`~pyproj.Geod` + geod : pyproj.Geod or None + PyProj Geod to use for forward azimuth and distance calculations. If None, use a + default spherical ellipsoid. Returns ------- @@ -804,11 +805,10 @@ def lat_lon_grid_deltas(longitude, latitude, y_dim=-2, x_dim=-1, **kwargs): take_y = make_take(latitude.ndim, y_dim) take_x = make_take(latitude.ndim, x_dim) - geod_args = {'ellps': 'sphere'} - if kwargs: - geod_args = kwargs - - g = Geod(**geod_args) + if geod is None: + g = Geod(ellps='sphere') + else: + g = geod forward_az, _, dy = g.inv(longitude[take_y(slice(None, -1))], latitude[take_y(slice(None, -1))], @@ -922,7 +922,7 @@ def grid_deltas_from_dataarray(f, kind='default'): (dx_var, dx_units), (dy_var, dy_units) = ( (xr.Variable(dims=latitude.dims, data=deltas.magnitude), deltas.units) for deltas in lat_lon_grid_deltas(longitude, latitude, y_dim=y_dim, x_dim=x_dim, - initstring=f.metpy.cartopy_crs.proj4_init)) + geod=f.metpy.pyproj_crs.get_geod())) else: # Obtain y/x coordinate differences y, x = f.metpy.coordinates('y', 'x') diff --git a/src/metpy/interpolate/slices.py b/src/metpy/interpolate/slices.py index 053644fb2ce..ca204b67530 100644 --- a/src/metpy/interpolate/slices.py +++ b/src/metpy/interpolate/slices.py @@ -74,8 +74,8 @@ def geodesic(crs, start, end, steps): Parameters ---------- - crs: `cartopy.crs` - Cartopy Coordinate Reference System to use for the output + crs: `pyproj.CRS` + PyProj Coordinate Reference System to use for the output start: (2, ) array_like A latitude-longitude pair designating the start point of the geodesic (units are degrees north and degrees east). @@ -96,18 +96,19 @@ def geodesic(crs, start, end, steps): cross_section """ - import cartopy.crs as ccrs - from pyproj import Geod + from pyproj import Proj + + g = crs.get_geod() + p = Proj(crs) # Geod.npts only gives points *in between* the start and end, and we want to include # the endpoints. - g = Geod(crs.proj4_init) geodesic = np.concatenate([ np.array(start[::-1])[None], np.array(g.npts(start[1], start[0], end[1], end[0], steps - 2)), np.array(end[::-1])[None] ]).transpose() - points = crs.transform_points(ccrs.Geodetic(), *geodesic)[:, :2] + points = np.stack(p(geodesic[0], geodesic[1], inverse=False, radians=False), axis=-1) return points @@ -162,7 +163,7 @@ def cross_section(data, start, end, steps=100, interp_type='linear'): # Get the projection and coordinates try: - crs_data = data.metpy.cartopy_crs + crs_data = data.metpy.pyproj_crs x = data.metpy.x except AttributeError: raise ValueError('Data missing required coordinate information. Verify that ' diff --git a/src/metpy/plots/mapping.py b/src/metpy/plots/mapping.py index 6d0d16ad661..1b13eca14a5 100644 --- a/src/metpy/plots/mapping.py +++ b/src/metpy/plots/mapping.py @@ -77,6 +77,12 @@ def to_cartopy(self): return proj_handler(self._attrs, globe) + def to_pyproj(self): + """Convert to a PyProj CRS.""" + import pyproj + + return pyproj.CRS.from_cf(self._attrs) + def to_dict(self): """Get the dictionary of metadata attributes.""" return self._attrs.copy() diff --git a/src/metpy/xarray.py b/src/metpy/xarray.py index 030ffb64a0b..bf88d30d41c 100644 --- a/src/metpy/xarray.py +++ b/src/metpy/xarray.py @@ -212,8 +212,8 @@ def dequantify(self): @property def crs(self): """Return the coordinate reference system (CRS) as a CFProjection object.""" - if 'crs' in self._data_array.coords: - return self._data_array.coords['crs'].item() + if 'metpy_crs' in self._data_array.coords: + return self._data_array.coords['metpy_crs'].item() raise AttributeError('crs attribute is not available.') @property @@ -231,6 +231,11 @@ def cartopy_geodetic(self): """Return the Geodetic CRS associated with the native CRS globe.""" return self.crs.cartopy_geodetic + @property + def pyproj_crs(self): + """Return the coordinate reference system (CRS) as a pyproj object.""" + return self.crs.to_pyproj() + def _fixup_coordinate_map(self, coord_map): """Ensure sure we have coordinate variables in map, not coordinate names.""" new_coord_map = {} @@ -569,7 +574,7 @@ def assign_latitude_longitude(self, force=False): Notes ----- - A valid CRS coordinate must be present. Cartopy is used for the coordinate + A valid CRS coordinate must be present. PyProj is used for the coordinate transformations. """ @@ -605,7 +610,7 @@ def assign_y_x(self, force=False, tolerance=None): Notes ----- - A valid CRS coordinate must be present. Cartopy is used for the coordinate + A valid CRS coordinate must be present. PyProj is used for the coordinate transformations. """ @@ -632,7 +637,7 @@ class MetPyDatasetAccessor: >>> import xarray as xr >>> from metpy.cbook import get_test_data >>> ds = xr.open_dataset(get_test_data('narr_example.nc', False)).metpy.parse_cf() - >>> print(ds['crs'].item()) + >>> print(ds['metpy_crs'].item()) Projection: lambert_conformal_conic """ @@ -675,6 +680,22 @@ def parse_cf(self, varname=None, coordinates=None): var = self._dataset[varname] + # Check for crs conflict + if varname == 'metpy_crs': + warnings.warn( + 'Attempting to parse metpy_crs as a data variable. Unexpected merge conflicts ' + 'may occur.' + ) + elif 'metpy_crs' in var.coords: + try: + assert isinstance(var.coords['metpy_crs'].item(), CFProjection) + except (ValueError, AssertionError): + # Catch non-scalar and non-CFProjection coordinates + warnings.warn( + 'metpy_crs already present as a non-CFProjection coordinate. Unexpected ' + 'merge conflicts may occur.' + ) + # Assign coordinates if the coordinates argument is given if coordinates is not None: var = var.metpy.assign_coordinates(coordinates) @@ -709,7 +730,7 @@ def _has_coord(coord_type): # Rebuild the coordinates of the dataarray, and return quantified DataArray var = self._rebuild_coords(var, crs) if crs is not None: - var = var.assign_coords(coords={'crs': crs}) + var = var.assign_coords(coords={'metpy_crs': crs}) return var def _rebuild_coords(self, var, crs): @@ -795,7 +816,7 @@ def assign_latitude_longitude(self, force=False): Notes ----- - A valid CRS coordinate must be present. Cartopy is used for the coordinate + A valid CRS coordinate must be present. PyProj is used for the coordinate transformations. """ @@ -840,7 +861,7 @@ def assign_y_x(self, force=False, tolerance=None): Notes ----- - A valid CRS coordinate must be present. Cartopy is used for the coordinate + A valid CRS coordinate must be present. PyProj is used for the coordinate transformations. """ @@ -994,14 +1015,16 @@ def _assign_crs(xarray_object, cf_attributes, cf_kwargs): attrs = cf_attributes if cf_attributes is not None else cf_kwargs # Assign crs coordinate to xarray object - return xarray_object.assign_coords(crs=CFProjection(attrs)) + return xarray_object.assign_coords(metpy_crs=CFProjection(attrs)) def _build_latitude_longitude(da): """Build latitude/longitude coordinates from DataArray's y/x coordinates.""" + from pyproj import Proj + y, x = da.metpy.coordinates('y', 'x') xx, yy = np.meshgrid(x.values, y.values) - lonlats = da.metpy.cartopy_geodetic.transform_points(da.metpy.cartopy_crs, xx, yy) + lonlats = np.stack(Proj(da.metpy.pyproj_crs)(xx, yy, inverse=True, radians=False), axis=-1) longitude = xr.DataArray(lonlats[..., 0], dims=(y.name, x.name), coords={y.name: y, x.name: x}, attrs={'units': 'degrees_east', 'standard_name': 'longitude'}) @@ -1013,6 +1036,8 @@ def _build_latitude_longitude(da): def _build_y_x(da, tolerance): """Build y/x coordinates from DataArray's latitude/longitude coordinates.""" + from pyproj import Proj + # Initial sanity checks latitude, longitude = da.metpy.coordinates('latitude', 'longitude') if latitude.dims != longitude.dims: @@ -1022,9 +1047,12 @@ def _build_y_x(da, tolerance): 'must be 2D') # Convert to projected y/x - xxyy = da.metpy.cartopy_crs.transform_points(da.metpy.cartopy_geodetic, - longitude.values, - latitude.values) + xxyy = np.stack(Proj(da.metpy.pyproj_crs)( + longitude.values, + latitude.values, + inverse=False, + radians=False + ), axis=-1) # Handle tolerance tolerance = 1 if tolerance is None else tolerance.m_as('m') @@ -1051,7 +1079,7 @@ def _build_y_x(da, tolerance): else: raise ValueError('Projected y and x coordinates cannot be collapsed to 1D within ' 'tolerance. Verify that your latitude and longitude coordinates ' - 'correpsond to your CRS coordinate.') + 'correspond to your CRS coordinate.') def preprocess_xarray(func): diff --git a/tests/calc/test_calc_tools.py b/tests/calc/test_calc_tools.py index 50f298b8fd5..e53a799d8e4 100644 --- a/tests/calc/test_calc_tools.py +++ b/tests/calc/test_calc_tools.py @@ -497,9 +497,11 @@ def test_lat_lon_grid_deltas_mismatched_shape(): @needs_pyproj def test_lat_lon_grid_deltas_geod_kwargs(): """Test that geod kwargs are overridden by users #774.""" + from pyproj import Geod + lat = np.arange(40, 50, 2.5) lon = np.arange(-100, -90, 2.5) - dx, dy = lat_lon_grid_deltas(lon, lat, a=4370997) + dx, dy = lat_lon_grid_deltas(lon, lat, geod=Geod(a=4370997)) dx_truth = np.array([[146095.76101984, 146095.76101984, 146095.76101984], [140608.9751528, 140608.9751528, 140608.9751528], [134854.56713287, 134854.56713287, 134854.56713287], @@ -1095,29 +1097,33 @@ def test_grid_deltas_from_dataarray_xy(test_da_xy): assert_array_almost_equal(dy, true_dy, 5) -def test_grid_deltas_from_dataarray_actual_xy(test_da_xy, ccrs): +@needs_pyproj +def test_grid_deltas_from_dataarray_actual_xy(test_da_xy): """Test grid_deltas_from_dataarray with a xy grid and kind='actual'.""" # Construct lon/lat coordinates + from pyproj import Proj y, x = xr.broadcast(*test_da_xy.metpy.coordinates('y', 'x')) - lonlat = (ccrs.Geodetic(test_da_xy.metpy.cartopy_globe) - .transform_points(test_da_xy.metpy.cartopy_crs, x.values, y.values)) - lon = lonlat[..., 0] - lat = lonlat[..., 1] + lon, lat = Proj(test_da_xy.metpy.pyproj_crs)( + x.values, + y.values, + inverse=True, + radians=False + ) test_da_xy = test_da_xy.assign_coords( longitude=xr.DataArray(lon, dims=('y', 'x'), attrs={'units': 'degrees_east'}), latitude=xr.DataArray(lat, dims=('y', 'x'), attrs={'units': 'degrees_north'})) # Actually test calculation dx, dy = grid_deltas_from_dataarray(test_da_xy, kind='actual') - true_dx = [[[[494426.3249766, 493977.6028005, 493044.0656467], - [498740.2046073, 498474.9771064, 497891.6588559], - [500276.2649627, 500256.3440237, 500139.9484845], - [498740.6956936, 499045.0391707, 499542.7244501]]]] * units.m - true_dy = [[[[496862.4106337, 496685.4729999, 496132.0732114, 495137.8882404], - [499774.9676486, 499706.3354977, 499467.5546773, 498965.2587818], - [499750.8962991, 499826.2263137, 500004.4977747, 500150.9897759]]]] * units.m - assert_array_almost_equal(dx, true_dx, 3) - assert_array_almost_equal(dy, true_dy, 3) + true_dx = [[[[494152.626, 493704.152, 492771.132], + [498464.118, 498199.037, 497616.042], + [499999.328, 499979.418, 499863.087], + [498464.608, 498768.783, 499266.193]]]] * units.m + true_dy = [[[[496587.363, 496410.523, 495857.430, 494863.795], + [499498.308, 499429.714, 499191.065, 498689.047], + [499474.250, 499549.538, 499727.711, 499874.122]]]] * units.m + assert_array_almost_equal(dx, true_dx, 2) + assert_array_almost_equal(dy, true_dy, 2) def test_grid_deltas_from_dataarray_nominal_lonlat(test_da_lonlat): @@ -1177,7 +1183,7 @@ def test_first_derivative_xarray_lonlat(test_da_lonlat): coords=(('lat', test_da_lonlat['lat']),) ) _, truth = xr.broadcast(test_da_lonlat, partial) - truth.coords['crs'] = test_da_lonlat['crs'] + truth.coords['metpy_crs'] = test_da_lonlat['metpy_crs'] truth.attrs['units'] = 'kelvin / meter' truth = truth.metpy.quantify() @@ -1231,7 +1237,7 @@ def test_second_derivative_xarray_lonlat(test_da_lonlat): coords=(('lat', test_da_lonlat['lat']),) ) _, truth = xr.broadcast(test_da_lonlat, partial) - truth.coords['crs'] = test_da_lonlat['crs'] + truth.coords['metpy_crs'] = test_da_lonlat['metpy_crs'] truth.attrs['units'] = 'kelvin / meter^2' truth = truth.metpy.quantify() @@ -1259,7 +1265,7 @@ def test_gradient_xarray(test_da_xy): coords=(('isobaric', test_da_xy['isobaric']),) ) _, truth_p = xr.broadcast(test_da_xy, partial) - truth_p.coords['crs'] = test_da_xy['crs'] + truth_p.coords['metpy_crs'] = test_da_xy['metpy_crs'] truth_p.attrs['units'] = 'kelvin / hectopascal' truth_p = truth_p.metpy.quantify() @@ -1343,7 +1349,7 @@ def test_laplacian_xarray_lonlat(test_da_lonlat): coords=(('lat', test_da_lonlat['lat']),) ) _, truth = xr.broadcast(test_da_lonlat, partial) - truth.coords['crs'] = test_da_lonlat['crs'] + truth.coords['metpy_crs'] = test_da_lonlat['metpy_crs'] truth.attrs['units'] = 'kelvin / meter^2' truth = truth.metpy.quantify() diff --git a/tests/calc/test_cross_sections.py b/tests/calc/test_cross_sections.py index 259a25d59d3..c2b09f477d9 100644 --- a/tests/calc/test_cross_sections.py +++ b/tests/calc/test_cross_sections.py @@ -119,7 +119,7 @@ def test_distances_from_cross_section_given_lonlat(test_cross_lonlat): true_x = xr.DataArray( true_x_values * units.meters, coords={ - 'crs': test_cross_lonlat['crs'], + 'metpy_crs': test_cross_lonlat['metpy_crs'], 'lat': test_cross_lonlat['lat'], 'lon': test_cross_lonlat['lon'], 'index': index, @@ -129,7 +129,7 @@ def test_distances_from_cross_section_given_lonlat(test_cross_lonlat): true_y = xr.DataArray( true_y_values * units.meters, coords={ - 'crs': test_cross_lonlat['crs'], + 'metpy_crs': test_cross_lonlat['metpy_crs'], 'lat': test_cross_lonlat['lat'], 'lon': test_cross_lonlat['lon'], 'index': index, @@ -168,7 +168,7 @@ def test_latitude_from_cross_section_given_y(test_cross_xy): true_latitude = xr.DataArray( true_latitude_values * units.degrees_north, coords={ - 'crs': test_cross_xy['crs'], + 'metpy_crs': test_cross_xy['metpy_crs'], 'y': test_cross_xy['y'], 'x': test_cross_xy['x'], 'index': index, diff --git a/tests/calc/test_kinematics.py b/tests/calc/test_kinematics.py index 0213dc645ec..ca661a04c54 100644 --- a/tests/calc/test_kinematics.py +++ b/tests/calc/test_kinematics.py @@ -18,7 +18,7 @@ total_deformation, vorticity, wind_components) from metpy.constants import g, omega, Re from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_array_equal, - get_test_data, needs_cartopy, needs_pyproj) + get_test_data, needs_pyproj) from metpy.units import concatenate, units @@ -1018,7 +1018,7 @@ def test_q_vector_with_static_stability(q_vector_data): @pytest.fixture -@needs_cartopy +@needs_pyproj def data_4d(): """Define 4D data (extracted from Irma GFS example) for testing kinematics functions.""" data = xr.open_dataset(get_test_data('irma_gfs_example.nc', False)) @@ -1035,7 +1035,7 @@ def data_4d(): ).isel(time1=[0, 1, 2]) dx, dy = lat_lon_grid_deltas(subset['longitude'].values, subset['latitude'].values, - initstring=subset['longitude'].metpy.cartopy_crs.proj4_init) + geod=subset['longitude'].metpy.pyproj_crs.get_geod()) return namedtuple('D_4D_Test_Data', 'height temperature pressure u v dx dy latitude')( subset['Geopotential_height_isobaric'].metpy.unit_array, diff --git a/tests/interpolate/test_slices.py b/tests/interpolate/test_slices.py index 5f51a9216a1..d992a928ff0 100644 --- a/tests/interpolate/test_slices.py +++ b/tests/interpolate/test_slices.py @@ -110,7 +110,7 @@ def test_interpolate_to_slice_against_selection(test_ds_lonlat): @needs_cartopy def test_geodesic(test_ds_xy): """Test the geodesic construction.""" - crs = test_ds_xy['temperature'].metpy.cartopy_crs + crs = test_ds_xy['temperature'].metpy.pyproj_crs path = geodesic(crs, (36.46, -112.45), (42.95, -68.74), 7) truth = np.array([[-4.99495719e+05, -1.49986599e+06], [9.84044354e+04, -1.26871737e+06], @@ -150,7 +150,7 @@ def test_cross_section_dataarray_and_linear_interp(test_ds_xy): truth_values_x, name='x', coords={ - 'crs': data['crs'], + 'metpy_crs': data['metpy_crs'], 'y': (['index'], truth_values_y), 'x': (['index'], truth_values_x), 'index': index, @@ -161,7 +161,7 @@ def test_cross_section_dataarray_and_linear_interp(test_ds_xy): truth_values_y, name='y', coords={ - 'crs': data['crs'], + 'metpy_crs': data['metpy_crs'], 'y': (['index'], truth_values_y), 'x': (['index'], truth_values_x), 'index': index, @@ -175,7 +175,7 @@ def test_cross_section_dataarray_and_linear_interp(test_ds_xy): 'time': data['time'], 'isobaric': data['isobaric'], 'index': index, - 'crs': data['crs'], + 'metpy_crs': data['metpy_crs'], 'y': data_truth_y, 'x': data_truth_x }, @@ -216,7 +216,7 @@ def test_cross_section_dataset_and_nearest_interp(test_ds_lonlat): coords={ 'isobaric': test_ds_lonlat['isobaric'], 'index': index, - 'crs': test_ds_lonlat['crs'], + 'metpy_crs': test_ds_lonlat['metpy_crs'], 'lat': (['index'], truth_values_lat), 'lon': (['index'], truth_values_lon) }, @@ -242,7 +242,7 @@ def test_cross_section_error_on_missing_coordinate(test_ds_lonlat): """Test that the proper error is raised with missing coordinate.""" # Use a variable with no crs coordinate data_bad = test_ds_lonlat['temperature'].copy() - del data_bad['crs'] + del data_bad['metpy_crs'] start, end = (30.5, 255.5), (44.5, 274.5) with pytest.raises(ValueError): diff --git a/tests/test_xarray.py b/tests/test_xarray.py index 98b3bd5811e..99676067459 100644 --- a/tests/test_xarray.py +++ b/tests/test_xarray.py @@ -10,7 +10,7 @@ from metpy.plots.mapping import CFProjection from metpy.testing import (assert_almost_equal, assert_array_almost_equal, assert_array_equal, - get_test_data, needs_cartopy) + get_test_data, needs_pyproj) from metpy.units import units from metpy.xarray import check_axis, check_matching_coordinates, preprocess_xarray @@ -212,7 +212,7 @@ def test_missing_grid_mapping(): ds = xr.Dataset({'data': data}) data_var = ds.metpy.parse_cf('data') - assert 'crs' in data_var.coords + assert 'metpy_crs' in data_var.coords def test_missing_grid_mapping_var(caplog): @@ -775,7 +775,7 @@ def test_assign_crs_dataarray_by_argument(test_ds_generic, ccrs): da = test_ds_generic['test'] new_da = da.metpy.assign_crs(sample_cf_attrs) assert isinstance(new_da.metpy.cartopy_crs, ccrs.LambertConformal) - assert new_da['crs'] == CFProjection(sample_cf_attrs) + assert new_da['metpy_crs'] == CFProjection(sample_cf_attrs) def test_assign_crs_dataarray_by_kwargs(test_ds_generic, ccrs): @@ -783,21 +783,21 @@ def test_assign_crs_dataarray_by_kwargs(test_ds_generic, ccrs): da = test_ds_generic['test'] new_da = da.metpy.assign_crs(**sample_cf_attrs) assert isinstance(new_da.metpy.cartopy_crs, ccrs.LambertConformal) - assert new_da['crs'] == CFProjection(sample_cf_attrs) + assert new_da['metpy_crs'] == CFProjection(sample_cf_attrs) def test_assign_crs_dataset_by_argument(test_ds_generic, ccrs): """Test assigning CRS to Dataset by projection dict.""" new_ds = test_ds_generic.metpy.assign_crs(sample_cf_attrs) assert isinstance(new_ds['test'].metpy.cartopy_crs, ccrs.LambertConformal) - assert new_ds['crs'] == CFProjection(sample_cf_attrs) + assert new_ds['metpy_crs'] == CFProjection(sample_cf_attrs) def test_assign_crs_dataset_by_kwargs(test_ds_generic, ccrs): """Test assigning CRS to Dataset by projection kwargs.""" new_ds = test_ds_generic.metpy.assign_crs(**sample_cf_attrs) assert isinstance(new_ds['test'].metpy.cartopy_crs, ccrs.LambertConformal) - assert new_ds['crs'] == CFProjection(sample_cf_attrs) + assert new_ds['metpy_crs'] == CFProjection(sample_cf_attrs) def test_assign_crs_error_with_both_attrs(test_ds_generic): @@ -835,7 +835,7 @@ def test_coord_helper_da_yx(): dims=('y', 'x'), coords={'y': np.linspace(0, 1e5, 3), 'x': np.linspace(-1e5, 0, 3), - 'crs': CFProjection(sample_cf_attrs)}) + 'metpy_crs': CFProjection(sample_cf_attrs)}) @pytest.fixture @@ -867,7 +867,7 @@ def test_coord_helper_da_latlon(): ), dims=('y', 'x') ), - 'crs': CFProjection(sample_cf_attrs) + 'metpy_crs': CFProjection(sample_cf_attrs) } ) @@ -878,7 +878,7 @@ def test_coord_helper_da_dummy_yx(test_coord_helper_da_latlon): return test_coord_helper_da_latlon.assign_coords(y=range(3), x=range(3)) -@needs_cartopy +@needs_pyproj def test_assign_latitude_longitude_basic_dataarray(test_coord_helper_da_yx, test_coord_helper_da_latlon): """Test assign_latitude_longitude in basic usage on DataArray.""" @@ -898,7 +898,7 @@ def test_assign_latitude_longitude_error_existing_dataarray( assert 'Latitude/longitude coordinate(s) are present' in str(exc) -@needs_cartopy +@needs_pyproj def test_assign_latitude_longitude_force_existing_dataarray( test_coord_helper_da_dummy_latlon, test_coord_helper_da_latlon): """Test assign_latitude_longitude with existing coordinates forcing new.""" @@ -910,7 +910,7 @@ def test_assign_latitude_longitude_force_existing_dataarray( lon.values, 3) -@needs_cartopy +@needs_pyproj def test_assign_latitude_longitude_basic_dataset(test_coord_helper_da_yx, test_coord_helper_da_latlon): """Test assign_latitude_longitude in basic usage on Dataset.""" @@ -922,7 +922,7 @@ def test_assign_latitude_longitude_basic_dataset(test_coord_helper_da_yx, lon.values, 3) -@needs_cartopy +@needs_pyproj def test_assign_y_x_basic_dataarray(test_coord_helper_da_yx, test_coord_helper_da_latlon): """Test assign_y_x in basic usage on DataArray.""" new_da = test_coord_helper_da_latlon.metpy.assign_y_x() @@ -939,7 +939,7 @@ def test_assign_y_x_error_existing_dataarray( assert 'y/x coordinate(s) are present' in str(exc) -@needs_cartopy +@needs_pyproj def test_assign_y_x_force_existing_dataarray( test_coord_helper_da_dummy_yx, test_coord_helper_da_yx): """Test assign_y_x with existing coordinates forcing new.""" @@ -949,7 +949,7 @@ def test_assign_y_x_force_existing_dataarray( np.testing.assert_array_almost_equal(test_coord_helper_da_yx['x'].values, x.values, 3) -@needs_cartopy +@needs_pyproj def test_assign_y_x_dataarray_outside_tolerance(test_coord_helper_da_latlon): """Test assign_y_x raises ValueError when tolerance is exceeded on DataArray.""" with pytest.raises(ValueError) as exc: @@ -957,7 +957,7 @@ def test_assign_y_x_dataarray_outside_tolerance(test_coord_helper_da_latlon): assert 'cannot be collapsed to 1D within tolerance' in str(exc) -@needs_cartopy +@needs_pyproj def test_assign_y_x_dataarray_transposed(test_coord_helper_da_yx, test_coord_helper_da_latlon): """Test assign_y_x on DataArray with transposed order.""" new_da = test_coord_helper_da_latlon.transpose(transpose_coords=True).metpy.assign_y_x() @@ -966,7 +966,7 @@ def test_assign_y_x_dataarray_transposed(test_coord_helper_da_yx, test_coord_hel np.testing.assert_array_almost_equal(test_coord_helper_da_yx['x'].values, x.values, 3) -@needs_cartopy +@needs_pyproj def test_assign_y_x_dataset_assumed_order(test_coord_helper_da_yx, test_coord_helper_da_latlon): """Test assign_y_x on Dataset where order must be assumed.""" diff --git a/tutorials/xarray_tutorial.py b/tutorials/xarray_tutorial.py index d64fcf7094d..ff716bb966d 100644 --- a/tutorials/xarray_tutorial.py +++ b/tutorials/xarray_tutorial.py @@ -150,15 +150,21 @@ # Getting the cartopy coordinate reference system (CRS) of the projection of a DataArray is as # straightforward as using the ``data_var.metpy.cartopy_crs`` property: -data_crs = data['temperature'].metpy.cartopy_crs -print(data_crs) +cartopy_crs = data['temperature'].metpy.cartopy_crs +print(cartopy_crs) + +######################################################################### +# Likewise, the PyProj CRS can be obtained with the ``.pyproj_crs`` property: + +pyproj_crs = data['temperature'].metpy.pyproj_crs +print(pyproj_crs) ######################################################################### # The cartopy ``Globe`` can similarly be accessed via the ``data_var.metpy.cartopy_globe`` # property: -data_globe = data['temperature'].metpy.cartopy_globe -print(data_globe) +cartopy_globe = data['temperature'].metpy.cartopy_globe +print(cartopy_globe) ######################################################################### # Calculations @@ -176,7 +182,7 @@ lat, lon = xr.broadcast(y, x) f = mpcalc.coriolis_parameter(lat) -dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat, initstring=data_crs.proj4_init) +dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat, geod=pyproj_crs.get_geod()) heights = data['height'].metpy.loc[{'time': time[0], 'vertical': 500. * units.hPa}] u_geo, v_geo = mpcalc.geostrophic_wind(heights, f, dx, dy) print(u_geo) @@ -240,7 +246,7 @@ # Let's add a projection and coastlines to it ax = plt.axes(projection=ccrs.LambertConformal()) data['height'].metpy.loc[{'time': time[0], - 'vertical': 500. * units.hPa}].plot(ax=ax, transform=data_crs) + 'vertical': 500. * units.hPa}].plot(ax=ax, transform=cartopy_crs) ax.coastlines() plt.show() @@ -252,7 +258,7 @@ data_level = data.metpy.loc[{time.name: time[0], vertical.name: 500. * units.hPa}] # Create the matplotlib figure and axis -fig, ax = plt.subplots(1, 1, figsize=(12, 8), subplot_kw={'projection': data_crs}) +fig, ax = plt.subplots(1, 1, figsize=(12, 8), subplot_kw={'projection': cartopy_crs}) # Plot RH as filled contours rh = ax.contourf(x, y, data_level['relative_humidity'], levels=[70, 80, 90, 100],