remove pyresample

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

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')
+    grid_mapping = pyproj.CRS.from_epsg(3411), 'absent'

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)):

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),

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__":