Merge pull request #786 from dopplershift/xarray-projections

Xarray + CF + CartoPy Projection Handling

.codeclimate.yml

@@ -9,6 +9,9 @@ checks:
   argument-count:
     config:
       threshold: 10
+  method-complexity:
+    config:
+      threshold: 15
 
 plugins:
   sonar-python:

.codecov.yml

@@ -22,7 +22,9 @@ coverage:
 
       tests:
         target: 100%
-        paths: "metpy/.*/tests/.*"
+        paths:
+          - "metpy/.*/tests/.*"
+          - "metpy/tests/.*"
 
   notify:
     gitter:

.travis.yml

@@ -48,7 +48,7 @@ matrix:
   include:
     - python: 2.7
       env:
-        - VERSIONS="numpy==1.10.0 matplotlib==1.4.0 scipy==0.14.0 pint==0.8"
+        - VERSIONS="numpy==1.10.0 matplotlib==1.4.0 scipy==0.14.0 pint==0.8 xarray==0.9.0"
         - TASK="coverage"
         - TEST_OUTPUT_CONTROL=""
     - python: 3.4

examples/XArray_Projections.py

@@ -0,0 +1,35 @@
+# Copyright (c) 2018 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""
+XArray Projection Handling
+==========================
+
+Use MetPy's XArray accessors to simplify opening a data file and plotting
+data on a map using CartoPy.
+"""
+import cartopy.feature as cfeature
+import matplotlib.pyplot as plt
+import xarray as xr
+
+# Ensure xarray accessors are available
+import metpy.io  # noqa: F401
+from metpy.testing import get_test_data
+
+ds = xr.open_dataset(get_test_data('narr_example.nc', as_file_obj=False))
+data_var = ds.metpy.parse_cf('Temperature')
+
+x = data_var.x
+y = data_var.y
+im_data = data_var.isel(time=0).sel(isobaric=1000.)
+
+fig = plt.figure(figsize=(14, 14))
+ax = fig.add_subplot(1, 1, 1, projection=data_var.metpy.cartopy_crs)
+
+ax.imshow(im_data, extent=(x.min(), x.max(), y.min(), y.max()),
+          cmap='RdBu', origin='lower' if y[0] < y[-1] else 'upper')
+ax.coastlines(color='tab:green', resolution='10m')
+ax.add_feature(cfeature.LAKES.with_scale('10m'), facecolor='none', edgecolor='tab:blue')
+ax.add_feature(cfeature.RIVERS.with_scale('10m'), edgecolor='tab:blue')
+
+plt.show()

examples/formats/GINI_Water_Vapor.py

@@ -8,7 +8,6 @@
 Use MetPy's support for GINI files to read in a water vapor satellite image and plot the
 data using CartoPy.
 """
-import cartopy.crs as ccrs
 import cartopy.feature as cfeature
 import matplotlib.pyplot as plt
 import xarray as xr
@@ -24,33 +23,21 @@
 print(f)
 
 ###########################################
-
 # Get a Dataset view of the data (essentially a NetCDF-like interface to the
-# underlying data). Pull out the data, (x, y) coordinates, and the projection
-# information.
+# underlying data). Pull out the data and (x, y) coordinates. We use `metpy.parse_cf` to
+# handle parsing some netCDF Climate and Forecasting (CF) metadata to simplify working with
+# projections.
 ds = xr.open_dataset(f)
 x = ds.variables['x'][:]
 y = ds.variables['y'][:]
-dat = ds.variables['WV']
-proj_var = ds.variables[dat.attrs['grid_mapping']]
-print(proj_var)
+dat = ds.metpy.parse_cf('WV')
 
 ###########################################
-
-# Create CartoPy projection information for the file
-globe = ccrs.Globe(ellipse='sphere', semimajor_axis=proj_var.attrs['earth_radius'],
-                   semiminor_axis=proj_var.attrs['earth_radius'])
-proj = ccrs.LambertConformal(central_longitude=proj_var.attrs['longitude_of_central_meridian'],
-                             central_latitude=proj_var.attrs['latitude_of_projection_origin'],
-                             standard_parallels=[proj_var.attrs['standard_parallel']],
-                             globe=globe)
-
-###########################################
-
-# Plot the image
+# Plot the image. We use MetPy's xarray/cartopy integration to automatically handle parsing
+# the projection information.
 fig = plt.figure(figsize=(10, 12))
 add_metpy_logo(fig, 125, 145)
