Skip to content

Commit

Permalink
Merge pull request #230 from vincelhx/branch_dev_vinc
Browse files Browse the repository at this point in the history 
Better gestion of error for recalibration
  • Loading branch information
@vincelhx
vincelhx authored Nov 25, 2024
2 parents 3412f47 + 3b4bd4c commit 2e87b95
Show file tree
Hide file tree
Showing 4 changed files with 107 additions and 83 deletions.
14 changes: 1 addition & 13 deletions src/xsar/config.yml
View file
Original file line number Diff line number Diff line change
@@ -1,16 +1,4 @@
# default data dir for tests and examples
data_dir: /tmp
auxiliary_dir:
path_auxiliary_df:

auxiliary_names:
S1B:
v_IPF_36:
AUX_PP1: 'S1B_AUX_PP1_V20160422T000000_G20220323T140710.SAFE'
AUX_CAL: 'S1B_AUX_CAL_V20160422T000000_G20210104T140113.SAFE'
#AUX_INS: 'S1B_AUX_INS_V20160422T000000_G20211027T134314.SAFE'
S1A:
v_IPF_36:
AUX_PP1: 'S1A_AUX_PP1_V20190228T092500_G20220323T153041.SAFE'
AUX_CAL: 'S1A_AUX_CAL_V20190228T092500_G20210104T141310.SAFE'
#AUX_INS: 'S1A_AUX_INS_V20190228T092500_G20211103T111906.SAFE'}
path_auxiliary_df:
21 changes: 10 additions & 11 deletions src/xsar/raster_readers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,12 +73,10 @@ def ecmwf_0100_1h(fname, **kwargs):
xarray.Dataset
"""
ecmwf_ds = xr.open_dataset(
fname, chunks={"Longitude": 1000, "Latitude": 1000}
).isel(time=0)
ecmwf_ds.attrs["time"] = datetime.datetime.fromtimestamp(
ecmwf_ds.time.item() // 1000000000
)
if "time" in ecmwf_ds:
fname, chunks={'time': 1, 'latitude': 901, 'longitude': 1800}).isel(time=0)
ecmwf_ds.attrs['time'] = datetime.datetime.fromtimestamp(
ecmwf_ds.time.item() // 1000000000)
if 'time' in ecmwf_ds:
ecmwf_ds = ecmwf_ds.drop("time")
ecmwf_ds = ecmwf_ds[["Longitude", "Latitude", "10U", "10V"]].rename(
{"Longitude": "x", "Latitude": "y", "10U": "U10", "10V": "V10"}
Expand Down Expand Up @@ -143,7 +141,7 @@ def hwrf_0015_3h(fname, **kwargs):
-------
xarray.Dataset
"""
hwrf_ds = xr.open_dataset(fname)
hwrf_ds = xr.open_dataset(fname, chunks={'t': 1, 'LAT': 700, 'LON': 700})
try:
hwrf_ds = hwrf_ds[["U", "V", "LON", "LAT"]]
hwrf_ds = hwrf_ds.squeeze("t", drop=True)
Expand Down Expand Up @@ -186,10 +184,11 @@ def era5_0250_1h(fname, **kwargs):
xarray.Dataset
"""

ds_era5 = xr.open_dataset(fname, chunks={"time": 1})
ds_era5 = ds_era5[["u10", "v10", "latitude025", "longitude025"]]
ds_era5 = ds_era5.sel(time=str(kwargs["date"]), method="nearest")
ds_era5 = ds_era5.drop("time")
ds_era5 = xr.open_dataset(
fname, chunks={'time': 6, 'latitude025': 721, 'longitude025': 1440})
ds_era5 = ds_era5[['u10', 'v10', 'latitude025', 'longitude025']]
ds_era5 = ds_era5.sel(time=str(kwargs['date']), method="nearest")
ds_era5 = ds_era5.drop('time')

ds_era5 = ds_era5.rename(
{"longitude025": "x", "latitude025": "y", "u10": "U10", "v10": "V10"}
Expand Down
91 changes: 60 additions & 31 deletions src/xsar/rcm_dataset.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,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
Expand Down Expand Up @@ -99,7 +100,8 @@ def __init__(
# assert isinstance(rs2meta.coords2ll(100, 100),tuple)
else:
# we want self.rs2meta to be a dask actor on a worker
self.sar_meta = BlockingActorProxy(RcmMeta.from_dict, dataset_id.dict)
self.sar_meta = BlockingActorProxy(
RcmMeta.from_dict, dataset_id.dict)
del dataset_id

if self.sar_meta.multidataset:
Expand Down Expand Up @@ -146,7 +148,8 @@ def __init__(
"beta0": "beta0",
}

geoloc_vars = ["latitude", "longitude", "altitude", "incidence", "elevation"]
geoloc_vars = ["latitude", "longitude",
"altitude", "incidence", "elevation"]
for vv in skip_variables:
if vv in geoloc_vars:
geoloc_vars.remove(vv)
Expand All @@ -168,8 +171,10 @@ def __init__(
value_res_sample = float(resolution.replace("m", ""))
value_res_line = value_res_sample
elif isinstance(resolution, dict):
value_res_sample = self.sar_meta.pixel_sample_m * resolution["sample"]
value_res_line = self.sar_meta.pixel_line_m * resolution["line"]
value_res_sample = self.sar_meta.pixel_sample_m * \
resolution["sample"]
value_res_line = self.sar_meta.pixel_line_m * \
resolution["line"]
else:
logger.warning(
"resolution type not handle (%s) should be str or dict -> sampleSpacing"
Expand Down Expand Up @@ -198,7 +203,8 @@ def __init__(
dask.array.empty_like(
self._dataset.digital_number.isel(pol=0).drop("pol"),
dtype=np.int8,
name="empty_var_tmpl-%s" % dask.base.tokenize(self.sar_meta.name),
name="empty_var_tmpl-%s" % dask.base.tokenize(
self.sar_meta.name),
),
dims=("line", "sample"),
coords={
Expand Down Expand Up @@ -234,7 +240,8 @@ def __init__(

self._luts = self.lazy_load_luts()
self._noise_luts = self.lazy_load_noise_luts()
self._noise_luts = self._noise_luts.drop(["pixelFirstNoiseValue", "stepSize"])
self._noise_luts = self._noise_luts.drop(
["pixelFirstNoiseValue", "stepSize"])
self.apply_calibration_and_denoising()
self._dataset = xr.merge(
[
Expand All @@ -249,9 +256,11 @@ def __init__(
if rasters is not None:
self._dataset = xr.merge([self._dataset, rasters])
if "ground_heading" not in skip_variables:
self._dataset = xr.merge([self.load_ground_heading(), self._dataset])
self._dataset = xr.merge(
[self.load_ground_heading(), self._dataset])
if "velocity" not in skip_variables:
self._dataset = xr.merge([self.get_sensor_velocity(), self._dataset])
self._dataset = xr.merge(
[self.get_sensor_velocity(), self._dataset])
self._rasterized_masks = self.load_rasterized_masks()
self._dataset = xr.merge([self._rasterized_masks, self._dataset])
self.datatree["measurement"] = self.datatree["measurement"].assign(
Expand Down Expand Up @@ -410,7 +419,6 @@ def _apply_calibration_lut(self, var_name):
lut = self._get_lut(var_name)
offset = lut.attrs["offset"]
res = ((self._dataset.digital_number**2.0) + offset) / lut
res = res.where(res > 0)
res.attrs.update(lut.attrs)
return res.to_dataset(name=var_name + "_raw")

Expand Down Expand Up @@ -464,10 +472,22 @@ def apply_calibration_and_denoising(self):
"Skipping variable '%s' ('%s' lut is missing)"
% (var_name, lut_name)
)

self._dataset = self._add_denoised(self._dataset)

for var_name, lut_name in self._map_var_lut.items():
var_name_raw = var_name + "_raw"
if var_name_raw in self._dataset:
self._dataset[var_name_raw] = self._dataset[var_name_raw].where(
self._dataset[var_name_raw] > 0, 0)
else:
logger.debug(
"Skipping variable '%s' ('%s' lut is missing)"
% (var_name, lut_name)
)
self.datatree["measurement"] = self.datatree["measurement"].assign(
self._dataset
)
self._dataset)

# self._dataset = self.datatree[
# 'measurement'].to_dataset() # test oct 22 to see if then I can modify variables of the dt
return
Expand Down Expand Up @@ -496,7 +516,6 @@ def _add_denoised(self, ds, clip=False, vars=None):
if vars is None:
vars = ["sigma0", "beta0", "gamma0"]
for varname in vars:

varname_raw = varname + "_raw"
noise = "ne%sz" % varname[0]
if varname_raw not in ds:
Expand All @@ -511,14 +530,18 @@ def _add_denoised(self, ds, clip=False, vars=None):
else:
denoised.attrs["comment"] = "not clipped, some values can be <0"
ds[varname] = denoised

ds[varname_raw].attrs[
"denoising information"
] = f"product was not denoised"

else:
logging.debug(
"product was already denoised (noiseSubtractionPerformed = True), noise added back"
)
denoised = ds[varname_raw]
denoised.attrs["denoising information"] = (
"already denoised by Canadian Agency"
"already denoised by Canadian Space Agency"
)
if clip:
denoised = denoised.clip(min=0)
Expand All @@ -530,7 +553,7 @@ def _add_denoised(self, ds, clip=False, vars=None):
ds[varname_raw] = denoised + ds[noise]
ds[varname_raw].attrs[
"denoising information"
] = "product was already denoised, noise added back"
] = "product was already denoised by Canadian Space Agency, noise added back"

return ds

Expand Down Expand Up @@ -561,7 +584,7 @@ def _load_incidence_from_lut(self):
# ds_lut_f_delayed.attrs = incidence.attrs
return da

@timing
@ timing
def _load_elevation_from_lut(self):
"""
Load elevation from lut.
Expand All @@ -583,7 +606,7 @@ def _load_elevation_from_lut(self):
inside = angle_rad * earth_radius / (earth_radius + satellite_height)
return np.degrees(np.arcsin(inside))

