From 579e22963df0803b629457c10a14be7457266ffd Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 23 Mar 2018 12:13:13 -0600 Subject: [PATCH 01/10] ENH: Add simple registry implementation This adds a class with a register decorator that can be used to register callables under a string. --- metpy/cbook.py | 39 +++++++++++++++++++++++++++++++++++++-- metpy/tests/test_cbook.py | 16 ++++++++++++++++ 2 files changed, 53 insertions(+), 2 deletions(-) create mode 100644 metpy/tests/test_cbook.py diff --git a/metpy/cbook.py b/metpy/cbook.py index 5c56f6a612d..d9a8622164f 100644 --- a/metpy/cbook.py +++ b/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') diff --git a/metpy/tests/test_cbook.py b/metpy/tests/test_cbook.py new file mode 100644 index 00000000000..c94713a4526 --- /dev/null +++ b/metpy/tests/test_cbook.py @@ -0,0 +1,16 @@ +# Copyright (c) 2018 MetPy Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +"""Test functionality of MetPy's utility code.""" + +from metpy.cbook import Registry + + +def test_registry(): + """Test that the registry properly registers things.""" + reg = Registry() + + a = 'foo' + reg.register('mine')(a) + + assert reg['mine'] is a From 3e0b99fa4ff892b12d522b5f2dd4f68ed7bd0cf4 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 23 Mar 2018 13:53:55 -0600 Subject: [PATCH 02/10] ENH: Add framework for converting CF attributes to CartoPy --- metpy/plots/mapping.py | 159 ++++++++++++++++++++++++++ metpy/plots/tests/test_mapping.py | 181 ++++++++++++++++++++++++++++++ 2 files changed, 340 insertions(+) create mode 100644 metpy/plots/mapping.py create mode 100644 metpy/plots/tests/test_mapping.py diff --git a/metpy/plots/mapping.py b/metpy/plots/mapping.py new file mode 100644 index 00000000000..c298f139671 --- /dev/null +++ b/metpy/plots/mapping.py @@ -0,0 +1,159 @@ +# Copyright (c) 2018 MetPy Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +"""Tools to help with mapping/geographic applications. + +Currently this includes tools for working with CartoPy projections. + +""" +import cartopy.crs as ccrs + +from ..cbook import Registry + + +class CFProjection(object): + """Handle parsing CF projection metadata.""" + + _default_attr_mapping = [('false_easting', 'false_easting'), + ('false_northing', 'false_northing'), + ('central_latitude', 'latitude_of_projection_origin'), + ('central_longitude', 'longitude_of_projection_origin')] + + projection_registry = Registry() + + def __init__(self, attrs): + """Initialize the CF Projection handler with a set of metadata attributes.""" + self._attrs = attrs + + @classmethod + def register(cls, name): + """Register a new projection to handle.""" + return cls.projection_registry.register(name) + + @classmethod + def build_projection_kwargs(cls, source, mapping): + """Handle mapping a dictionary of metadata to keyword arguments.""" + return cls._map_arg_names(source, cls._default_attr_mapping + mapping) + + @staticmethod + def _map_arg_names(source, mapping): + """Map one set of keys to another.""" + return {cartopy_name: source[cf_name] for cartopy_name, cf_name in mapping + if cf_name in source} + + def _make_cartopy_globe(self): + """Initialize a `cartopy.crs.Globe` from the metadata.""" + if 'earth_radius' in self._attrs: + kwargs = {'ellipse': 'sphere', 'semimajor_axis': self._attrs['earth_radius'], + 'semiminor_axis': self._attrs['earth_radius']} + else: + attr_mapping = [('semimajor_axis', 'semi_major_axis'), + ('semiminor_axis', 'semi_minor_axis'), + ('inverse_flattening', 'inverse_flattening')] + kwargs = self._map_arg_names(self._attrs, attr_mapping) + + # WGS84 with semi_major==semi_minor is NOT the same as spherical Earth + # Also need to handle the case where we're not given any spheroid + kwargs['ellipse'] = None if kwargs else 'sphere' + + return ccrs.Globe(**kwargs) + + def to_cartopy(self): + """Convert to a CartoPy projection.""" + globe = self._make_cartopy_globe() + proj_name = self._attrs['grid_mapping_name'] + try: + proj_handler = self.projection_registry[proj_name] + except KeyError: + raise ValueError('Unhandled projection: {}'.format(proj_name)) + + return proj_handler(self._attrs, globe) + + def to_dict(self): + """Get the dictionary of metadata attributes.""" + return self._attrs.copy() + + def __str__(self): + """Get a string representation of the projection.""" + return 'Projection: ' + self._attrs['grid_mapping_name'] + + def __getitem__(self, item): + """Return a given attribute.""" + return self._attrs[item] + + +@CFProjection.register('geostationary') +def make_geo(attrs_dict, globe): + """Handle geostationary projection.""" + attr_mapping = [('satellite_height', 'perspective_point_height'), + ('sweep_axis', 'sweep_angle_axis')] + kwargs = CFProjection.build_projection_kwargs(attrs_dict, attr_mapping) + + # CartoPy can't handle central latitude for Geostationary (nor should it) + # Just remove it if it's 0. + if not kwargs.get('central_latitude'): + kwargs.pop('central_latitude', None) + + # If sweep_angle_axis is not present, we should look for fixed_angle_axis and adjust + if 'sweep_axis' not in kwargs: + kwargs['sweep_axis'] = 'x' if attrs_dict['fixed_angle_axis'] == 'y' else 'y' + + return ccrs.Geostationary(globe=globe, **kwargs) + + +@CFProjection.register('lambert_conformal_conic') +def make_lcc(attrs_dict, globe): + """Handle Lambert conformal conic projection.""" + attr_mapping = [('central_longitude', 'longitude_of_central_meridian'), + ('standard_parallels', 'standard_parallel')] + kwargs = CFProjection.build_projection_kwargs(attrs_dict, attr_mapping) + if 'standard_parallels' in kwargs: + try: + len(kwargs['standard_parallels']) + except TypeError: + kwargs['standard_parallels'] = [kwargs['standard_parallels']] + return ccrs.LambertConformal(globe=globe, **kwargs) + + +@CFProjection.register('latitude_longitude') +def make_latlon(attrs_dict, globe): + """Handle plain latitude/longitude mapping.""" + # TODO: Really need to use Geodetic to pass the proper globe + return ccrs.PlateCarree() + + +@CFProjection.register('mercator') +def make_mercator(attrs_dict, globe): + """Handle Mercator projection.""" + attr_mapping = [('latitude_true_scale', 'standard_parallel'), + ('scale_factor', 'scale_factor_at_projection_origin')] + kwargs = CFProjection.build_projection_kwargs(attrs_dict, attr_mapping) + + # Work around the fact that in CartoPy <= 0.16 can't handle the easting/northing + # in Mercator + if not kwargs.get('false_easting'): + kwargs.pop('false_easting', None) + if not kwargs.get('false_northing'): + kwargs.pop('false_northing', None) + + return ccrs.Mercator(globe=globe, **kwargs) + + +@CFProjection.register('stereographic') +def make_stereo(attrs_dict, globe): + """Handle generic stereographic projection.""" + attr_mapping = [('scale_factor', 'scale_factor_at_projection_origin')] + kwargs = CFProjection.build_projection_kwargs(attrs_dict, attr_mapping) + + return ccrs.Stereographic(globe=globe, **kwargs) + + +@CFProjection.register('polar_stereographic') +def make_polar_stereo(attrs_dict, globe): + """Handle polar stereographic projection.""" + attr_mapping = [('central_longitude', 'straight_vertical_longitude_from_pole'), + ('true_scale_latitude', 'standard_parallel'), + ('scale_factor', 'scale_factor_at_projection_origin')] + kwargs = CFProjection.build_projection_kwargs(attrs_dict, attr_mapping) + + return ccrs.Stereographic(globe=globe, **kwargs) diff --git a/metpy/plots/tests/test_mapping.py b/metpy/plots/tests/test_mapping.py new file mode 100644 index 00000000000..4b46399a126 --- /dev/null +++ b/metpy/plots/tests/test_mapping.py @@ -0,0 +1,181 @@ +# Copyright (c) 2018 MetPy Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +"""Test the handling of various mapping tasks.""" + +import cartopy.crs as ccrs +import pytest + +from metpy.plots.mapping import CFProjection + + +def test_cfprojection_arg_mapping(): + """Test the projection mapping arguments.""" + source = {'source': 'a', 'longitude_of_projection_origin': -100} + + # 'dest' should be the argument in the output, with the value from source + mapping = [('dest', 'source')] + + kwargs = CFProjection.build_projection_kwargs(source, mapping) + assert kwargs == {'dest': 'a', 'central_longitude': -100} + + +def test_cfprojection_api(): + """Test the basic API of the projection interface.""" + attrs = {'grid_mapping_name': 'lambert_conformal_conic', 'earth_radius': 6367000} + proj = CFProjection(attrs) + + assert proj['earth_radius'] == 6367000 + assert proj.to_dict() == attrs + assert str(proj) == 'Projection: lambert_conformal_conic' + + +def test_bad_projection_raises(): + """Test behavior when given an unknown projection.""" + attrs = {'grid_mapping_name': 'unknown'} + with pytest.raises(ValueError) as exc: + CFProjection(attrs).to_cartopy() + + assert 'Unhandled projection' in str(exc.value) + + +def test_globe(): + """Test handling building a cartopy globe.""" + attrs = {'grid_mapping_name': 'lambert_conformal_conic', 'earth_radius': 6367000, + 'standard_parallel': 25} + proj = CFProjection(attrs) + + crs = proj.to_cartopy() + globe_params = crs.globe.to_proj4_params() + + assert globe_params['ellps'] == 'sphere' + assert globe_params['a'] == 6367000 + assert globe_params['b'] == 6367000 + + +def test_globe_spheroid(): + """Test handling building a cartopy globe that is not spherical.""" + attrs = {'grid_mapping_name': 'lambert_conformal_conic', 'semi_major_axis': 6367000, + 'semi_minor_axis': 6360000} + proj = CFProjection(attrs) + + crs = proj.to_cartopy() + globe_params = crs.globe.to_proj4_params() + + assert 'ellps' not in globe_params + assert globe_params['a'] == 6367000 + assert globe_params['b'] == 6360000 + + +def test_lcc(): + """Test handling lambert conformal conic projection.""" + attrs = {'grid_mapping_name': 'lambert_conformal_conic', 'earth_radius': 6367000, + 'standard_parallel': [25, 30]} + proj = CFProjection(attrs) + + crs = proj.to_cartopy() + assert isinstance(crs, ccrs.LambertConformal) + assert crs.proj4_params['lat_1'] == 25 + assert crs.proj4_params['lat_2'] == 30 + assert crs.globe.to_proj4_params()['ellps'] == 'sphere' + + +def test_lcc_minimal(): + """Test handling lambert conformal conic projection with minimal attributes.""" + attrs = {'grid_mapping_name': 'lambert_conformal_conic'} + crs = CFProjection(attrs).to_cartopy() + assert isinstance(crs, ccrs.LambertConformal) + + +def test_lcc_single_std_parallel(): + """Test lambert conformal projection with one standard parallel.""" + attrs = {'grid_mapping_name': 'lambert_conformal_conic', 'standard_parallel': 25} + crs = CFProjection(attrs).to_cartopy() + assert isinstance(crs, ccrs.LambertConformal) + assert crs.proj4_params['lat_1'] == 25 + + +def test_mercator(): + """Test handling a mercator projection.""" + attrs = {'grid_mapping_name': 'mercator', 'standard_parallel': 25, + 'longitude_of_projection_origin': -100, 'false_easting': 0, 'false_westing': 0} + crs = CFProjection(attrs).to_cartopy() + + assert isinstance(crs, ccrs.Mercator) + assert crs.proj4_params['lat_ts'] == 25 + assert crs.proj4_params['lon_0'] == -100 + + +# This won't work until at least CartoPy > 0.16.0 +# def test_mercator_scale_factor(): +# """Test handling a mercator projection with a scale factor.""" +# attrs = {'grid_mapping_name': 'mercator', 'scale_factor_at_projection_origin': 0.9} +# crs = CFProjection(attrs).to_cartopy() +# +# assert isinstance(crs, ccrs.Mercator) +# assert crs.proj4_params['k_0'] == 0.9 + + +def test_geostationary(): + """Test handling a geostationary projection.""" + attrs = {'grid_mapping_name': 'geostationary', 'perspective_point_height': 35000000, + 'longitude_of_projection_origin': -100, 'sweep_angle_axis': 'x', + 'latitude_of_projection_origin': 0} + crs = CFProjection(attrs).to_cartopy() + + assert isinstance(crs, ccrs.Geostationary) + assert crs.proj4_params['h'] == 35000000 + assert crs.proj4_params['lon_0'] == -100 + assert crs.proj4_params['sweep'] == 'x' + + +def test_geostationary_fixed_angle(): + """Test handling geostationary information that gives fixed angle instead of sweep.""" + attrs = {'grid_mapping_name': 'geostationary', 'fixed_angle_axis': 'y'} + crs = CFProjection(attrs).to_cartopy() + + assert isinstance(crs, ccrs.Geostationary) + assert crs.proj4_params['sweep'] == 'x' + + +def test_stereographic(): + """Test handling a stereographic projection.""" + attrs = {'grid_mapping_name': 'stereographic', 'scale_factor_at_projection_origin': 0.9, + 'longitude_of_projection_origin': -100, 'latitude_of_projection_origin': 60} + crs = CFProjection(attrs).to_cartopy() + + assert isinstance(crs, ccrs.Stereographic) + assert crs.proj4_params['lon_0'] == -100 + assert crs.proj4_params['lat_0'] == 60 + assert crs.proj4_params['k_0'] == 0.9 + + +def test_polar_stereographic(): + """Test handling a polar stereographic projection.""" + attrs = {'grid_mapping_name': 'polar_stereographic', 'latitude_of_projection_origin': 90, + 'scale_factor_at_projection_origin': 0.9, + 'straight_vertical_longitude_from_pole': -100, } + crs = CFProjection(attrs).to_cartopy() + + assert isinstance(crs, ccrs.Stereographic) + assert crs.proj4_params['lon_0'] == -100 + assert crs.proj4_params['lat_0'] == 90 + assert crs.proj4_params['k_0'] == 0.9 + + +def test_polar_stereographic_std_parallel(): + """Test handling a polar stereographic projection that gives a standard parallel.""" + attrs = {'grid_mapping_name': 'polar_stereographic', 'latitude_of_projection_origin': -90, + 'standard_parallel': 60} + crs = CFProjection(attrs).to_cartopy() + + assert isinstance(crs, ccrs.Stereographic) + assert crs.proj4_params['lat_0'] == -90 + assert crs.proj4_params['lat_ts'] == 60 + + +def test_lat_lon(): + """Test handling basic lat/lon projection.""" + attrs = {'grid_mapping_name': 'latitude_longitude'} + crs = CFProjection(attrs).to_cartopy() + assert isinstance(crs, ccrs.PlateCarree) From 68c86c82783dcb18e80ab6c531b5be90d563950d Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 23 Mar 2018 13:54:15 -0600 Subject: [PATCH 03/10] ENH: Add xarray accessors to simplify getting projection --- metpy/io/__init__.py | 4 +- metpy/io/tests/test_xarray.py | 108 ++++++++++++++++++++++++++++++ metpy/io/xarray.py | 120 ++++++++++++++++++++++++++++++++++ 3 files changed, 231 insertions(+), 1 deletion(-) create mode 100644 metpy/io/tests/test_xarray.py create mode 100644 metpy/io/xarray.py diff --git a/metpy/io/__init__.py b/metpy/io/__init__.py index 09405cd9a4c..be88b402f78 100644 --- a/metpy/io/__init__.py +++ b/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 diff --git a/metpy/io/tests/test_xarray.py b/metpy/io/tests/test_xarray.py new file mode 100644 index 00000000000..558a74098f0 --- /dev/null +++ b/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') diff --git a/metpy/io/xarray.py b/metpy/io/xarray.py new file mode 100644 index 00000000000..fb191bbc5a1 --- /dev/null +++ b/metpy/io/xarray.py @@ -0,0 +1,120 @@ +# Copyright (c) 2018 MetPy Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +"""Provide accessors to enhance interoperability between XArray and MetPy.""" +import xarray as xr + +from ..units import DimensionalityError, units + +__all__ = [] + + +# TODO: Should we be providing renaming/mapping for easy access to x/y coordinates? +@xr.register_dataarray_accessor('metpy') +class MetPyAccessor(object): + """Provide custom attributes and methods on XArray DataArray for MetPy functionality.""" + + def __init__(self, data_array): + """Initialize accessor with a DataArray.""" + self._data_array = data_array + self._units = self._data_array.attrs.get('units', 'dimensionless') + + @property + def unit_array(self): + """Return data values as a `pint.Quantity`.""" + return self._data_array.values * units(self._units) + + @unit_array.setter + def unit_array(self, values): + """Set data values as a `pint.Quantity`.""" + self._data_array.values = values + self._units = self._data_array.attrs['units'] = str(values.units) + + def convert_units(self, units): + """Convert the data values to different units in-place.""" + self.unit_array = self.unit_array.to(units) + + @property + def crs(self): + """Provide easy access to the `crs` coordinate.""" + if 'crs' in self._data_array.coords: + return self._data_array.coords['crs'].item() + raise AttributeError('crs attribute is not available due to lack of crs coordinate.') + + @property + def cartopy_crs(self): + """Return the coordinate reference system (CRS) as a cartopy object.""" + return self.crs.to_cartopy() + + +@xr.register_dataset_accessor('metpy') +class CFConventionHandler(object): + """Provide custom attributes and methods on XArray Dataset for MetPy functionality.""" + + def __init__(self, dataset): + """Initialize accessor with a Dataset.""" + self._dataset = dataset + + def parse_cf(self, varname): + """Parse Climate and Forecasting (CF) convention metadata.""" + from ..plots.mapping import CFProjection + + var = self._dataset[varname] + if 'grid_mapping' in var.attrs: + proj_name = var.attrs['grid_mapping'] + try: + proj_var = self._dataset.variables[proj_name] + except KeyError: + import warnings + warnings.warn( + 'Could not find variable corresponding to the value of ' + 'grid_mapping: {}'.format(proj_name)) + else: + var.coords['crs'] = CFProjection(proj_var.attrs) + var.attrs.pop('grid_mapping') + self._fixup_coords(var) + + # Trying to guess whether we should be adding a crs to this variable's coordinates + # First make sure it's missing CRS but isn't lat/lon itself + if not (self._check_lat(var) or self._check_lon(var)) and 'crs' not in var.coords: + # Look for both lat/lon in the coordinates + has_lat = has_lon = False + for coord_var in var.coords.values(): + has_lat = has_lat or self._check_lat(coord_var) + has_lon = has_lon or self._check_lon(coord_var) + + # If we found them, create a lat/lon projection as the crs coord + if has_lat and has_lon: + var.coords['crs'] = CFProjection({'grid_mapping_name': 'latitude_longitude'}) + + return var + + @staticmethod + def _check_lat(var): + if var.attrs.get('standard_name') == 'latitude': + return True + + units = var.attrs.get('units', '').replace('degrees', 'degree') + return units in {'degree_north', 'degree_N', 'degreeN'} + + @staticmethod + def _check_lon(var): + if var.attrs.get('standard_name') == 'longitude': + return True + + units = var.attrs.get('units', '').replace('degrees', 'degree') + return units in {'degree_east', 'degree_E', 'degreeE'} + + def _fixup_coords(self, var): + """Clean up the units on the coordinate variables.""" + for coord_name, data_array in var.coords.items(): + if data_array.attrs.get('standard_name') in ('projection_x_coordinate', + 'projection_y_coordinate'): + try: + var.coords[coord_name].metpy.convert_units('meters') + except DimensionalityError: # Radians! + new_data_array = data_array.copy() + height = var.coords['crs'].item()['perspective_point_height'] + scaled_vals = new_data_array.metpy.unit_array * (height * units.meters) + new_data_array.metpy.unit_array = scaled_vals.to('meters') + var.coords[coord_name] = new_data_array From e100e6be59d5c78a600637d89399c55f39242918 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 23 Mar 2018 13:54:38 -0600 Subject: [PATCH 04/10] ENH: Add example of using xarray projection info --- examples/XArray_Projections.py | 35 ++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 examples/XArray_Projections.py diff --git a/examples/XArray_Projections.py b/examples/XArray_Projections.py new file mode 100644 index 00000000000..8cd9a88a30d --- /dev/null +++ b/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() From 7b92265cd42fcaf43c5b2ec8afbc8e99d731e98f Mon Sep 17 00:00:00 2001 From: Ryan May Date: Mon, 26 Mar 2018 10:58:14 -0600 Subject: [PATCH 05/10] MNT: Bump up complexity threshold for CodeClimate --- .codeclimate.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.codeclimate.yml b/.codeclimate.yml index 44e0cd15c4c..9c5c77ddd57 100644 --- a/.codeclimate.yml +++ b/.codeclimate.yml @@ -9,6 +9,9 @@ checks: argument-count: config: threshold: 10 + method-complexity: + config: + threshold: 15 plugins: sonar-python: From fb6e7aff9a30f0a1efbf7bf28735e985aaf50fe9 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 11 May 2018 14:54:15 -0600 Subject: [PATCH 06/10] ENH: Use xarray support to simplify GINI example --- examples/formats/GINI_Water_Vapor.py | 27 +++++++-------------------- 1 file changed, 7 insertions(+), 20 deletions(-) diff --git a/examples/formats/GINI_Water_Vapor.py b/examples/formats/GINI_Water_Vapor.py index f67580f483f..2088918af7b 100644 --- a/examples/formats/GINI_Water_Vapor.py +++ b/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, From 045f2ff949f7140136f0fb54f883fce062e8b6a3 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 11 May 2018 16:21:42 -0600 Subject: [PATCH 07/10] MNT: Fix import on Python 2.7 Default behavior of Python 2.7 gets confused by our having our own xarray module and trying to import the main library. Some future imports solve this... --- metpy/io/gini.py | 1 + metpy/io/xarray.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/metpy/io/gini.py b/metpy/io/gini.py index f541a9e9cd1..1dcc8eda33e 100644 --- a/metpy/io/gini.py +++ b/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 diff --git a/metpy/io/xarray.py b/metpy/io/xarray.py index fb191bbc5a1..88760eced6c 100644 --- a/metpy/io/xarray.py +++ b/metpy/io/xarray.py @@ -2,6 +2,8 @@ # Distributed under the terms of the BSD 3-Clause License. # SPDX-License-Identifier: BSD-3-Clause """Provide accessors to enhance interoperability between XArray and MetPy.""" +from __future__ import absolute_import + import xarray as xr from ..units import DimensionalityError, units From 2d500b49d71a88091a691a92ef64653517fc2de7 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 11 May 2018 16:47:25 -0600 Subject: [PATCH 08/10] MNT: Move pyproj import Right now metpy.calc won't import without having pyproj installed, which is a bit more of a pain in the but than I want for a hard dependency right now. --- metpy/calc/kinematics.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/metpy/calc/kinematics.py b/metpy/calc/kinematics.py index 5c11491233c..9b51cc0634f 100644 --- a/metpy/calc/kinematics.py +++ b/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.') From 888d4294621a5dc14247404b31b4fb9dcb9ee158 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Fri, 11 May 2018 16:46:58 -0600 Subject: [PATCH 09/10] MNT: Add xarray as a hard dependency --- .travis.yml | 2 +- setup.py | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.travis.yml b/.travis.yml index d085e07f958..130cb042c5e 100644 --- a/.travis.yml +++ b/.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 diff --git a/setup.py b/setup.py index 3c3197e6b36..774b24aa248 100644 --- a/setup.py +++ b/setup.py @@ -51,15 +51,15 @@ python_requires='>=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*', install_requires=['matplotlib>=1.4', 'numpy>=1.10.0', 'scipy>=0.14', - 'pint>=0.8', 'enum34;python_version<"3.4"'], + 'pint>=0.8', 'xarray>=0.9.0', 'enum34;python_version<"3.4"'], extras_require={ 'cdm': ['pyproj>=1.9.4'], 'dev': ['ipython[all]>=3.1'], 'doc': ['sphinx>=1.4', 'sphinx-gallery', 'doc8', 'recommonmark', - 'netCDF4!=1.4', 'pandas'], - 'examples': ['cartopy>=0.13.1', 'pandas', 'xarray'], - 'test': ['pandas', 'pytest>=2.4', 'pytest-runner', 'pytest-mpl', 'pytest-flake8', - 'cartopy>=0.13.1', 'xarray', 'flake8>3.2.0', 'flake8-builtins!=1.4.0', + 'netCDF4!=1.4'], + 'examples': ['cartopy>=0.13.1'], + 'test': ['pytest>=2.4', 'pytest-runner', 'pytest-mpl', 'pytest-flake8', + 'cartopy>=0.13.1', 'flake8>3.2.0', 'flake8-builtins!=1.4.0', 'flake8-comprehensions', 'flake8-copyright', 'flake8-docstrings', 'flake8-import-order', 'flake8-mutable', 'flake8-pep3101', 'flake8-print', 'flake8-quotes', From d9d7f006abb9e362a8976c8b02e84b31ecc562d5 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Mon, 14 May 2018 16:17:56 -0600 Subject: [PATCH 10/10] MNT: Update codecov config Looks like we were missing some test lines --- .codecov.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.codecov.yml b/.codecov.yml index f252536cbe1..f3e2fb2ef3d 100644 --- a/.codecov.yml +++ b/.codecov.yml @@ -22,7 +22,9 @@ coverage: tests: target: 100% - paths: "metpy/.*/tests/.*" + paths: + - "metpy/.*/tests/.*" + - "metpy/tests/.*" notify: gitter: