Merge pull request #224 from vincelhx/branch_dev_vinc

check on rasters + fix RCM bug

src/xsar/base_meta.py

@@ -22,7 +22,8 @@
 logger.addHandler(logging.NullHandler())
 
 # we know tiff as no geotransform : ignore warning
-warnings.filterwarnings("ignore", category=rasterio.errors.NotGeoreferencedWarning)
+warnings.filterwarnings(
+    "ignore", category=rasterio.errors.NotGeoreferencedWarning)
 
 # allow nan without warnings
 # some dask warnings are still non filtered: https://github.com/dask/dask/issues/3245
@@ -41,8 +42,6 @@ class BaseMeta(BaseDataset):
         'land': cartopy.feature.NaturalEarthFeature('physical', 'land', '10m')
     }
 
-    rasters = available_rasters.iloc[0:0].copy()
-
     _mask_features = {}
     _mask_intersecting_geometries = {}
     _mask_geometry = {}
@@ -61,6 +60,9 @@ class BaseMeta(BaseDataset):
     safe = None
     geoloc = None
 
+    def __init__(self):
+        self.rasters = available_rasters.iloc[0:0].copy()
+
     def _get_mask_feature(self, name):
         # internal method that returns a cartopy feature from a mask name
         if self._mask_features[name] is None:
@@ -80,14 +82,16 @@ def _get_mask_feature(self, name):
                     except KeyError:
                         crs_in = fshp.crs
                     crs_in = pyproj.CRS(crs_in)
-                proj_transform = pyproj.Transformer.from_crs(pyproj.CRS('EPSG:4326'), crs_in, always_xy=True).transform
+                proj_transform = pyproj.Transformer.from_crs(
+                    pyproj.CRS('EPSG:4326'), crs_in, always_xy=True).transform
                 footprint_crs = transform(proj_transform, self.footprint)
 
                 with warnings.catch_warnings():
                     # ignore "RuntimeWarning: Sequential read of iterator was interrupted. Resetting iterator."
                     warnings.simplefilter("ignore", RuntimeWarning)
                     feature = cartopy.feature.ShapelyFeature(
-                        gpd.read_file(feature, mask=footprint_crs).to_crs(epsg=4326).geometry,
+                        gpd.read_file(feature, mask=footprint_crs).to_crs(
+                            epsg=4326).geometry,
                         cartopy.crs.PlateCarree()
                     )
             if not isinstance(feature, cartopy.feature.Feature):
@@ -160,7 +164,8 @@ def get_mask(self, name, describe=False):
             descr = self._mask_features_raw[name]
             try:
                 # nice repr for a class (like 'cartopy.feature.NaturalEarthFeature land')
-                descr = '%s.%s %s' % (descr.__module__, descr.__class__.__name__, descr.name)
+                descr = '%s.%s %s' % (
+                    descr.__module__, descr.__class__.__name__, descr.name)
             except AttributeError:
                 pass
             return descr
@@ -343,11 +348,14 @@ def ll2coords(self, *args):
         lon, lat = args
 
         # approximation with global inaccurate transform
-        line_approx, sample_approx = ~self.approx_transform * (np.asarray(lon), np.asarray(lat))
+        line_approx, sample_approx = ~self.approx_transform * \
+            (np.asarray(lon), np.asarray(lat))
 
         # Theoretical identity. It should be the same, but the difference show the error.
-        lon_identity, lat_identity = self.coords2ll(line_approx, sample_approx, to_grid=False)
-        line_identity, sample_identity = ~self.approx_transform * (lon_identity, lat_identity)
+        lon_identity, lat_identity = self.coords2ll(
+            line_approx, sample_approx, to_grid=False)
+        line_identity, sample_identity = ~self.approx_transform * \
+            (lon_identity, lat_identity)
 
         # we are now able to compute the error, and make a correction
         line_error = line_identity - line_approx