@timing
@ timing
def _get_offboresight_from_elevation(self):
"""
Compute offboresight angle.
Expand All @@ -609,7 +632,7 @@ def _get_offboresight_from_elevation(self):
"comment"
] = "built from elevation angle and latitude"

@timing
@ timing
def load_from_geoloc(self, varnames, lazy_loading=True):
"""
Interpolate (with RectBiVariateSpline) variables from `self.sar_meta.geoloc` to `self._dataset`
Expand Down Expand Up @@ -663,7 +686,8 @@ def load_from_geoloc(self, varnames, lazy_loading=True):
)
typee = self.sar_meta.geoloc[varname_in_geoloc].dtype
if lazy_loading:
da_var = map_blocks_coords(self._da_tmpl.astype(typee), interp_func)
da_var = map_blocks_coords(
self._da_tmpl.astype(typee), interp_func)
else:
da_val = interp_func(
self._dataset.digital_number.line,
Expand Down Expand Up @@ -695,7 +719,7 @@ def load_from_geoloc(self, varnames, lazy_loading=True):
ds = xr.merge(da_list)
return ds

@property
@ property
def interpolate_times(self):
"""
Apply interpolation with RectBivariateSpline to the azimuth time extracted from `self.sar_meta.geoloc`
Expand All @@ -712,7 +736,8 @@ def interpolate_times(self):
interp_func = RectBivariateSpline(
x=lines, y=samples, z=time_values_2d.astype(float), kx=1, ky=1
)
da_var = map_blocks_coords(self._da_tmpl.astype("datetime64[ns]"), interp_func)
da_var = map_blocks_coords(
self._da_tmpl.astype("datetime64[ns]"), interp_func)
return da_var.isel(sample=0).to_dataset(name="time")

def get_sensor_velocity(self):
Expand All @@ -735,10 +760,11 @@ def get_sensor_velocity(self):
vels = np.sqrt(np.sum(velos, axis=0))
interp_f = interp1d(azimuth_times.astype(float), vels)
_vels = interp_f(self.interpolate_times["time"].astype(float))
res = xr.DataArray(_vels, dims=["line"], coords={"line": self.dataset.line})
res = xr.DataArray(_vels, dims=["line"], coords={
"line": self.dataset.line})
return xr.Dataset({"velocity": res})

@timing
@ timing
def load_digital_number(
self, resolution=None, chunks=None, resampling=rasterio.enums.Resampling.rms
):
Expand Down Expand Up @@ -898,8 +924,10 @@ def _sort_list_files_and_get_pols(list_tiff):
(resolution["sample"] - 1) / 2, (resolution["line"] - 1) / 2
)
scale = Affine.scale(
rio.width // resolution["sample"] * resolution["sample"] / out_shape[1],
rio.height // resolution["line"] * resolution["line"] / out_shape[0],
rio.width // resolution["sample"] *
resolution["sample"] / out_shape[1],
rio.height // resolution["line"] *
resolution["line"] / out_shape[0],
)
sample, _ = translate * scale * (dn.sample, 0)
_, line = translate * scale * (0, dn.line)
Expand All @@ -916,7 +944,8 @@ def _sort_list_files_and_get_pols(list_tiff):
"history": yaml.safe_dump(
{
var_name: get_glob(
[p.replace(self.sar_meta.path + "/", "") for p in tiff_files]
[p.replace(self.sar_meta.path + "/", "")
for p in tiff_files]
)
}
),
Expand All @@ -935,7 +964,7 @@ def __repr__(self):
intro = "full coverage"
return "<RcmDataset %s object>" % intro

@timing
@ timing
def flip_sample_da(self, ds):
"""
When a product is flipped, flip back data arrays (from a dataset) sample dimensions to respect the xsar
Expand Down Expand Up @@ -969,7 +998,7 @@ def flip_sample_da(self, ds):
new_ds = ds
return new_ds

@timing
@ timing
def flip_line_da(self, ds):
"""
Flip dataArrays (from a dataset) that depend on line dimension when a product is ascending, in order to
Expand Down Expand Up @@ -1004,7 +1033,7 @@ def reconfigure_reader_datatree(self):
self.datatree.attrs |= self.sar_meta.dt.attrs
return

@property
@ property
def dataset(self):
"""
`xarray.Dataset` representation of this `xsar.RcmDataset` object.
Expand All @@ -1015,7 +1044,7 @@ def dataset(self):
res.attrs = self.datatree.attrs
return res

@dataset.setter
@ dataset.setter
def dataset(self, ds):
if self.sar_meta.name == ds.attrs["name"]:
# check if new ds has changed coordinates
Expand All @@ -1032,6 +1061,6 @@ def dataset(self, ds):
else:
raise ValueError("dataset must be same kind as original one.")

@dataset.deleter
@ dataset.deleter
def dataset(self):
logger.debug("deleter dataset")
Loading

0 comments on commit 2e87b95

Please sign in to comment.