Skip to content

Commit bb625d5

Browse files
agrouazeagrouaze
authored
Merge pull request #67 from umr-lops/newWW3hindcast
Fix L1C coloc with WW3spectra for tiles on land (i.e. subswath without SAR xspectra)
2 parents 72aa907 + 89d658b commit bb625d5

File tree

4 files changed

+245
-305
lines changed

4 files changed

+245
-305
lines changed

‎slcl1butils/coloc/coloc_IW_WW3spectra.py

+97-95
Original file line numberDiff line numberDiff line change
@@ -35,105 +35,107 @@ def resampleWW3spectra_on_TOPS_SAR_cartesian_grid(dsar, xspeckind):
3535
for d, k in dsar.sizes.items()
3636
if d in ["burst", "tile_sample", "tile_line"]
3737
}
38+
if (xspeckind == "intra" and "xspectra_2tau_Re" in dsar) or (xspeckind == "inter" and "xspectra_Re" in dsar): # in a future version of L1B xspectra variable could be always present (even on land) but filled by NaN
39+
# symmetrize and combine Re+1j*Im for all the xspectra SAR
40+
if xspeckind == "intra":
41+
xsSAR = dsar["xspectra_2tau_Re"] + 1j * dsar["xspectra_2tau_Im"]
42+
elif xspeckind == "inter":
43+
xsSAR = dsar["xspectra_Re"] + 1j * dsar["xspectra_Im"]
44+
else:
45+
raise ValueError("%s" % xspeckind)
46+
xsSAR = xsSAR.assign_coords(
47+
{
48+
"k_rg": xsSAR["k_rg"].mean(
49+
dim=set(xsSAR["k_rg"].dims) - set(["freq_sample"]),
50+
keep_attrs=True,
51+
)
52+
}
53+
).swap_dims({"freq_sample": "k_rg", "freq_line": "k_az"})
54+
if xspeckind == "intra":
55+
xsSAR = symmetrize_xspectrum(xsSAR).squeeze(dim="2tau")
56+
xsSAR.attrs = dsar["xspectra_2tau_Re"].attrs
57+
elif xspeckind == "inter":
58+
xsSAR = symmetrize_xspectrum(xsSAR)
59+
xsSAR.attrs = dsar["xspectra_Re"].attrs
60+
else:
61+
raise ValueError("%s" % xspeckind)
62+
# replace the half spectrum by a single variable, to save complexe values only possibility is zarr
3863

39-
# symmetrize and combine Re+1j*Im for all the xspectra SAR
40-
if xspeckind == "intra":
41-
xsSAR = dsar["xspectra_2tau_Re"] + 1j * dsar["xspectra_2tau_Im"]
42-
elif xspeckind == "inter":
43-
xsSAR = dsar["xspectra_Re"] + 1j * dsar["xspectra_Im"]
44-
else:
45-
raise ValueError("%s" % xspeckind)
46-
xsSAR = xsSAR.assign_coords(
47-
{
48-
"k_rg": xsSAR["k_rg"].mean(
49-
dim=set(xsSAR["k_rg"].dims) - set(["freq_sample"]),
50-
keep_attrs=True,
51-
)
52-
}
53-
).swap_dims({"freq_sample": "k_rg", "freq_line": "k_az"})
54-
if xspeckind == "intra":
55-
xsSAR = symmetrize_xspectrum(xsSAR).squeeze(dim="2tau")
56-
xsSAR.attrs = dsar["xspectra_2tau_Re"].attrs
57-
elif xspeckind == "inter":
58-
xsSAR = symmetrize_xspectrum(xsSAR)
59-
xsSAR.attrs = dsar["xspectra_Re"].attrs
60-
else:
61-
raise ValueError("%s" % xspeckind)
62-
# replace the half spectrum by a single variable, to save complexe values only possibility is zarr
63-
64-
list_ww3_cart_sp = []
65-
list_ww3_efth_sp = []
66-
dk_az = np.diff(xsSAR["k_az"])
67-
dk_rg = np.diff(xsSAR["k_rg"])
68-
dk = (dk_rg.mean(), dk_az.mean())
69-
# xsSAR["k_rg"].attrs["dkx"] = dk[0] # temporarily add dkx attrs
70-
# xsSAR["k_az"].attrs["dky"] = dk[1] # temporarily add dky attrs
71-
kmax = (
72-
np.abs(xsSAR["k_rg"]).max().item(),
73-
np.abs(xsSAR["k_az"]).max().item(),
74-
) # FN reviewed
64+
list_ww3_cart_sp = []
65+
list_ww3_efth_sp = []
66+
dk_az = np.diff(xsSAR["k_az"])
67+
dk_rg = np.diff(xsSAR["k_rg"])
68+
dk = (dk_rg.mean(), dk_az.mean())
69+
# xsSAR["k_rg"].attrs["dkx"] = dk[0] # temporarily add dkx attrs
70+
# xsSAR["k_az"].attrs["dky"] = dk[1] # temporarily add dky attrs
71+
kmax = (
72+
np.abs(xsSAR["k_rg"]).max().item(),
73+
np.abs(xsSAR["k_az"]).max().item(),
74+
) # FN reviewed
7575