@@ -381,8 +389,10 @@ def coords2heading(self, lines, samples, to_grid=False, approx=True):
 
         """
 
-        lon1, lat1 = self.coords2ll(lines - 1, samples, to_grid=to_grid, approx=approx)
-        lon2, lat2 = self.coords2ll(lines + 1, samples, to_grid=to_grid, approx=approx)
+        lon1, lat1 = self.coords2ll(
+            lines - 1, samples, to_grid=to_grid, approx=approx)
+        lon2, lat2 = self.coords2ll(
+            lines + 1, samples, to_grid=to_grid, approx=approx)
         _, heading = haversine(lon1, lat1, lon2, lat2)
         return heading
 
@@ -424,9 +434,12 @@ def set_raster(self_or_cls, name, resource, read_function=None, get_function=Non
         default = available_rasters.loc[name:name]
 
         # set from params, or from default
-        self_or_cls.rasters.loc[name, 'resource'] = resource or default.loc[name, 'resource']
-        self_or_cls.rasters.loc[name, 'read_function'] = read_function or default.loc[name, 'read_function']
-        self_or_cls.rasters.loc[name, 'get_function'] = get_function or default.loc[name, 'get_function']
+        self_or_cls.rasters.loc[name,
+                                'resource'] = resource or default.loc[name, 'resource']
+        self_or_cls.rasters.loc[name,
+                                'read_function'] = read_function or default.loc[name, 'read_function']
+        self_or_cls.rasters.loc[name,
+                                'get_function'] = get_function or default.loc[name, 'get_function']
 
         return
 
@@ -459,4 +472,4 @@ def from_dict(cls, minidict):
         minidict = copy.copy(minidict)
         new = cls(minidict['name'])
         new.__dict__.update(minidict)
-        return new
+        return new

src/xsar/radarsat2_meta.py

@@ -30,6 +30,8 @@ class RadarSat2Meta(BaseMeta):
 
     @timing
     def __init__(self, name):
+        super().__init__()
+
         from xradarsat2 import rs2_reader
         if ':' in name:
             self.dt = rs2_reader(name.split(':')[1])
@@ -95,7 +97,8 @@ def _create_manifest_attrs(self):
             footprint_dict[ll] = [
                 self.geoloc[ll].isel(line=a, pixel=x).values for a, x in [(0, 0), (0, -1), (-1, -1), (-1, 0)]
             ]
-        corners = list(zip(footprint_dict['longitude'], footprint_dict['latitude']))
+        corners = list(
+            zip(footprint_dict['longitude'], footprint_dict['latitude']))
         p = Polygon(corners)
         self.geoloc.attrs['footprint'] = p
         dic["footprints"] = p
@@ -158,7 +161,7 @@ def to_dict(self, keys='minimal'):
 
         info_keys = {
             'minimal': [
-                #'platform',
+                # 'platform',
                 'swath', 'product', 'pols']
         }
         info_keys['all'] = info_keys['minimal'] + ['name', 'start_date', 'stop_date',
@@ -167,8 +170,8 @@ def to_dict(self, keys='minimal'):
                                                    'pixel_line_m', 'pixel_sample_m',
                                                    'approx_transform',
 
-                                                   #'orbit_pass',
-                                                   #'platform_heading'
+                                                   # 'orbit_pass',
+                                                   # 'platform_heading'
                                                    ]
 
         if isinstance(keys, str):
@@ -181,7 +184,8 @@ def to_dict(self, keys='minimal'):
             elif k in self.manifest_attrs.keys():
                 res_dict[k] = self.manifest_attrs[k]
             else:
-                raise KeyError('Unable to find key/attr "%s" in RadarSat2Meta' % k)
+                raise KeyError(
+                    'Unable to find key/attr "%s" in RadarSat2Meta' % k)
         return res_dict
 
     @property
@@ -214,7 +218,8 @@ def pixel_sample_m(self):
 
     def _get_time_range(self):
         if self.multidataset:
-            time_range = [self.manifest_attrs['start_date'], self.manifest_attrs['stop_date']]
+            time_range = [self.manifest_attrs['start_date'],
+                          self.manifest_attrs['stop_date']]
         else:
             time_range = self.orbit_and_attitude.timeStamp
         return pd.Interval(left=pd.Timestamp(time_range.values[0]), right=pd.Timestamp(time_range.values[-1]), closed='both')
@@ -275,7 +280,8 @@ def _dict_coords2ll(self):
         idx_line = np.array(geoloc.line)
 
         for ll in ['longitude', 'latitude']:
-            resdict[ll] = RectBivariateSpline(idx_line, idx_sample, np.asarray(geoloc[ll]), kx=1, ky=1)
+            resdict[ll] = RectBivariateSpline(
+                idx_line, idx_sample, np.asarray(geoloc[ll]), kx=1, ky=1)
 
         return resdict
 
@@ -291,7 +297,8 @@ def flip_sample_da(self):
         if (antenna_pointing, pass_direction) in flipped_cases:
             for ds_name in samples_depending_ds:
                 if 'radar' in ds_name:
-                    self.dt[ds_name] = self.dt[ds_name].rename({'NbOfNoiseLevelValues': 'pixel'})
+                    self.dt[ds_name] = self.dt[ds_name].rename(
+                        {'NbOfNoiseLevelValues': 'pixel'})
                 self.dt[ds_name] = self.dt[ds_name].copy().isel(pixel=slice(None, None, -1))\
                     .assign_coords(pixel=self.dt[ds_name].ds.pixel)
             self.samples_flipped = True
@@ -313,6 +320,3 @@ def flip_line_da(self):
                     .assign_coords(line=self.dt[ds_name].ds.line)
             self.lines_flipped = True
         return
-
-
-

src/xsar/rcm_meta.py

@@ -27,6 +27,8 @@ class RcmMeta(BaseMeta):
 
     @timing
     def __init__(self, name):
+        super().__init__()
+
         from safe_rcm import api
         if ':' in name:
             self.dt = api.open_rcm(name.split(':')[1])
@@ -68,7 +70,8 @@ def __init__(self, name):
         # replace `coefficients` dim by `pixel` in the datatree
         var_with_coeff = ['noise_lut', 'lut', 'incidence']
         for var in var_with_coeff:
-            self.dt[self._xpath[var]] = self.dt[self._xpath[var]].rename({'coefficients': 'pixel'})
+            self.dt[self._xpath[var]] = self.dt[self._xpath[var]].rename(
+                {'coefficients': 'pixel'})
         # build pixels coords
         self.assign_index(var_with_coeff)
 
@@ -103,7 +106,7 @@ def flip_sample_da(self):
         antenna_pointing = self.dt['sourceAttributes/radarParameters'].attrs['antennaPointing']
         pass_direction = self.dt['sourceAttributes/orbitAndAttitude/orbitInformation'].attrs['passDirection']
         flipped_cases = [('Left', 'Ascending'), ('Right', 'Descending')]
-        samples_depending_ds = ['geolocationGrid', 'lut', 'noise_lut']
+        samples_depending_ds = ['geolocationGrid']
         if (antenna_pointing, pass_direction) in flipped_cases:
             for ds_name in samples_depending_ds:
                 self.dt[self._xpath[ds_name]] = self.dt[self._xpath[ds_name]].copy().isel(pixel=slice(None, None, -1))\
@@ -182,7 +185,8 @@ def _assign_index_incidence():
             return
 
         if not all(item in ['lut', 'noise_lut', 'incidence'] for item in ds_list):
-            raise ValueError("Please enter accepted dataset names ('lut', 'noise_lut', 'incidence')")
+            raise ValueError(
+                "Please enter accepted dataset names ('lut', 'noise_lut', 'incidence')")
         if 'lut' in ds_list:
             _assign_index_lut()
         if 'noise_lut' in ds_list:
@@ -205,7 +209,8 @@ def _create_manifest_attrs(self):
             footprint_dict[ll] = [
                 self.geoloc[ll].isel(line=a, pixel=x).values for a, x in [(0, 0), (0, -1), (-1, -1), (-1, 0)]
             ]
-        corners = list(zip(footprint_dict['longitude'], footprint_dict['latitude']))
+        corners = list(
+            zip(footprint_dict['longitude'], footprint_dict['latitude']))
         p = Polygon(corners)
         self.geoloc.attrs['footprint'] = p
         dic["footprints"] = p
@@ -319,7 +324,8 @@ def __repr__(self):
 
     def _get_time_range(self):
         if self.multidataset:
-            time_range = [self.manifest_attrs['start_date'], self.manifest_attrs['stop_date']]
+            time_range = [self.manifest_attrs['start_date'],
+                          self.manifest_attrs['stop_date']]
         else:
             time_range = self.orbit.timeStamp
         return pd.Interval(left=pd.Timestamp(time_range.values[0]), right=pd.Timestamp(time_range.values[-1]), closed='both')
@@ -349,15 +355,16 @@ def _dict_coords2ll(self):
         idx_line = np.array(geoloc.line)
 
         for ll in ['longitude', 'latitude']:
-            resdict[ll] = RectBivariateSpline(idx_line, idx_sample, np.asarray(geoloc[ll]), kx=1, ky=1)
+            resdict[ll] = RectBivariateSpline(
+                idx_line, idx_sample, np.asarray(geoloc[ll]), kx=1, ky=1)
 
         return resdict
 
     def to_dict(self, keys='minimal'):
 
         info_keys = {
             'minimal': [
-                #'platform',
+                # 'platform',
                 'swath', 'product', 'pols']
         }
         info_keys['all'] = info_keys['minimal'] + ['name', 'start_date', 'stop_date',
@@ -366,8 +373,8 @@ def to_dict(self, keys='minimal'):
                                                    'pixel_line_m', 'pixel_sample_m',
                                                    'approx_transform',
 
-                                                   #'orbit_pass',
-                                                   #'platform_heading'
+                                                   # 'orbit_pass',
+                                                   # 'platform_heading'
                                                    ]
 
         if isinstance(keys, str):
@@ -381,4 +388,4 @@ def to_dict(self, keys='minimal'):
                 res_dict[k] = self.manifest_attrs[k]
             else:
                 raise KeyError('Unable to find key/attr "%s" in RcmMeta' % k)
-        return res_dict
+        return res_dict