-ax = fig.add_subplot(1, 1, 1, projection=proj)
+ax = fig.add_subplot(1, 1, 1, projection=dat.metpy.cartopy_crs)
 wv_norm, wv_cmap = colortables.get_with_range('WVCIMSS', 100, 260)
 wv_cmap.set_under('k')
 im = ax.imshow(dat[:], cmap=wv_cmap, norm=wv_norm,

metpy/calc/kinematics.py

@@ -8,7 +8,6 @@
 import warnings
 
 import numpy as np
-from pyproj import Geod
 
 from ..calc import coriolis_parameter
 from ..calc.tools import first_derivative, get_layer_heights, gradient
@@ -698,6 +697,8 @@ def lat_lon_grid_deltas(longitude, latitude, **kwargs):
     Assumes [Y, X] for 2D arrays
 
     """
+    from pyproj import Geod
+
     # Inputs must be the same number of dimensions
     if latitude.ndim != longitude.ndim:
         raise ValueError('Latitude and longitude must have the same number of dimensions.')

metpy/cbook.py

@@ -1,4 +1,4 @@
-# Copyright (c) 2008,2015 MetPy Developers.
+# Copyright (c) 2008,2015,2018 MetPy Developers.
 # Distributed under the terms of the BSD 3-Clause License.
 # SPDX-License-Identifier: BSD-3-Clause
 """Collection of generally useful utility code from the cookbook."""
@@ -39,4 +39,39 @@ def get_test_data(fname, as_file_obj=True):
     return path
 
 
-__all__ = ('get_test_data', 'is_string_like', 'iterable')
+class Registry(object):
+    """Provide a generic function registry.
+
+    This provides a class to instantiate, which then has a `register` method that can
+    be used as a decorator on functions to register them under a particular name.
+    """
+
+    def __init__(self):
+        """Initialize an empty registry."""
+        self._registry = {}
+
+    def register(self, name):
+        """Register a callable with the registry under a particular name.
+
+        Parameters
+        ----------
+        name : str
+            The name under which to register a function
+
+        Returns
+        -------
+        dec : callable
+            A decorator that takes a function and will register it under the name.
+
+        """
+        def dec(func):
+            self._registry[name] = func
+            return func
+        return dec
+
+    def __getitem__(self, name):
+        """Return any callable registered under name."""
+        return self._registry[name]
+
+
+__all__ = ('Registry', 'get_test_data', 'is_string_like', 'iterable')

metpy/io/__init__.py

@@ -1,4 +1,4 @@
-# Copyright (c) 2015,2016 MetPy Developers.
+# Copyright (c) 2015,2016,2018 MetPy Developers.
 # Distributed under the terms of the BSD 3-Clause License.
 # SPDX-License-Identifier: BSD-3-Clause
 """MetPy's IO module contains classes for reading files. These classes are written
@@ -9,6 +9,8 @@
 
 from .gini import *  # noqa: F403
 from .nexrad import *  # noqa: F403
+from .xarray import *  # noqa: F403
 
 __all__ = gini.__all__[:]  # pylint: disable=undefined-variable
 __all__.extend(nexrad.__all__)  # pylint: disable=undefined-variable
+__all__.extend(xarray.__all__)  # pylint: disable=undefined-variable

metpy/io/gini.py

@@ -2,6 +2,7 @@
 # Distributed under the terms of the BSD 3-Clause License.
 # SPDX-License-Identifier: BSD-3-Clause
 """Tools to process GINI-formatted products."""
+from __future__ import absolute_import
 
 import contextlib
 from datetime import datetime

metpy/io/tests/test_xarray.py

@@ -0,0 +1,108 @@
+# Copyright (c) 2018 MetPy Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+"""Test the operation of MetPy's XArray accessors."""
+
+import cartopy.crs as ccrs
+import numpy as np
+import pytest
+import xarray as xr
+
+# Get activate the xarray accessors
+import metpy.io  # noqa: F401
+from metpy.testing import assert_almost_equal, get_test_data
+from metpy.units import units
+
+
+@pytest.fixture
+def test_ds():
+    """Provide an xarray dataset for testing."""
+    return xr.open_dataset(get_test_data('narr_example.nc', as_file_obj=False))
+
+
+@pytest.fixture
+def test_var():
+    """Provide a standard, parsed, variable for testing."""
+    ds = test_ds()
+    return ds.metpy.parse_cf('Temperature')
+
+
+def test_projection(test_var):
+    """Test getting the proper projection out of the variable."""
+    crs = test_var.metpy.crs
+    assert crs['grid_mapping_name'] == 'lambert_conformal_conic'
+
+    assert isinstance(test_var.metpy.cartopy_crs, ccrs.LambertConformal)
+
+
+def test_no_projection(test_ds):
+    """Test getting the crs attribute when not available produces a sensible error."""
+    var = test_ds.lat
+    with pytest.raises(AttributeError) as exc:
+        var.metpy.crs
+
+    assert 'not available' in str(exc.value)
+
+
+def test_units(test_var):
+    """Test unit handling through the accessor."""
+    arr = test_var.metpy.unit_array
+    assert isinstance(arr, units.Quantity)
+    assert arr.units == units.kelvin
+
+
+def test_convert_units(test_var):
+    """Test in-place conversion of units."""
+    test_var.metpy.convert_units('degC')
+
+    # Check that variable metadata is updated
+    assert test_var.attrs['units'] == 'degC'
+
+    # Make sure we now get an array back with properly converted values
+    assert test_var.metpy.unit_array.units == units.degC
+    assert_almost_equal(test_var[0, 0, 0, 0], 18.44 * units.degC, 2)
+
+
+def test_radian_projection_coords():
+    """Test fallback code for radian units in projection coordinate variables."""
+    proj = xr.DataArray(0, attrs={'grid_mapping_name': 'geostationary',
+                                  'perspective_point_height': 3})
+    x = xr.DataArray(np.arange(3),
+                     attrs={'standard_name': 'projection_x_coordinate', 'units': 'radians'})
+    y = xr.DataArray(np.arange(2),
+                     attrs={'standard_name': 'projection_y_coordinate', 'units': 'radians'})
+    data = xr.DataArray(np.arange(6).reshape(2, 3), coords=(y, x), dims=('y', 'x'),
+                        attrs={'grid_mapping': 'fixedgrid_projection'})
+    ds = xr.Dataset({'data': data, 'fixedgrid_projection': proj})
+
+    # Check that the coordinates in this case are properly converted
+    data_var = ds.metpy.parse_cf('data')
+    assert data_var.coords['x'].metpy.unit_array[1] == 3 * units.meter
+    assert data_var.coords['y'].metpy.unit_array[1] == 3 * units.meter
+
+
+def test_missing_grid_mapping():
+    """Test falling back to implicit lat/lon projection."""
+    lon = xr.DataArray(-np.arange(3),
+                       attrs={'standard_name': 'longitude', 'units': 'degrees_east'})
+    lat = xr.DataArray(np.arange(2),
+                       attrs={'standard_name': 'latitude', 'units': 'degrees_north'})
+    data = xr.DataArray(np.arange(6).reshape(2, 3), coords=(lat, lon), dims=('y', 'x'))
+    ds = xr.Dataset({'data': data})
+
+    data_var = ds.metpy.parse_cf('data')
+    assert 'crs' in data_var.coords
+
+
+def test_missing_grid_mapping_var():
+    """Test behavior when we can't find the variable pointed to by grid_mapping."""
+    x = xr.DataArray(np.arange(3),
+                     attrs={'standard_name': 'projection_x_coordinate', 'units': 'radians'})
+    y = xr.DataArray(np.arange(2),
+                     attrs={'standard_name': 'projection_y_coordinate', 'units': 'radians'})
+    data = xr.DataArray(np.arange(6).reshape(2, 3), coords=(y, x), dims=('y', 'x'),
+                        attrs={'grid_mapping': 'fixedgrid_projection'})
+    ds = xr.Dataset({'data': data})
+
+    with pytest.warns(UserWarning, match='Could not find'):
+        ds.metpy.parse_cf('data')