76-
if pathww3sp:
77-
if os.path.exists(pathww3sp):
78-
flag_ww3spectra_found = True
79-
dsww3raw = xr.open_dataset(pathww3sp)
80-
for i in xndindex(gridsar):
81-
lonsar = dsar["longitude"][i].values
82-
latsar = dsar["latitude"][i].values
83-
heading = dsar["ground_heading"][i].values
84-
logging.debug("heading %s", heading)
85-
logging.debug("timesar :%s", start_date_dt)
86-
rotate = 90 + heading # deg clockwise wrt North
87-
logging.debug("rotate:%s", rotate)
76+
if pathww3sp:
77+
if os.path.exists(pathww3sp):
78+
flag_ww3spectra_found = True
79+
dsww3raw = xr.open_dataset(pathww3sp)
80+
for i in xndindex(gridsar):
81+
lonsar = dsar["longitude"][i].values
82+
latsar = dsar["latitude"][i].values
83+
heading = dsar["ground_heading"][i].values
84+
logging.debug("heading %s", heading)
85+
logging.debug("timesar :%s", start_date_dt)
86+
rotate = 90 + heading # deg clockwise wrt North
87+
logging.debug("rotate:%s", rotate)
8888

89-
# add the interpolated cartesian EFTH(kx,ky) spectra from WW3
90-
ds_ww3_cartesian = from_ww3(
91-
pathww3sp,
92-
dk=dk,
93-
kmax=kmax,
94-
strict="dk",
95-
rotate=rotate,
96-
clockwise_to_trigo=True,
97-
lon=lonsar,
98-
lat=latsar,
99-
time=start_date_dt,
100-
) # TODO use sensingTime
101-
ds_ww3_cartesian.attrs["source"] = "ww3"
102-
# TODO: check kx ky names to be same as the one from intra burst ds
103-
indiceww3spectra = find_closest_ww3(
104-
ww3_path=pathww3sp, lon=lonsar, lat=latsar, time=start_date_dt
105-
)
106-
# add the raw EFTH(f,dir) spectra from WW3
107-
rawspww3 = (
108-
dsww3raw["efth"].isel(time=indiceww3spectra).rename("ww3EFTHraw")
109-
)
110-
rawspww3.attrs["description"] = "raw EFTH(f,dir) spectra"
111-
ds_ww3_cartesian = ds_ww3_cartesian.swap_dims(
112-
{"kx": "k_rg", "ky": "k_az"}
113-
).T
114-
rawspww3 = rawspww3.assign_coords(i)
115-
rawspww3 = rawspww3.expand_dims(["burst", "tile_line", "tile_sample"])
116-
ds_ww3_cartesian = ds_ww3_cartesian.assign_coords(i)
117-
ds_ww3_cartesian = ds_ww3_cartesian.expand_dims(
118-
["burst", "tile_line", "tile_sample"]
119-
)
120-
list_ww3_cart_sp.append(ds_ww3_cartesian)
121-
list_ww3_efth_sp.append(rawspww3)
122-
ds_ww3_cartesian_merged = xr.merge(list_ww3_cart_sp)
123-
ds_ww3_efth_merged = xr.merge(list_ww3_efth_sp)
124-
dsar = xr.merge([dsar, ds_ww3_cartesian_merged, ds_ww3_efth_merged])
125-
flag_ww3spectra_added = True
126-
if xspeckind == "intra":
127-
dsar = dsar.drop_vars(["xspectra_2tau_Re", "xspectra_2tau_Im"])
128-
elif xspeckind == "inter":
129-
dsar = dsar.drop_vars(["xspectra_Re", "xspectra_Im"])
130-
dsar = dsar.drop_dims(["freq_sample", "freq_line"])
131-
if xspeckind == "intra":
132-
dsar["xspectra_2tau"] = xsSAR
133-
elif xspeckind == "inter":
134-
dsar["xspectra"] = xsSAR
89+
# add the interpolated cartesian EFTH(kx,ky) spectra from WW3
90+
ds_ww3_cartesian = from_ww3(
91+
pathww3sp,
92+
dk=dk,
93+
kmax=kmax,
94+
strict="dk",
95+
rotate=rotate,
96+
clockwise_to_trigo=True,
97+
lon=lonsar,
98+
lat=latsar,
99+
time=start_date_dt,
100+
) # TODO use sensingTime
101+
ds_ww3_cartesian.attrs["source"] = "ww3"
102+
# TODO: check kx ky names to be same as the one from intra burst ds
103+
indiceww3spectra = find_closest_ww3(
104+
ww3_path=pathww3sp, lon=lonsar, lat=latsar, time=start_date_dt
105+
)
106+
# add the raw EFTH(f,dir) spectra from WW3
107+
rawspww3 = (
108+
dsww3raw["efth"].isel(time=indiceww3spectra).rename("ww3EFTHraw")
109+
)
110+
rawspww3.attrs["description"] = "raw EFTH(f,dir) spectra"
111+
ds_ww3_cartesian = ds_ww3_cartesian.swap_dims(
112+
{"kx": "k_rg", "ky": "k_az"}
113+
).T
114+
rawspww3 = rawspww3.assign_coords(i)
115+
rawspww3 = rawspww3.expand_dims(["burst", "tile_line", "tile_sample"])
116+
ds_ww3_cartesian = ds_ww3_cartesian.assign_coords(i)
117+
ds_ww3_cartesian = ds_ww3_cartesian.expand_dims(
118+
["burst", "tile_line", "tile_sample"]
119+
)
120+
list_ww3_cart_sp.append(ds_ww3_cartesian)
121+
list_ww3_efth_sp.append(rawspww3)
122+
ds_ww3_cartesian_merged = xr.merge(list_ww3_cart_sp)
123+
ds_ww3_efth_merged = xr.merge(list_ww3_efth_sp)
124+
dsar = xr.merge([dsar, ds_ww3_cartesian_merged, ds_ww3_efth_merged])
125+
flag_ww3spectra_added = True
126+
if xspeckind == "intra":
127+
dsar = dsar.drop_vars(["xspectra_2tau_Re", "xspectra_2tau_Im"])
128+
elif xspeckind == "inter":
129+
dsar = dsar.drop_vars(["xspectra_Re", "xspectra_Im"])
130+
dsar = dsar.drop_dims(["freq_sample", "freq_line"])
131+
if xspeckind == "intra":
132+
dsar["xspectra_2tau"] = xsSAR
133+
elif xspeckind == "inter":
134+
dsar["xspectra"] = xsSAR
135+
else:
136+
raise ValueError("%s" % xspeckind)
135137
else:
136-
raise ValueError("%s" % xspeckind)
138+
logging.info('there is no xspectra in this subswath')
137139
return dsar, flag_ww3spectra_added, flag_ww3spectra_found
138140

139141

0 commit comments

Comments
 (0)