From 3610c26f29f2c83893113c0a6154233dafe1f36a Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 8 Jan 2020 21:05:26 -0600 Subject: [PATCH 01/11] Switch to pyproj for projection to CF NetCDF grid mapping --- satpy/tests/writer_tests/test_cf.py | 4 ++-- satpy/writers/cf_writer.py | 14 +++++++++++--- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/satpy/tests/writer_tests/test_cf.py b/satpy/tests/writer_tests/test_cf.py index cab2981b26..345050b691 100644 --- a/satpy/tests/writer_tests/test_cf.py +++ b/satpy/tests/writer_tests/test_cf.py @@ -780,7 +780,7 @@ def _gm_matches(gmapping, expected): ds.attrs['area'] = geos res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'geos') + self.assertEqual(res.attrs['grid_mapping'], 'geostationary') _gm_matches(grid_mapping, geos_expected) # b) Projection does not have a corresponding CF representation (COSMO) @@ -834,7 +834,7 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = tmerc res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'tmerc') + self.assertEqual(res.attrs['grid_mapping'], 'transverse_mercator') _gm_matches(grid_mapping, tmerc_expected) # d) Projection that has a representation but no explicit a/b diff --git a/satpy/writers/cf_writer.py b/satpy/writers/cf_writer.py index e5bfe555da..1542b21ad5 100644 --- a/satpy/writers/cf_writer.py +++ b/satpy/writers/cf_writer.py @@ -215,8 +215,15 @@ def laea2cf(area): def create_grid_mapping(area): """Create the grid mapping instance for `area`.""" try: - grid_mapping = mappings[area.proj_dict['proj']](area) - grid_mapping['name'] = area.proj_dict['proj'] + try: + # let pyproj do the heavily lifting + # pyproj 2.0+ required + grid_mapping = area.crs.to_cf() + # not sure there is a better "standard" way of doing this + grid_mapping['name'] = grid_mapping['grid_mapping_name'] + except AttributeError: + grid_mapping = mappings[area.proj_dict['proj']](area) + grid_mapping['name'] = area.proj_dict['proj'] except KeyError: warnings.warn('The projection "{}" is either not CF compliant or not implemented yet. ' 'Using the proj4 string instead.'.format(area.proj_str)) @@ -258,7 +265,8 @@ def area2gridmapping(dataarray): """Convert an area to at CF grid mapping.""" area = dataarray.attrs['area'] attrs = create_grid_mapping(area) - if attrs is not None and 'name' in attrs.keys() and attrs['name'] != "proj4": + if (attrs is not None and 'name' in attrs and + attrs['name'] not in ("proj4", "unknown")): dataarray.attrs['grid_mapping'] = attrs['name'] name = attrs['name'] else: From 3a897a825b87c332818188475220c18d5b43a950 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 8 Jan 2020 21:09:18 -0600 Subject: [PATCH 02/11] Add cf test for projection with no a/b --- satpy/tests/writer_tests/test_cf.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/satpy/tests/writer_tests/test_cf.py b/satpy/tests/writer_tests/test_cf.py index 345050b691..1cb18bccc7 100644 --- a/satpy/tests/writer_tests/test_cf.py +++ b/satpy/tests/writer_tests/test_cf.py @@ -897,6 +897,31 @@ def _gm_matches(gmapping, expected): self.assertEqual(res.attrs['grid_mapping'], 'omerc') _gm_matches(grid_mapping, omerc_expected) + # d) Projection that has a representation but no explicit a/b + h = 35785831. + geos = pyresample.geometry.AreaDefinition( + area_id='geos', + description='geos', + proj_id='geos', + projection={'proj': 'geos', 'h': h, 'datum': 'WGS84', 'ellps': 'GRS80'}, + width=2, height=2, + area_extent=[-1, -1, 1, 1]) + geos_expected = xr.DataArray(data=0, + attrs={'perspective_point_height': h, + 'latitude_of_projection_origin': None, + 'longitude_of_projection_origin': None, + 'grid_mapping_name': 'geostationary', + 'reference_ellipsoid_name': 'GS80', + 'sweep_axis': None, + 'name': 'geos'}) + + ds = ds_base.copy() + ds.attrs['area'] = geos + res, grid_mapping = area2gridmapping(ds) + + self.assertEqual(res.attrs['grid_mapping'], 'geostationary') + self.assertEqual(grid_mapping, geos_expected) + def test_area2lonlat(self): """Test the conversion from areas to lon/lat.""" import pyresample.geometry From 5d4d265784fd1a6b518721693ab558e7dc0128ad Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Thu, 9 Jan 2020 15:43:49 +0100 Subject: [PATCH 03/11] Add cf test for omerc grid mapping --- satpy/tests/writer_tests/test_cf.py | 38 +++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/satpy/tests/writer_tests/test_cf.py b/satpy/tests/writer_tests/test_cf.py index 1cb18bccc7..2a2c8c5201 100644 --- a/satpy/tests/writer_tests/test_cf.py +++ b/satpy/tests/writer_tests/test_cf.py @@ -922,6 +922,44 @@ def _gm_matches(gmapping, expected): self.assertEqual(res.attrs['grid_mapping'], 'geostationary') self.assertEqual(grid_mapping, geos_expected) + # e) oblique Mercator + area = pyresample.geometry.AreaDefinition( + area_id='omerc_otf', + description='On-the-fly omerc area', + proj_id='omerc', + projection={'alpha': '9.02638777018478', 'ellps': 'WGS84', 'gamma': '0', 'k': '1', + 'lat_0': '-0.256794486098476', 'lonc': '13.7888658224205', + 'proj': 'omerc', 'units': 'm'}, + width=2837, + height=5940, + area_extent=[-1460463.0893, 3455291.3877, 1538407.1158, 9615788.8787] + ) + + omerc_expected = xr.DataArray(data=0, + attrs={}) + + omerc_dict = {'name': 'omerc', + 'azimuth_of_central_line': 9.02638777018478, + 'false_easting': 0., + 'false_northing': 0., + 'gamma': 0, + 'geographic_crs_name': "unknown", + 'grid_mapping_name': "oblique_mercator", + 'horizontal_datum_name': "unknown", + 'latitude_of_projection_origin': -0.256794486098476, + 'long_name': "omerc", + 'longitude_of_projection_origin': 13.7888658224205, + 'prime_meridian_name': "Greenwich", + 'reference_ellipsoid_name': "WGS84"} + omerc_expected = xr.DataArray(data=0, attrs=omerc_dict) + + ds = ds_base.copy() + ds.attrs['area'] = area + res, grid_mapping = area2gridmapping(ds) + + self.assertEqual(res.attrs['grid_mapping'], 'oblique_mercator') + assert(set(omerc_expected.attrs.items()) <= set(grid_mapping.attrs.items())) + def test_area2lonlat(self): """Test the conversion from areas to lon/lat.""" import pyresample.geometry From 10337c1c9b1c4924f0a1b2f9b0b7a2745b72df5f Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 23 Jan 2020 11:29:14 -0600 Subject: [PATCH 04/11] Update CF conversion to use pyproj 2.0+ --- .travis.yml | 2 +- appveyor.yml | 2 +- satpy/tests/writer_tests/test_cf.py | 74 ++++++++++++----------------- 3 files changed, 33 insertions(+), 45 deletions(-) diff --git a/.travis.yml b/.travis.yml index 485de13f4b..976a3acc10 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,7 +5,7 @@ env: - PYTHON_VERSION=$TRAVIS_PYTHON_VERSION - NUMPY_VERSION=1.17 - MAIN_CMD='python setup.py' - - CONDA_DEPENDENCIES='xarray!=0.13.0 dask distributed toolz Cython sphinx cartopy pillow matplotlib scipy pyyaml pyproj=2.4.1 pyresample coveralls coverage codecov behave netcdf4 h5py h5netcdf gdal rasterio imageio pyhdf mock libtiff geoviews zarr six python-eccodes' + - CONDA_DEPENDENCIES='xarray!=0.13.0 dask distributed toolz Cython sphinx cartopy pillow matplotlib scipy pyyaml pyproj=2.4.1 pyresample coveralls coverage codecov behave netcdf4 h5py h5netcdf gdal rasterio imageio pyhdf libtiff geoviews zarr six python-eccodes' - PIP_DEPENDENCIES='trollsift trollimage pyspectral pyorbital libtiff' - SETUP_XVFB=False - EVENT_TYPE='push pull_request' diff --git a/appveyor.yml b/appveyor.yml index 9fca571320..6dc0b7ee71 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -3,7 +3,7 @@ environment: PYTHON: "C:\\conda" MINICONDA_VERSION: "latest" CMD_IN_ENV: "cmd /E:ON /V:ON /C .\\ci-helpers\\appveyor\\windows_sdk.cmd" - CONDA_DEPENDENCIES: "xarray dask distributed toolz Cython sphinx cartopy pillow matplotlib scipy pyyaml pyproj pyresample coverage netcdf4 h5py h5netcdf gdal rasterio imageio pyhdf mock libtiff zarr six" + CONDA_DEPENDENCIES: "xarray dask distributed toolz Cython sphinx cartopy pillow matplotlib scipy pyyaml pyproj pyresample coverage netcdf4 h5py h5netcdf gdal rasterio imageio pyhdf libtiff zarr six" PIP_DEPENDENCIES: "trollsift trollimage pyspectral pyorbital libtiff" CONDA_CHANNELS: "conda-forge" diff --git a/satpy/tests/writer_tests/test_cf.py b/satpy/tests/writer_tests/test_cf.py index 2a2c8c5201..a0cad8e6f1 100644 --- a/satpy/tests/writer_tests/test_cf.py +++ b/satpy/tests/writer_tests/test_cf.py @@ -23,18 +23,11 @@ from datetime import datetime import tempfile from satpy import DatasetID +import unittest +from unittest import mock import numpy as np -if sys.version_info < (2, 7): - import unittest2 as unittest -else: - import unittest - -try: - from unittest import mock -except ImportError: - import mock try: from pyproj import CRS @@ -763,7 +756,8 @@ def _gm_matches(gmapping, expected): area_id='geos', description='geos', proj_id='geos', - projection={'proj': 'geos', 'h': h, 'a': a, 'b': b}, + projection={'proj': 'geos', 'h': h, 'a': a, 'b': b, + 'lat_0': 0, 'lon_0': 0}, width=2, height=2, area_extent=[-1, -1, 1, 1]) geos_expected = xr.DataArray(data=0, @@ -773,8 +767,8 @@ def _gm_matches(gmapping, expected): 'grid_mapping_name': 'geostationary', 'semi_major_axis': a, 'semi_minor_axis': b, - 'sweep_axis': None, - 'name': 'geos'}) + # 'sweep_angle_axis': None, + 'name': 'geostationary'}) ds = ds_base.copy() ds.attrs['area'] = geos @@ -808,7 +802,7 @@ def _gm_matches(gmapping, expected): self.assertEqual(proj_dict['proj'], 'ob_tran') self.assertEqual(proj_dict['o_proj'], 'stere') self.assertEqual(proj_dict['ellps'], 'WGS84') - self.assertEqual(grid_mapping.attrs['name'], 'proj4') + self.assertEqual(grid_mapping.attrs['name'], 'unknown') # c) Projection Transverse Mercator lat_0 = 36.5 @@ -829,7 +823,7 @@ def _gm_matches(gmapping, expected): 'reference_ellipsoid_name': 'WGS84', 'false_easting': 0., 'false_northing': 0., - 'name': 'tmerc'}) + 'name': 'transverse_mercator'}) ds = ds_base.copy() ds.attrs['area'] = tmerc @@ -843,7 +837,8 @@ def _gm_matches(gmapping, expected): area_id='geos', description='geos', proj_id='geos', - projection={'proj': 'geos', 'h': h, 'datum': 'WGS84', 'ellps': 'GRS80'}, + projection={'proj': 'geos', 'h': h, 'datum': 'WGS84', 'ellps': 'GRS80', + 'lat_0': 0, 'lon_0': 0}, width=2, height=2, area_extent=[-1, -1, 1, 1]) geos_expected = xr.DataArray(data=0, @@ -851,16 +846,16 @@ def _gm_matches(gmapping, expected): 'latitude_of_projection_origin': 0, 'longitude_of_projection_origin': 0, 'grid_mapping_name': 'geostationary', - 'semi_major_axis': 6378137.0, - 'semi_minor_axis': 6356752.314, - 'sweep_axis': None, - 'name': 'geos'}) + # 'semi_major_axis': 6378137.0, + # 'semi_minor_axis': 6356752.314, + # 'sweep_angle_axis': None, + 'name': 'geostationary'}) ds = ds_base.copy() ds.attrs['area'] = geos res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'geos') + self.assertEqual(res.attrs['grid_mapping'], 'geostationary') _gm_matches(grid_mapping, geos_expected) # e) oblique Mercator @@ -876,17 +871,15 @@ def _gm_matches(gmapping, expected): area_extent=[-1460463.0893, 3455291.3877, 1538407.1158, 9615788.8787] ) - omerc_dict = {'name': 'omerc', + omerc_dict = {'name': 'oblique_mercator', 'azimuth_of_central_line': 9.02638777018478, 'false_easting': 0., 'false_northing': 0., - 'gamma': 0, - 'geographic_crs_name': "unknown", + # 'gamma': 0, # this is not CF compliant 'grid_mapping_name': "oblique_mercator", - 'horizontal_datum_name': "unknown", 'latitude_of_projection_origin': -0.256794486098476, 'longitude_of_projection_origin': 13.7888658224205, - 'prime_meridian_name': "Greenwich", + # 'prime_meridian_name': "Greenwich", 'reference_ellipsoid_name': "WGS84"} omerc_expected = xr.DataArray(data=0, attrs=omerc_dict) @@ -894,7 +887,7 @@ def _gm_matches(gmapping, expected): ds.attrs['area'] = area res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'omerc') + self.assertEqual(res.attrs['grid_mapping'], 'oblique_mercator') _gm_matches(grid_mapping, omerc_expected) # d) Projection that has a representation but no explicit a/b @@ -903,24 +896,24 @@ def _gm_matches(gmapping, expected): area_id='geos', description='geos', proj_id='geos', - projection={'proj': 'geos', 'h': h, 'datum': 'WGS84', 'ellps': 'GRS80'}, + projection={'proj': 'geos', 'h': h, 'datum': 'WGS84', 'ellps': 'GRS80', + 'lat_0': 0, 'lon_0': 0}, width=2, height=2, area_extent=[-1, -1, 1, 1]) geos_expected = xr.DataArray(data=0, attrs={'perspective_point_height': h, - 'latitude_of_projection_origin': None, - 'longitude_of_projection_origin': None, + 'latitude_of_projection_origin': 0, + 'longitude_of_projection_origin': 0, 'grid_mapping_name': 'geostationary', - 'reference_ellipsoid_name': 'GS80', - 'sweep_axis': None, - 'name': 'geos'}) + 'reference_ellipsoid_name': 'GRS80', + 'name': 'geostationary'}) ds = ds_base.copy() ds.attrs['area'] = geos res, grid_mapping = area2gridmapping(ds) self.assertEqual(res.attrs['grid_mapping'], 'geostationary') - self.assertEqual(grid_mapping, geos_expected) + _gm_matches(grid_mapping, geos_expected) # e) oblique Mercator area = pyresample.geometry.AreaDefinition( @@ -935,21 +928,16 @@ def _gm_matches(gmapping, expected): area_extent=[-1460463.0893, 3455291.3877, 1538407.1158, 9615788.8787] ) - omerc_expected = xr.DataArray(data=0, - attrs={}) - - omerc_dict = {'name': 'omerc', + omerc_dict = {'name': 'oblique_mercator', 'azimuth_of_central_line': 9.02638777018478, 'false_easting': 0., 'false_northing': 0., - 'gamma': 0, - 'geographic_crs_name': "unknown", + # 'gamma': 0, # not CF compliant 'grid_mapping_name': "oblique_mercator", - 'horizontal_datum_name': "unknown", 'latitude_of_projection_origin': -0.256794486098476, - 'long_name': "omerc", + # 'long_name': "omerc", 'longitude_of_projection_origin': 13.7888658224205, - 'prime_meridian_name': "Greenwich", + # 'prime_meridian_name': "Greenwich", 'reference_ellipsoid_name': "WGS84"} omerc_expected = xr.DataArray(data=0, attrs=omerc_dict) @@ -958,7 +946,7 @@ def _gm_matches(gmapping, expected): res, grid_mapping = area2gridmapping(ds) self.assertEqual(res.attrs['grid_mapping'], 'oblique_mercator') - assert(set(omerc_expected.attrs.items()) <= set(grid_mapping.attrs.items())) + _gm_matches(grid_mapping, omerc_expected) def test_area2lonlat(self): """Test the conversion from areas to lon/lat.""" From f92fbee46e7d19207b9069672b4062d2df5bd62f Mon Sep 17 00:00:00 2001 From: David Hoese Date: Sun, 29 Mar 2020 21:13:54 -0500 Subject: [PATCH 05/11] Switch to only using pyproj for CF CRS conversion --- satpy/tests/writer_tests/test_cf.py | 59 +++++++------- satpy/writers/cf_writer.py | 116 ++++------------------------ 2 files changed, 44 insertions(+), 131 deletions(-) diff --git a/satpy/tests/writer_tests/test_cf.py b/satpy/tests/writer_tests/test_cf.py index cdb0ebc1f0..78d1574b42 100644 --- a/satpy/tests/writer_tests/test_cf.py +++ b/satpy/tests/writer_tests/test_cf.py @@ -761,13 +761,13 @@ def _gm_matches(gmapping, expected): 'semi_major_axis': a, 'semi_minor_axis': b, # 'sweep_angle_axis': None, - 'name': 'geostationary'}) + }) ds = ds_base.copy() ds.attrs['area'] = geos res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'geostationary') + self.assertEqual(res.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) # b) Projection does not have a corresponding CF representation (COSMO) @@ -784,18 +784,19 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = cosmo7 - with mock.patch('satpy.writers.cf_writer.warnings.warn') as warn: - res, grid_mapping = area2gridmapping(ds) - warn.assert_called() - proj_dict = pyresample.geometry.proj4_str_to_dict(res.attrs['grid_proj4']) - self.assertEqual(proj_dict['lon_0'], 4.535) - self.assertEqual(proj_dict['lat_0'], 46.0) - self.assertEqual(proj_dict['o_lon_p'], -5.465) - self.assertEqual(proj_dict['o_lat_p'], 90.0) - self.assertEqual(proj_dict['proj'], 'ob_tran') - self.assertEqual(proj_dict['o_proj'], 'stere') - self.assertEqual(proj_dict['ellps'], 'WGS84') - self.assertEqual(grid_mapping.attrs['name'], 'unknown') + # with mock.patch('satpy.writers.cf_writer.warnings.warn') as warn: + res, grid_mapping = area2gridmapping(ds) + # warn.assert_called() + # proj_dict = pyresample.geometry.proj4_str_to_dict(res.attrs['grid_proj4']) + # self.assertEqual(proj_dict['lon_0'], 4.535) + # self.assertEqual(proj_dict['lat_0'], 46.0) + # self.assertEqual(proj_dict['o_lon_p'], -5.465) + # self.assertEqual(proj_dict['o_lat_p'], 90.0) + # self.assertEqual(proj_dict['proj'], 'ob_tran') + # self.assertEqual(proj_dict['o_proj'], 'stere') + # self.assertEqual(proj_dict['ellps'], 'WGS84') + self.assertIn('crs_wkt', grid_mapping.attrs) + self.assertEqual(res.attrs['grid_mapping'], 'cosmo7') # c) Projection Transverse Mercator lat_0 = 36.5 @@ -813,15 +814,15 @@ def _gm_matches(gmapping, expected): attrs={'latitude_of_projection_origin': lat_0, 'longitude_of_central_meridian': lon_0, 'grid_mapping_name': 'transverse_mercator', - 'reference_ellipsoid_name': 'WGS84', + 'reference_ellipsoid_name': 'WGS 84', 'false_easting': 0., 'false_northing': 0., - 'name': 'transverse_mercator'}) + }) ds = ds_base.copy() ds.attrs['area'] = tmerc res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'transverse_mercator') + self.assertEqual(res.attrs['grid_mapping'], 'tmerc') _gm_matches(grid_mapping, tmerc_expected) # d) Projection that has a representation but no explicit a/b @@ -842,13 +843,13 @@ def _gm_matches(gmapping, expected): # 'semi_major_axis': 6378137.0, # 'semi_minor_axis': 6356752.314, # 'sweep_angle_axis': None, - 'name': 'geostationary'}) + }) ds = ds_base.copy() ds.attrs['area'] = geos res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'geostationary') + self.assertEqual(res.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) # e) oblique Mercator @@ -864,8 +865,7 @@ def _gm_matches(gmapping, expected): area_extent=[-1460463.0893, 3455291.3877, 1538407.1158, 9615788.8787] ) - omerc_dict = {'name': 'oblique_mercator', - 'azimuth_of_central_line': 9.02638777018478, + omerc_dict = {'azimuth_of_central_line': 9.02638777018478, 'false_easting': 0., 'false_northing': 0., # 'gamma': 0, # this is not CF compliant @@ -873,14 +873,14 @@ def _gm_matches(gmapping, expected): 'latitude_of_projection_origin': -0.256794486098476, 'longitude_of_projection_origin': 13.7888658224205, # 'prime_meridian_name': "Greenwich", - 'reference_ellipsoid_name': "WGS84"} + 'reference_ellipsoid_name': "WGS 84"} omerc_expected = xr.DataArray(data=0, attrs=omerc_dict) ds = ds_base.copy() ds.attrs['area'] = area res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'oblique_mercator') + self.assertEqual(res.attrs['grid_mapping'], 'omerc_otf') _gm_matches(grid_mapping, omerc_expected) # d) Projection that has a representation but no explicit a/b @@ -898,14 +898,14 @@ def _gm_matches(gmapping, expected): 'latitude_of_projection_origin': 0, 'longitude_of_projection_origin': 0, 'grid_mapping_name': 'geostationary', - 'reference_ellipsoid_name': 'GRS80', - 'name': 'geostationary'}) + 'reference_ellipsoid_name': 'WGS 84', + }) ds = ds_base.copy() ds.attrs['area'] = geos res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'geostationary') + self.assertEqual(res.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) # e) oblique Mercator @@ -921,8 +921,7 @@ def _gm_matches(gmapping, expected): area_extent=[-1460463.0893, 3455291.3877, 1538407.1158, 9615788.8787] ) - omerc_dict = {'name': 'oblique_mercator', - 'azimuth_of_central_line': 9.02638777018478, + omerc_dict = {'azimuth_of_central_line': 9.02638777018478, 'false_easting': 0., 'false_northing': 0., # 'gamma': 0, # not CF compliant @@ -931,14 +930,14 @@ def _gm_matches(gmapping, expected): # 'long_name': "omerc", 'longitude_of_projection_origin': 13.7888658224205, # 'prime_meridian_name': "Greenwich", - 'reference_ellipsoid_name': "WGS84"} + 'reference_ellipsoid_name': "WGS 84"} omerc_expected = xr.DataArray(data=0, attrs=omerc_dict) ds = ds_base.copy() ds.attrs['area'] = area res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'oblique_mercator') + self.assertEqual(res.attrs['grid_mapping'], 'omerc_otf') _gm_matches(grid_mapping, omerc_expected) def test_area2lonlat(self): diff --git a/satpy/writers/cf_writer.py b/satpy/writers/cf_writer.py index 1542b21ad5..50ed1fead6 100644 --- a/satpy/writers/cf_writer.py +++ b/satpy/writers/cf_writer.py @@ -113,6 +113,12 @@ from satpy.writers import Writer from satpy.writers.utils import flatten_dict +from distutils.version import LooseVersion +import pyproj +if LooseVersion(pyproj.__version__) < LooseVersion('2.4.1'): + # technically 2.2, but important bug fixes in 2.4.1 + raise ImportError("'cf' writer requires pyproj 2.4.1 or greater") + logger = logging.getLogger(__name__) @@ -139,97 +145,12 @@ CF_VERSION = 'CF-1.7' -def tmerc2cf(area): - """Return the cf grid mapping for the tmerc projection.""" - proj_dict = area.proj_dict - args = dict(latitude_of_projection_origin=proj_dict.get('lat_0'), - longitude_of_central_meridian=proj_dict.get('lon_0'), - grid_mapping_name='transverse_mercator', - reference_ellipsoid_name=proj_dict.get('ellps', 'WGS84'), - false_easting=0., - false_northing=0. - ) - if "no_rot" in proj_dict: - args['no_rotation'] = 1 - if "gamma" in proj_dict: - args['gamma'] = proj_dict['gamma'] - return args - - -def omerc2cf(area): - """Return the cf grid mapping for the omerc projection.""" - proj_dict = area.proj_dict - - args = dict(azimuth_of_central_line=proj_dict.get('alpha'), - latitude_of_projection_origin=proj_dict.get('lat_0'), - longitude_of_projection_origin=proj_dict.get('lonc'), - grid_mapping_name='oblique_mercator', - reference_ellipsoid_name=proj_dict.get('ellps', 'WGS84'), - prime_meridian_name=proj_dict.get('pm', 'Greenwich'), - horizontal_datum_name=proj_dict.get('datum', 'unknown'), - geographic_crs_name='unknown', - false_easting=0., - false_northing=0. - ) - if "no_rot" in proj_dict: - args['no_rotation'] = 1 - if "gamma" in proj_dict: - args['gamma'] = proj_dict['gamma'] - return args - - -def geos2cf(area): - """Return the cf grid mapping for the geos projection.""" - from pyresample.utils import proj4_radius_parameters - proj_dict = area.proj_dict - a, b = proj4_radius_parameters(proj_dict) - args = dict(perspective_point_height=proj_dict.get('h'), - latitude_of_projection_origin=proj_dict.get('lat_0', 0), - longitude_of_projection_origin=proj_dict.get('lon_0', 0), - grid_mapping_name='geostationary', - semi_major_axis=a, - semi_minor_axis=b, - # semi_major_axis=proj_dict.get('a'), - # semi_minor_axis=proj_dict.get('b'), - sweep_axis=proj_dict.get('sweep'), - ) - return args - - -def laea2cf(area): - """Return the cf grid mapping for the laea projection.""" - proj_dict = area.proj_dict - args = dict(latitude_of_projection_origin=proj_dict.get('lat_0'), - longitude_of_projection_origin=proj_dict.get('lon_0'), - grid_mapping_name='lambert_azimuthal_equal_area', - ) - return args - - -mappings = {'omerc': omerc2cf, - 'laea': laea2cf, - 'geos': geos2cf, - 'tmerc': tmerc2cf} - - def create_grid_mapping(area): """Create the grid mapping instance for `area`.""" - try: - try: - # let pyproj do the heavily lifting - # pyproj 2.0+ required - grid_mapping = area.crs.to_cf() - # not sure there is a better "standard" way of doing this - grid_mapping['name'] = grid_mapping['grid_mapping_name'] - except AttributeError: - grid_mapping = mappings[area.proj_dict['proj']](area) - grid_mapping['name'] = area.proj_dict['proj'] - except KeyError: - warnings.warn('The projection "{}" is either not CF compliant or not implemented yet. ' - 'Using the proj4 string instead.'.format(area.proj_str)) - grid_mapping = {'name': 'proj4', 'proj4': area.proj_str} - - return grid_mapping + # let pyproj do the heavily lifting + # pyproj 2.0+ required + grid_mapping = area.crs.to_cf() + return area.area_id, grid_mapping def get_extra_ds(dataset): @@ -245,7 +166,8 @@ def get_extra_ds(dataset): def area2lonlat(dataarray): """Convert an area to longitudes and latitudes.""" area = dataarray.attrs['area'] - lons, lats = area.get_lonlats_dask() + chunks = getattr(dataarray.data, 'chunks', None) + lons, lats = area.get_lonlats(chunks=chunks) lons = xr.DataArray(lons, dims=['y', 'x'], attrs={'name': "longitude", 'standard_name': "longitude", @@ -264,17 +186,9 @@ def area2lonlat(dataarray): def area2gridmapping(dataarray): """Convert an area to at CF grid mapping.""" area = dataarray.attrs['area'] - attrs = create_grid_mapping(area) - if (attrs is not None and 'name' in attrs and - attrs['name'] not in ("proj4", "unknown")): - dataarray.attrs['grid_mapping'] = attrs['name'] - name = attrs['name'] - else: - # Handle the case when the projection cannot be converted to a standard CF representation or this has not - # been implemented yet. - dataarray.attrs['grid_proj4'] = area.proj4_string - name = "proj4" - return [dataarray, xr.DataArray(0, attrs=attrs, name=name)] + gmapping_var_name, attrs = create_grid_mapping(area) + dataarray.attrs['grid_mapping'] = gmapping_var_name + return [dataarray, xr.DataArray(0, attrs=attrs, name=gmapping_var_name)] def area2cf(dataarray, strict=False): From 1078deae53558cc10ce4d58b03230bb1cff081e5 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Mon, 30 Mar 2020 08:09:56 -0500 Subject: [PATCH 06/11] Remove redundant CF omerc test --- satpy/tests/writer_tests/test_cf.py | 34 +---------------------------- 1 file changed, 1 insertion(+), 33 deletions(-) diff --git a/satpy/tests/writer_tests/test_cf.py b/satpy/tests/writer_tests/test_cf.py index 78d1574b42..b3a85ae89d 100644 --- a/satpy/tests/writer_tests/test_cf.py +++ b/satpy/tests/writer_tests/test_cf.py @@ -883,7 +883,7 @@ def _gm_matches(gmapping, expected): self.assertEqual(res.attrs['grid_mapping'], 'omerc_otf') _gm_matches(grid_mapping, omerc_expected) - # d) Projection that has a representation but no explicit a/b + # f) Projection that has a representation but no explicit a/b h = 35785831. geos = pyresample.geometry.AreaDefinition( area_id='geos', @@ -908,38 +908,6 @@ def _gm_matches(gmapping, expected): self.assertEqual(res.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) - # e) oblique Mercator - area = pyresample.geometry.AreaDefinition( - area_id='omerc_otf', - description='On-the-fly omerc area', - proj_id='omerc', - projection={'alpha': '9.02638777018478', 'ellps': 'WGS84', 'gamma': '0', 'k': '1', - 'lat_0': '-0.256794486098476', 'lonc': '13.7888658224205', - 'proj': 'omerc', 'units': 'm'}, - width=2837, - height=5940, - area_extent=[-1460463.0893, 3455291.3877, 1538407.1158, 9615788.8787] - ) - - omerc_dict = {'azimuth_of_central_line': 9.02638777018478, - 'false_easting': 0., - 'false_northing': 0., - # 'gamma': 0, # not CF compliant - 'grid_mapping_name': "oblique_mercator", - 'latitude_of_projection_origin': -0.256794486098476, - # 'long_name': "omerc", - 'longitude_of_projection_origin': 13.7888658224205, - # 'prime_meridian_name': "Greenwich", - 'reference_ellipsoid_name': "WGS 84"} - omerc_expected = xr.DataArray(data=0, attrs=omerc_dict) - - ds = ds_base.copy() - ds.attrs['area'] = area - res, grid_mapping = area2gridmapping(ds) - - self.assertEqual(res.attrs['grid_mapping'], 'omerc_otf') - _gm_matches(grid_mapping, omerc_expected) - def test_area2lonlat(self): """Test the conversion from areas to lon/lat.""" import pyresample.geometry From 1218d219e74b2654ce6d180ca7e2114be2c042ad Mon Sep 17 00:00:00 2001 From: David Hoese Date: Mon, 30 Mar 2020 10:45:24 -0500 Subject: [PATCH 07/11] Remove pyproj restriction on travis environments --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 3efc2e8f88..88a9f6889d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -6,7 +6,7 @@ env: - NUMPY_VERSION=stable - MAIN_CMD='python setup.py' - CONDA_DEPENDENCIES='xarray dask distributed toolz Cython sphinx cartopy pillow matplotlib scipy pyyaml - pyproj=2.4.1 pyresample coveralls coverage codecov behave netcdf4 h5py h5netcdf gdal rasterio imageio pyhdf + pyproj pyresample coveralls coverage codecov behave netcdf4 h5py h5netcdf gdal rasterio imageio pyhdf mock libtiff geoviews zarr python-eccodes geoviews pytest pytest-cov' - PIP_DEPENDENCIES='trollsift trollimage pyspectral pyorbital libtiff' - SETUP_XVFB=False From be98fe2ab5b82ca4e585bb17ee0470f794015a95 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Mon, 30 Mar 2020 12:30:09 -0500 Subject: [PATCH 08/11] Add travis workaround for long logs See https://github.com/travis-ci/travis-ci/issues/8920 --- .travis.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.travis.yml b/.travis.yml index 88a9f6889d..10dcf14ea2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -41,6 +41,8 @@ matrix: install: - git clone --depth 1 git://github.com/astropy/ci-helpers.git - source ci-helpers/travis/setup_conda.sh + # See https://github.com/travis-ci/travis-ci/issues/8920 + - python -c "import fcntl; fcntl.fcntl(1, fcntl.F_SETFL, 0)" - if [ "$UNSTABLE_DEPS" == "True" ]; then python -m pip install -f https://7933911d6844c6c53a7d-47bd50c35cd79bd838daf386af554a83.ssl.cf2.rackcdn.com From 26a93a484e4bb613406d3d6dd7a9f3ebf7136bdf Mon Sep 17 00:00:00 2001 From: David Hoese Date: Mon, 30 Mar 2020 15:39:41 -0500 Subject: [PATCH 09/11] Fix CF writer tests to be more thorough and remove extra name attr --- satpy/tests/writer_tests/test_cf.py | 119 +++++++++++++--------------- satpy/writers/cf_writer.py | 19 +++-- 2 files changed, 69 insertions(+), 69 deletions(-) diff --git a/satpy/tests/writer_tests/test_cf.py b/satpy/tests/writer_tests/test_cf.py index b3a85ae89d..4b34cf1d8f 100644 --- a/satpy/tests/writer_tests/test_cf.py +++ b/satpy/tests/writer_tests/test_cf.py @@ -534,26 +534,18 @@ def test_da2cf(self): self.assertDictWithArraysEqual(res_flat.attrs, attrs_expected_flat) @mock.patch('satpy.writers.cf_writer.CFWriter.__init__', return_value=None) - @mock.patch('satpy.writers.cf_writer.area2cf') - @mock.patch('satpy.writers.cf_writer.CFWriter.da2cf') - @mock.patch('satpy.writers.cf_writer.make_alt_coords_unique') - @mock.patch('satpy.writers.cf_writer.assert_xy_unique') - @mock.patch('satpy.writers.cf_writer.link_coords') - def test_collect_datasets(self, link_coords, assert_xy_unique, make_alt_coords_unique, da2cf, area2cf, *mocks): + def test_collect_datasets(self, *mocks): """Test collecting CF datasets from a DataArray objects.""" from satpy.writers.cf_writer import CFWriter import xarray as xr - - # Patch methods - def identity(arg, **kwargs): - return arg - - def raise_key_error(arg, **kwargs): - raise KeyError - - da2cf.side_effect = identity - area2cf.side_effect = raise_key_error - make_alt_coords_unique.return_value = 'unique_coords' + import pyresample.geometry + geos = pyresample.geometry.AreaDefinition( + area_id='geos', + description='geos', + proj_id='geos', + projection={'proj': 'geos', 'h': 35785831., 'a': 6378169., 'b': 6356583.8}, + width=2, height=2, + area_extent=[-1, -1, 1, 1]) # Define test datasets data = [[1, 2], [3, 4]] @@ -563,32 +555,29 @@ def raise_key_error(arg, **kwargs): tstart = datetime(2019, 4, 1, 12, 0) tend = datetime(2019, 4, 1, 12, 15) datasets = [xr.DataArray(data=data, dims=('y', 'x'), coords={'y': y, 'x': x, 'acq_time': ('y', time)}, - attrs={'name': 'var1', 'start_time': tstart, 'end_time': tend}), + attrs={'name': 'var1', 'start_time': tstart, 'end_time': tend, 'area': geos}), xr.DataArray(data=data, dims=('y', 'x'), coords={'y': y, 'x': x, 'acq_time': ('y', time)}, - attrs={'name': 'var2'})] - expected = {'var1': datasets[0], 'var2': datasets[1]} + attrs={'name': 'var2', 'long_name': 'variable 2'})] # Collect datasets writer = CFWriter() datas, start_times, end_times = writer._collect_datasets(datasets, include_lonlats=True) # Test results - self.assertEqual(datas, 'unique_coords') - self.assertListEqual(start_times, [tstart, None]) - self.assertListEqual(end_times, [tend, None]) - - # Test method calls - self.assertEqual(len(area2cf.call_args_list), 2) - for call_args, ds in zip(area2cf.call_args_list, datasets): - self.assertEqual(call_args, mock.call(ds, strict=True)) - - for func in (assert_xy_unique, link_coords, make_alt_coords_unique): - func.assert_called() - call_arg = func.call_args[0][0] - self.assertIsInstance(call_arg, dict) - self.assertSetEqual(set(call_arg.keys()), {'var1', 'var2'}) - for key, ds in expected.items(): - self.assertTrue(call_arg[key].identical(ds)) + self.assertEqual(len(datas), 3) + self.assertEqual(set(datas.keys()), {'var1', 'var2', 'geos'}) + self.assertListEqual(start_times, [None, tstart, None]) + self.assertListEqual(end_times, [None, tend, None]) + var1 = datas['var1'] + var2 = datas['var2'] + self.assertEqual(var1.name, 'var1') + self.assertEqual(var1.attrs['grid_mapping'], 'geos') + self.assertEqual(var1.attrs['start_time'], '2019-04-01 12:00:00') + self.assertEqual(var1.attrs['end_time'], '2019-04-01 12:15:00') + self.assertEqual(var1.attrs['long_name'], 'var1') + # variable 2 + self.assertNotIn('grid_mapping', var2.attrs) + self.assertEqual(var2.attrs['long_name'], 'variable 2') def test_assert_xy_unique(self): """Test that the x and y coordinates are unique.""" @@ -680,16 +669,12 @@ def test_make_alt_coords_unique(self): self.assertNotIn('var1_acq_time', res['var1'].coords) self.assertNotIn('var2_acq_time', res['var2'].coords) - @mock.patch('satpy.writers.cf_writer.area2lonlat') - @mock.patch('satpy.writers.cf_writer.area2gridmapping') - def test_area2cf(self, area2gridmapping, area2lonlat): + def test_area2cf(self): """Test the conversion of an area to CF standards.""" import xarray as xr import pyresample.geometry from satpy.writers.cf_writer import area2cf - area2gridmapping.side_effect = lambda x: [1, 2, 3] - area2lonlat.side_effect = lambda x: [4, 5, 6] ds_base = xr.DataArray(data=[[1, 2], [3, 4]], dims=('y', 'x'), coords={'y': [1, 2], 'x': [3, 4]}, attrs={'name': 'var1'}) @@ -705,13 +690,21 @@ def test_area2cf(self, area2gridmapping, area2lonlat): ds.attrs['area'] = geos res = area2cf(ds) - self.assertEqual(len(res), 4) - self.assertListEqual(res[0:3], [1, 2, 3]) - self.assertTrue(ds.identical(res[3])) + self.assertEqual(len(res), 2) + self.assertEqual(res[0].size, 1) # grid mapping variable + self.assertEqual(res[0].name, res[1].attrs['grid_mapping']) # b) Area Definition and strict=False - area2cf(ds, strict=True) - area2lonlat.assert_called() + ds = ds_base.copy(deep=True) + ds.attrs['area'] = geos + res = area2cf(ds, strict=True) + # same as above + self.assertEqual(len(res), 2) + self.assertEqual(res[0].size, 1) # grid mapping variable + self.assertEqual(res[0].name, res[1].attrs['grid_mapping']) + # but now also have the lon/lats + self.assertIn('longitude', res[1].coords) + self.assertIn('latitude', res[1].coords) # c) Swath Definition swath = pyresample.geometry.SwathDefinition(lons=[[1, 1], [2, 2]], lats=[[1, 2], [1, 2]]) @@ -719,9 +712,10 @@ def test_area2cf(self, area2gridmapping, area2lonlat): ds.attrs['area'] = swath res = area2cf(ds) - self.assertEqual(len(res), 4) - self.assertListEqual(res[0:3], [4, 5, 6]) - self.assertTrue(ds.identical(res[3])) + self.assertEqual(len(res), 1) + self.assertIn('longitude', res[0].coords) + self.assertIn('latitude', res[0].coords) + self.assertNotIn('grid_mapping', res[0].attrs) def test_area2gridmapping(self): """Test the conversion from pyresample area object to CF grid mapping.""" @@ -765,9 +759,9 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = geos - res, grid_mapping = area2gridmapping(ds) + grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'geos') + self.assertEqual(ds.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) # b) Projection does not have a corresponding CF representation (COSMO) @@ -785,7 +779,7 @@ def _gm_matches(gmapping, expected): ds.attrs['area'] = cosmo7 # with mock.patch('satpy.writers.cf_writer.warnings.warn') as warn: - res, grid_mapping = area2gridmapping(ds) + grid_mapping = area2gridmapping(ds) # warn.assert_called() # proj_dict = pyresample.geometry.proj4_str_to_dict(res.attrs['grid_proj4']) # self.assertEqual(proj_dict['lon_0'], 4.535) @@ -796,7 +790,7 @@ def _gm_matches(gmapping, expected): # self.assertEqual(proj_dict['o_proj'], 'stere') # self.assertEqual(proj_dict['ellps'], 'WGS84') self.assertIn('crs_wkt', grid_mapping.attrs) - self.assertEqual(res.attrs['grid_mapping'], 'cosmo7') + self.assertEqual(ds.attrs['grid_mapping'], 'cosmo7') # c) Projection Transverse Mercator lat_0 = 36.5 @@ -821,8 +815,8 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = tmerc - res, grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'tmerc') + grid_mapping = area2gridmapping(ds) + self.assertEqual(ds.attrs['grid_mapping'], 'tmerc') _gm_matches(grid_mapping, tmerc_expected) # d) Projection that has a representation but no explicit a/b @@ -847,9 +841,9 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = geos - res, grid_mapping = area2gridmapping(ds) + grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'geos') + self.assertEqual(ds.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) # e) oblique Mercator @@ -878,9 +872,9 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = area - res, grid_mapping = area2gridmapping(ds) + grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'omerc_otf') + self.assertEqual(ds.attrs['grid_mapping'], 'omerc_otf') _gm_matches(grid_mapping, omerc_expected) # f) Projection that has a representation but no explicit a/b @@ -903,9 +897,9 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = geos - res, grid_mapping = area2gridmapping(ds) + grid_mapping = area2gridmapping(ds) - self.assertEqual(res.attrs['grid_mapping'], 'geos') + self.assertEqual(ds.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) def test_area2lonlat(self): @@ -925,7 +919,8 @@ def test_area2lonlat(self): lons_ref, lats_ref = area.get_lonlats() dataarray = xr.DataArray(data=[[1, 2], [3, 4]], dims=('y', 'x'), attrs={'area': area}) - res = area2lonlat(dataarray) + area2lonlat(dataarray) + res = [dataarray] # modified in-place self.assertEqual(len(res), 1) self.assertEqual(set(res[0].coords), {'longitude', 'latitude'}) diff --git a/satpy/writers/cf_writer.py b/satpy/writers/cf_writer.py index 50ed1fead6..0c659911df 100644 --- a/satpy/writers/cf_writer.py +++ b/satpy/writers/cf_writer.py @@ -180,7 +180,6 @@ def area2lonlat(dataarray): name='latitude') dataarray['longitude'] = lons dataarray['latitude'] = lats - return [dataarray] def area2gridmapping(dataarray): @@ -188,7 +187,7 @@ def area2gridmapping(dataarray): area = dataarray.attrs['area'] gmapping_var_name, attrs = create_grid_mapping(area) dataarray.attrs['grid_mapping'] = gmapping_var_name - return [dataarray, xr.DataArray(0, attrs=attrs, name=gmapping_var_name)] + return xr.DataArray(0, attrs=attrs, name=gmapping_var_name) def area2cf(dataarray, strict=False): @@ -196,9 +195,10 @@ def area2cf(dataarray, strict=False): res = [] dataarray = dataarray.copy(deep=True) if isinstance(dataarray.attrs['area'], SwathDefinition) or strict: - res = area2lonlat(dataarray) + # modifies dataarray in-place + area2lonlat(dataarray) if isinstance(dataarray.attrs['area'], AreaDefinition): - res.extend(area2gridmapping(dataarray)) + res.append(area2gridmapping(dataarray)) res.append(dataarray) return res @@ -418,6 +418,9 @@ def da2cf(dataarray, epoch=EPOCH, flatten_attrs=False, exclude_attrs=None, compr exclude_attrs = [] new_data = dataarray.copy() + if 'name' in new_data.attrs: + name = new_data.attrs.pop('name') + new_data = new_data.rename(name) # Remove area as well as user-defined attributes for key in ['area'] + exclude_attrs: @@ -456,7 +459,8 @@ def da2cf(dataarray, epoch=EPOCH, flatten_attrs=False, exclude_attrs=None, compr if 'crs' in new_data.coords: new_data = new_data.drop('crs') - new_data.attrs.setdefault('long_name', new_data.attrs.pop('name')) + if 'long_name' not in new_data.attrs and 'standard_name' not in new_data.attrs: + new_data.attrs['long_name'] = new_data.name if 'prerequisites' in new_data.attrs: new_data.attrs['prerequisites'] = [np.string_(str(prereq)) for prereq in new_data.attrs['prerequisites']] @@ -493,8 +497,9 @@ def _collect_datasets(self, datasets, epoch=EPOCH, flatten_attrs=False, exclude_ for new_ds in new_datasets: start_times.append(new_ds.attrs.get("start_time", None)) end_times.append(new_ds.attrs.get("end_time", None)) - datas[new_ds.attrs['name']] = self.da2cf(new_ds, epoch=epoch, flatten_attrs=flatten_attrs, - exclude_attrs=exclude_attrs, compression=compression) + new_var = self.da2cf(new_ds, epoch=epoch, flatten_attrs=flatten_attrs, + exclude_attrs=exclude_attrs, compression=compression) + datas[new_var.name] = new_var # Check and prepare coordinates assert_xy_unique(datas) From 641bc862fdbb498c9e1f5832f2e77ecd7e277d81 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 31 Mar 2020 09:29:34 -0500 Subject: [PATCH 10/11] Refactor area2X functions to copy instead of doing things in-place --- satpy/tests/writer_tests/test_cf.py | 57 +++++++++++++++-------------- satpy/writers/cf_writer.py | 12 +++--- 2 files changed, 36 insertions(+), 33 deletions(-) diff --git a/satpy/tests/writer_tests/test_cf.py b/satpy/tests/writer_tests/test_cf.py index 4b34cf1d8f..001d18714f 100644 --- a/satpy/tests/writer_tests/test_cf.py +++ b/satpy/tests/writer_tests/test_cf.py @@ -759,10 +759,15 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = geos - grid_mapping = area2gridmapping(ds) + new_ds, grid_mapping = area2gridmapping(ds) + if 'sweep_angle_axis' in grid_mapping.attrs: + # older versions of pyproj might not include this + self.assertEqual(grid_mapping.attrs['sweep_angle_axis'], 'y') - self.assertEqual(ds.attrs['grid_mapping'], 'geos') + self.assertEqual(new_ds.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) + # should not have been modified + self.assertNotIn('grid_mapping', ds.attrs) # b) Projection does not have a corresponding CF representation (COSMO) cosmo7 = pyresample.geometry.AreaDefinition( @@ -778,19 +783,15 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = cosmo7 - # with mock.patch('satpy.writers.cf_writer.warnings.warn') as warn: - grid_mapping = area2gridmapping(ds) - # warn.assert_called() - # proj_dict = pyresample.geometry.proj4_str_to_dict(res.attrs['grid_proj4']) - # self.assertEqual(proj_dict['lon_0'], 4.535) - # self.assertEqual(proj_dict['lat_0'], 46.0) - # self.assertEqual(proj_dict['o_lon_p'], -5.465) - # self.assertEqual(proj_dict['o_lat_p'], 90.0) - # self.assertEqual(proj_dict['proj'], 'ob_tran') - # self.assertEqual(proj_dict['o_proj'], 'stere') - # self.assertEqual(proj_dict['ellps'], 'WGS84') + new_ds, grid_mapping = area2gridmapping(ds) self.assertIn('crs_wkt', grid_mapping.attrs) - self.assertEqual(ds.attrs['grid_mapping'], 'cosmo7') + wkt = grid_mapping.attrs['crs_wkt'] + self.assertIn('ELLIPSOID["WGS 84"', wkt) + self.assertIn('PARAMETER["lat_0",46', wkt) + self.assertIn('PARAMETER["lon_0",4.535', wkt) + self.assertIn('PARAMETER["o_lat_p",90', wkt) + self.assertIn('PARAMETER["o_lon_p",-5.465', wkt) + self.assertEqual(new_ds.attrs['grid_mapping'], 'cosmo7') # c) Projection Transverse Mercator lat_0 = 36.5 @@ -815,8 +816,8 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = tmerc - grid_mapping = area2gridmapping(ds) - self.assertEqual(ds.attrs['grid_mapping'], 'tmerc') + new_ds, grid_mapping = area2gridmapping(ds) + self.assertEqual(new_ds.attrs['grid_mapping'], 'tmerc') _gm_matches(grid_mapping, tmerc_expected) # d) Projection that has a representation but no explicit a/b @@ -841,9 +842,9 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = geos - grid_mapping = area2gridmapping(ds) + new_ds, grid_mapping = area2gridmapping(ds) - self.assertEqual(ds.attrs['grid_mapping'], 'geos') + self.assertEqual(new_ds.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) # e) oblique Mercator @@ -872,9 +873,9 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = area - grid_mapping = area2gridmapping(ds) + new_ds, grid_mapping = area2gridmapping(ds) - self.assertEqual(ds.attrs['grid_mapping'], 'omerc_otf') + self.assertEqual(new_ds.attrs['grid_mapping'], 'omerc_otf') _gm_matches(grid_mapping, omerc_expected) # f) Projection that has a representation but no explicit a/b @@ -897,9 +898,9 @@ def _gm_matches(gmapping, expected): ds = ds_base.copy() ds.attrs['area'] = geos - grid_mapping = area2gridmapping(ds) + new_ds, grid_mapping = area2gridmapping(ds) - self.assertEqual(ds.attrs['grid_mapping'], 'geos') + self.assertEqual(new_ds.attrs['grid_mapping'], 'geos') _gm_matches(grid_mapping, geos_expected) def test_area2lonlat(self): @@ -919,13 +920,13 @@ def test_area2lonlat(self): lons_ref, lats_ref = area.get_lonlats() dataarray = xr.DataArray(data=[[1, 2], [3, 4]], dims=('y', 'x'), attrs={'area': area}) - area2lonlat(dataarray) - res = [dataarray] # modified in-place + res = area2lonlat(dataarray) - self.assertEqual(len(res), 1) - self.assertEqual(set(res[0].coords), {'longitude', 'latitude'}) - lat = res[0]['latitude'] - lon = res[0]['longitude'] + # original should be unmodified + self.assertNotIn('longitude', dataarray.coords) + self.assertEqual(set(res.coords), {'longitude', 'latitude'}) + lat = res['latitude'] + lon = res['longitude'] self.assertTrue(np.all(lat.data == lats_ref)) self.assertTrue(np.all(lon.data == lons_ref)) self.assertDictContainsSubset({'name': 'latitude', 'standard_name': 'latitude', 'units': 'degrees_north'}, diff --git a/satpy/writers/cf_writer.py b/satpy/writers/cf_writer.py index 0c659911df..aad9f04fba 100644 --- a/satpy/writers/cf_writer.py +++ b/satpy/writers/cf_writer.py @@ -165,6 +165,7 @@ def get_extra_ds(dataset): def area2lonlat(dataarray): """Convert an area to longitudes and latitudes.""" + dataarray = dataarray.copy() area = dataarray.attrs['area'] chunks = getattr(dataarray.data, 'chunks', None) lons, lats = area.get_lonlats(chunks=chunks) @@ -180,14 +181,16 @@ def area2lonlat(dataarray): name='latitude') dataarray['longitude'] = lons dataarray['latitude'] = lats + return dataarray def area2gridmapping(dataarray): """Convert an area to at CF grid mapping.""" + dataarray = dataarray.copy() area = dataarray.attrs['area'] gmapping_var_name, attrs = create_grid_mapping(area) dataarray.attrs['grid_mapping'] = gmapping_var_name - return xr.DataArray(0, attrs=attrs, name=gmapping_var_name) + return dataarray, xr.DataArray(0, attrs=attrs, name=gmapping_var_name) def area2cf(dataarray, strict=False): @@ -195,11 +198,10 @@ def area2cf(dataarray, strict=False): res = [] dataarray = dataarray.copy(deep=True) if isinstance(dataarray.attrs['area'], SwathDefinition) or strict: - # modifies dataarray in-place - area2lonlat(dataarray) + dataarray = area2lonlat(dataarray) if isinstance(dataarray.attrs['area'], AreaDefinition): - res.append(area2gridmapping(dataarray)) - + dataarray, gmapping = area2gridmapping(dataarray) + res.append(gmapping) res.append(dataarray) return res From 11ae1123513b5aa2e2c806fe5f7fb78d8a2511ab Mon Sep 17 00:00:00 2001 From: David Hoese Date: Tue, 31 Mar 2020 09:51:23 -0500 Subject: [PATCH 11/11] Move CF DataArray deep copy to be more obvious --- satpy/writers/cf_writer.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/satpy/writers/cf_writer.py b/satpy/writers/cf_writer.py index aad9f04fba..10ae9d8c69 100644 --- a/satpy/writers/cf_writer.py +++ b/satpy/writers/cf_writer.py @@ -196,7 +196,6 @@ def area2gridmapping(dataarray): def area2cf(dataarray, strict=False): """Convert an area to at CF grid mapping or lon and lats.""" res = [] - dataarray = dataarray.copy(deep=True) if isinstance(dataarray.attrs['area'], SwathDefinition) or strict: dataarray = area2lonlat(dataarray) if isinstance(dataarray.attrs['area'], AreaDefinition): @@ -391,6 +390,7 @@ def encode_attrs_nc(attrs): Attributes to be encoded Returns: dict: Encoded (and sorted) attributes + """ encoded_attrs = [] for key, val in sorted(attrs.items()): @@ -415,6 +415,7 @@ def da2cf(dataarray, epoch=EPOCH, flatten_attrs=False, exclude_attrs=None, compr If True, flatten dict-type attributes exclude_attrs (list): List of dataset attributes to be excluded + """ if exclude_attrs is None: exclude_attrs = [] @@ -489,13 +490,17 @@ def _collect_datasets(self, datasets, epoch=EPOCH, flatten_attrs=False, exclude_ datas = {} start_times = [] end_times = [] - for ds_name, ds in sorted(ds_collection.items()): + # sort by name, but don't use the name + for _, ds in sorted(ds_collection.items()): if ds.dtype not in CF_DTYPES: warnings.warn('Dtype {} not compatible with {}.'.format(str(ds.dtype), CF_VERSION)) + # we may be adding attributes, coordinates, or modifying the + # structure of attributes + ds = ds.copy(deep=True) try: new_datasets = area2cf(ds, strict=include_lonlats) except KeyError: - new_datasets = [ds.copy(deep=True)] + new_datasets = [ds] for new_ds in new_datasets: start_times.append(new_ds.attrs.get("start_time", None)) end_times.append(new_ds.attrs.get("end_time", None)) @@ -518,7 +523,7 @@ def update_encoding(self, dataset, to_netcdf_kwargs): other_to_netcdf_kwargs = to_netcdf_kwargs.copy() encoding = other_to_netcdf_kwargs.pop('encoding', {}).copy() coord_vars = [] - for name, data_array in dataset.items(): + for data_array in dataset.values(): coord_vars.extend(set(data_array.dims).intersection(data_array.coords)) for coord_var in coord_vars: encoding.setdefault(coord_var, {}) @@ -579,6 +584,7 @@ def save_datasets(self, datasets, filename=None, groups=None, header_attrs=None, Compression to use on the datasets before saving, for example {'zlib': True, 'complevel': 9}. This is in turn passed the xarray's `to_netcdf` method: http://xarray.pydata.org/en/stable/generated/xarray.Dataset.to_netcdf.html for more possibilities. + """ logger.info('Saving datasets to NetCDF4/CF.')