From f61a028de500d5cefbe6e2238875ae2bbe68c604 Mon Sep 17 00:00:00 2001 From: akorosov Date: Mon, 17 Jan 2022 14:46:16 +0100 Subject: [PATCH] remove pyresample --- environment.yml | 2 - geodataset/custom_geodataset.py | 49 ++------------- geodataset/geodataset.py | 96 ++++++++++++++++++++--------- geodataset/tests/test_geodataset.py | 90 ++++++++++++++++++++++----- geodataset/tests/test_tools.py | 15 ++--- geodataset/tools.py | 12 ++-- 6 files changed, 162 insertions(+), 102 deletions(-) diff --git a/environment.yml b/environment.yml index e34245e..0f79c37 100644 --- a/environment.yml +++ b/environment.yml @@ -5,7 +5,6 @@ dependencies: - cartopy - cython - pip - - pykdtree - pip: - cmocean - gdown @@ -15,7 +14,6 @@ dependencies: - netcdftime - numpy - pyproj - - pyresample - pytest - pyyaml - scipy diff --git a/geodataset/custom_geodataset.py b/geodataset/custom_geodataset.py index 4f9defe..efed448 100644 --- a/geodataset/custom_geodataset.py +++ b/geodataset/custom_geodataset.py @@ -4,7 +4,7 @@ import numpy as np import pyproj -from geodataset.geodataset import GeoDatasetRead, GeoDatasetWrite +from geodataset.geodataset import GeoDatasetRead from geodataset.utils import InvalidDatasetError class CustomDatasetRead(GeoDatasetRead): @@ -18,11 +18,6 @@ def _check_input_file(self): class CmemsMetIceChart(CustomDatasetRead): pattern = re.compile(r'ice_conc_svalbard_\d{12}.nc') lonlat_names = 'lon', 'lat' - grid_mapping_variable = 'crs' - - @property - def projection(self): - return pyproj.Proj(self.variables['crs'].proj4_string) class Dist2Coast(CustomDatasetRead): @@ -42,24 +37,7 @@ def get_lonlat_arrays(self): class JaxaAmsr2IceConc(CustomDatasetRead): pattern = re.compile(r'Arc_\d{8}_res3.125_pyres.nc') lonlat_names = 'longitude', 'latitude' - projection = pyproj.Proj(3411) - grid_mapping_variable = 'absent' - - -class MooringsNextsim(CustomDatasetRead): - pattern = re.compile(r'Moorings.*.nc') - projection = pyproj.Proj( - '+proj=stere +a=6378273.0 +b=6356889.448910593 ' - '+lon_0=-45.0 +lat_0=90.0 +lat_ts=60.0') - grid_mapping_variable = 'Polar_Stereographic_Grid' - - -class MooringsArcMfc(CustomDatasetRead): - pattern = re.compile(r'Moorings.*.nc') - projection = pyproj.Proj( - '+proj=stere +a=6378273.0 +b=6378273.0 ' - '+lon_0=-45.0 +lat_0=90.0 +lat_ts=90.0') - grid_mapping_variable = 'absent' + grid_mapping = pyproj.CRS.from_epsg(3411), 'absent' class NerscSarProducts(CustomDatasetRead): @@ -79,27 +57,12 @@ class NerscIceType(NerscSarProducts): class OsisafDriftersNextsim(CustomDatasetRead): pattern = re.compile(r'OSISAF_Drifters_.*.nc') - projection = pyproj.Proj("+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 " - " +a=6378273 +b=6356889.44891 ") - grid_mapping_variable = 'absent' + grid_mapping = pyproj.CRS.from_proj4( + " +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 " + " +a=6378273 +b=6356889.44891 "), 'absent' is_lonlat_2d = False class SmosIceThickness(CustomDatasetRead): pattern = re.compile(r'SMOS_Icethickness_v3.2_north_\d{8}.nc') - projection = pyproj.Proj(3411) - grid_mapping_variable = 'absent' - - -class Topaz4Forecast(CustomDatasetRead): - pattern = re.compile(r'\d{8}_dm-metno-MODEL-topaz4-ARC-b\d{8}-fv02.0.nc') - projection = pyproj.Proj("+proj=stere +lat_0=90 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") - grid_mapping_variable = 'stereographic' - - -class NetcdfArcMFC(GeoDatasetWrite): - """ wrapper for netCDF4.Dataset with info about ArcMFC products """ - grid_mapping_variable = 'stereographic' - projection = pyproj.Proj( - '+proj=stere +a=6378273 +b=6378273.0 ' - ' +lon_0=-45 +lat_0=90 +lat_ts=90') \ No newline at end of file + grid_mapping = pyproj.CRS.from_epsg(3411), 'absent' diff --git a/geodataset/geodataset.py b/geodataset/geodataset.py index 60b2a3f..c295dcf 100644 --- a/geodataset/geodataset.py +++ b/geodataset/geodataset.py @@ -6,7 +6,6 @@ import numpy as np import pyproj from pyproj.exceptions import CRSError -from pyresample.utils import load_cf_area from scipy.interpolate import RegularGridInterpolator from xarray.core.variable import MissingDimensionsError @@ -292,54 +291,93 @@ def variable_names(self): @property def projection(self): - """ Read projection of the dataset from self.area_definition + """ Read projection of the dataset from self.grid_mapping Returns ------- projection : pyproj.Proj """ - return pyproj.Proj(self.area_definition.crs) - - @property - def area_definition(self): - """ Read area definition of the dataset from self._area_def_cf_info - - Returns - ------- - area_definition : pyresample.AreaDefinition - - """ - return self._area_def_cf_info[0] + return pyproj.Proj(self.grid_mapping[0]) @property def grid_mapping_variable(self): - """ Read name of the grid mapping variable from self._area_def_cf_info + """ Read name of the grid mapping variable from self.grid_mapping Returns ------- grid_mapping_variable : str """ - return self._area_def_cf_info[1]['grid_mapping_variable'] - + return self.grid_mapping[1] + @property @lru_cache(1) - def _area_def_cf_info(self): - """ Read area definition and grid mapping information using pyresample.load_cf_area + def grid_mapping(self): + """ Load CRS and grid mapping variable name from CF-attrinbutes OR from lon/lat + If grid mapping cannot be loaded from file, InvalidDatasetError is raised Returns ------- - area_def_cf_info : tuple - area_definition and grid_mapping_information - - """ - try: - area_def_cf_info = load_cf_area(self.filename) - except (MissingDimensionsError, CRSError, KeyError, ValueError) as e: - raise InvalidDatasetError - else: - return area_def_cf_info + csr : pyproj.CRS + coordinate reference system + v : str + name of grid_mapping_variable or "absent" + + """ + crs, v = self.get_grid_mapping_from_cf_attrs() + if not crs: + crs, v = self.get_grid_mapping_from_lonlat() + if not crs: + raise InvalidDatasetError + return crs, v + + def get_grid_mapping_from_cf_attrs(self): + """ Load CRS and grid mapping var name from CF-attributes + + Returns + ------- + csr : pyproj.CRS or None + coordinate reference system + v : str or None + name of grid mapping variable + + """ + for var_name, variable in self.variables.items(): + attrs = {attr:variable.getncattr(attr) for attr in variable.ncattrs()} + try: + crs = pyproj.CRS.from_cf(attrs) + except CRSError: + pass + else: + return crs, var_name + return None, None + + def get_grid_mapping_from_lonlat(self): + """ Check if longitude and latitude are dimentions and return longlat CRS and "absent", + otherwise return None, None + + Returns + ------- + csr : pyproj.CRS or None + coordinate reference system + v : str + name of grid mapping variable + + """ + lon_is_dim = False + lat_is_dim = False + for d in self.dimensions: + if d in self.variables: + if 'standard_name' in self.variables[d].ncattrs(): + if self.variables[d].standard_name == 'longitude': + lon_is_dim = True + elif self.variables[d].standard_name == 'latitude': + lat_is_dim = True + if lon_is_dim and lat_is_dim: + return pyproj.CRS( + '+proj=longlat +datum=WGS84 +no_defs +type=crs'), 'absent' + return None, None def get_variable_array( self, var_name, time_index=0, ij_range=(None, None, None, None)): diff --git a/geodataset/tests/test_geodataset.py b/geodataset/tests/test_geodataset.py index 070acb5..5100f18 100644 --- a/geodataset/tests/test_geodataset.py +++ b/geodataset/tests/test_geodataset.py @@ -8,7 +8,7 @@ from netCDF4 import Dataset import numpy as np import pyproj -from pyresample.utils import load_cf_area +from pyproj.exceptions import CRSError from geodataset.geodataset import GeoDatasetBase, GeoDatasetWrite, GeoDatasetRead from geodataset.utils import InvalidDatasetError @@ -279,11 +279,9 @@ def test_get_lonlat_names(self, **kwargs): variables['lon'].standard_name = 'longitude' variables['lat'].standard_name = 'latitude' variables['sic'].standard_name = 'sea_ice_concentration' - with GeoDatasetRead() as ds: ds.variables = variables lon_name, lat_name = ds.lonlat_names - self.assertEqual(lon_name, 'lon') self.assertEqual(lat_name, 'lat') @@ -298,27 +296,87 @@ def test_get_lonlat_names_raises(self, **kwargs): } variables['lon'].ncattrs.return_value = ['standard_name', 'a', 'b'] variables['lon'].standard_name = 'longitude' - with GeoDatasetRead() as ds: ds.variables = variables with self.assertRaises(InvalidDatasetError): lon_name, lat_name = ds.lonlat_names - @patch('geodataset.geodataset.load_cf_area') - def test_area_def_cf_info(self, mock_lca): - with GeoDatasetRead(self.osisaf_filename) as ds: - a = ds._area_def_cf_info - b = ds._area_def_cf_info - mock_lca.assert_called_once_with(ds.filename) - def test_grid_mapping_variable(self): with GeoDatasetRead(self.osisaf_filename) as ds: self.assertEqual(ds.grid_mapping_variable, 'Polar_Stereographic_Grid') - - def test_projection(self): - p = pyproj.Proj(load_cf_area(self.osisaf_filename)[0].crs) - with GeoDatasetRead(self.osisaf_filename) as ds: - self.assertEqual(ds.projection, p) + + @patch.multiple(GeoDatasetRead, + __init__=MagicMock(return_value=None), + grid_mapping=["gm_crs", "gm_var"], + ) + @patch('geodataset.geodataset.pyproj.Proj') + def test_projection(self, mock_Proj): + ds = GeoDatasetRead() + p = ds.projection + mock_Proj.assert_called_once_with('gm_crs') + + @patch.multiple(GeoDatasetRead, + __init__=MagicMock(return_value=None), + grid_mapping=["gm_crs", "gm_var"], + ) + def test_grid_mapping_variable(self): + ds = GeoDatasetRead() + self.assertEqual(ds.grid_mapping_variable, 'gm_var') + + @patch.multiple(GeoDatasetRead, + __init__=MagicMock(return_value=None), + get_grid_mapping_from_cf_attrs=DEFAULT, + get_grid_mapping_from_lonlat=DEFAULT, + ) + def test_grid_mapping(self, **kwargs): + GeoDatasetRead.get_grid_mapping_from_cf_attrs.return_value = None, None + GeoDatasetRead.get_grid_mapping_from_lonlat.return_value = None, None + ds = GeoDatasetRead() + with self.assertRaises(InvalidDatasetError): + gm = ds.grid_mapping + GeoDatasetRead.get_grid_mapping_from_lonlat.return_value = "longlat_crs", "absent" + ds = GeoDatasetRead() + self.assertEqual(ds.grid_mapping, ("longlat_crs", "absent")) + GeoDatasetRead.get_grid_mapping_from_cf_attrs.return_value = 'crs', 'gm_var_name' + ds = GeoDatasetRead() + self.assertEqual(ds.grid_mapping, ('crs', 'gm_var_name')) + + @patch.multiple(GeoDatasetRead, + __init__=MagicMock(return_value=None), + variables=DEFAULT) + @patch('geodataset.geodataset.pyproj.CRS') + def test_get_grid_mapping_from_cf_attrs(self, mock_CRS, **kwargs): + GeoDatasetRead.variables = { + 'var_name': MagicMock(**{ + 'ncattrs.return_value':['attr_name'], + 'getncattr.return_value': 'attr_val'})} + ds = GeoDatasetRead() + mock_CRS.from_cf.side_effect = CRSError('error message') + crs, varname = ds.get_grid_mapping_from_cf_attrs() + self.assertEqual((crs, varname), (None, None)) + mock_CRS.from_cf.return_value = 'crs' + mock_CRS.from_cf.side_effect = None + crs, varname = ds.get_grid_mapping_from_cf_attrs() + self.assertEqual((crs, varname), ('crs', 'var_name')) + + @patch.multiple(GeoDatasetRead, __init__=MagicMock(return_value=None), dimensions=DEFAULT, variables=DEFAULT) + def test_get_grid_mapping_from_lonlat(self, **kwargs): + GeoDatasetRead.dimensions = ['lon', 'lat'] + variables = {'lon': Mock(), 'lat': Mock()} + variables['lon'].ncattrs.return_value = ['standard_name'] + variables['lon'].standard_name = 'longitude' + variables['lat'].ncattrs.return_value = ['standard_name'] + variables['lat'].standard_name = 'latitude' + ds = GeoDatasetRead() + ds.variables = variables + crs, varname = ds.get_grid_mapping_from_lonlat() + self.assertEqual( + (crs, varname), + (pyproj.CRS('+proj=longlat +datum=WGS84 +no_defs +type=crs'), 'absent')) + ds = GeoDatasetRead() + ds.variables = {'bla': Mock(), 'blo': Mock()} + crs, varname = ds.get_grid_mapping_from_lonlat() + self.assertEqual((crs, varname), (None, None)) @patch.multiple(GeoDatasetRead, __init__=MagicMock(return_value=None), diff --git a/geodataset/tests/test_tools.py b/geodataset/tests/test_tools.py index 3cf1be4..82aae39 100644 --- a/geodataset/tests/test_tools.py +++ b/geodataset/tests/test_tools.py @@ -5,7 +5,6 @@ import pyproj from geodataset.tools import open_netcdf - from geodataset.tests.base_for_tests import BaseForTests @@ -16,12 +15,13 @@ def setUp(self): def test_open_netcdf(self): for nc_file in self.nc_files: - ds = open_netcdf(nc_file) - print(nc_file, ds.lonlat_names) - self.assertIsInstance(ds.lonlat_names[0], str) - self.assertIsInstance(ds.lonlat_names[1], str) - self.assertIsInstance(ds.variable_names, list) - self.assertIsInstance(ds.variable_names[0], str) + with self.subTest(nc_file=nc_file): + ds = open_netcdf(nc_file) + print(nc_file, ds.lonlat_names) + self.assertIsInstance(ds.lonlat_names[0], str) + self.assertIsInstance(ds.lonlat_names[1], str) + self.assertIsInstance(ds.variable_names, list) + self.assertIsInstance(ds.variable_names[0], str) def test_get_lonlat_arrays(self): for nc_file in self.nc_files: @@ -42,6 +42,7 @@ def test_projection(self): for nc_file in self.nc_files: with self.subTest(nc_file=nc_file): ds = open_netcdf(nc_file) + self.assertIsInstance(ds.grid_mapping_variable, str) self.assertIsInstance(ds.projection, pyproj.Proj) if __name__ == "__main__": diff --git a/geodataset/tools.py b/geodataset/tools.py index 2877b0d..3142fa4 100644 --- a/geodataset/tools.py +++ b/geodataset/tools.py @@ -5,12 +5,10 @@ Dist2Coast, Etopo, JaxaAmsr2IceConc, - MooringsNextsim, NerscDeformation, NerscIceType, OsisafDriftersNextsim, SmosIceThickness, - Topaz4Forecast, ) custom_read_classes = [ @@ -18,18 +16,22 @@ Dist2Coast, Etopo, JaxaAmsr2IceConc, - MooringsNextsim, NerscDeformation, NerscIceType, OsisafDriftersNextsim, SmosIceThickness, - Topaz4Forecast, # always last: GeoDatasetRead, ] def open_netcdf(file_address): - + """ Open NetCDF with read access and add geospatial metadata + + Returns + ------- + ds : GeoDataset or custom children + similar to netCDF4.Dataset with geospatial metadata and methods + """ for class_ in custom_read_classes: try: obj = class_(file_address)