Skip to content

Commit f575d1c

Browse files
agrouazeagrouaze
authored
Merge pull request #91 from umr-lops/security-ww3spectra-colocs
Add WW3 spectra co-location security and debug information
2 parents 56708ea + 42f3877 commit f575d1c

File tree

8 files changed

+263
-35
lines changed

8 files changed

+263
-35
lines changed

‎assests/LOPS_WW3-GLOB-30M_202302_trck.nc

264 KB
Binary file not shown.

‎slcl1butils/coloc/coloc.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -76,6 +76,7 @@ def coloc_tiles_from_l1bgroup_with_raster(l1b_ds, raster_bb_ds, apply_merging=Tr
7676
projected_field = upscaled_da.interp(
7777
x=lonsar, y=latsar, assume_sorted=False
7878
).drop_vars(["x", "y"])
79+
projected_field.attrs['source'] = raster_bb_ds.attrs.get('name', 'unknown')
7980
mapped_ds_list.append(projected_field)
8081
raster_mapped = xr.merge(mapped_ds_list)
8182

‎slcl1butils/coloc/coloc_IW_WW3spectra.py

Lines changed: 174 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66
import numpy as np
77
import xarray as xr
88

9-
from slcl1butils.legacy_ocean.ocean.xPolarSpectrum import find_closest_ww3
9+
from slcl1butils.legacy_ocean.ocean.xPolarSpectrum import haversine
1010

1111
# from ocean.xspectrum import from_ww3
1212
# from ocean.xPolarSpectrum import find_closest_ww3
@@ -15,6 +15,125 @@
1515
from slcl1butils.symmetrize_l1b_spectra import symmetrize_xspectrum
1616
from slcl1butils.utils import xndindex
1717

18+
COLOC_LIMIT_SPACE = 100 # km
19+
COLOC_LIMIT_TIME = 3 # hours
20+
INDEX_WW3_FILL_VALUE = -999
21+
22+
23+
def check_colocation(cartesianspWW3, lon, lat, time):
24+
"""
25+
from information associated to the closest WW3 spectra found, check whether it is a usable co-location or not.
26+
and add co-location score information
27+
28+
Args:
29+
cartesianspWW3 (xr.DataArray): closest cartesian wave spectra WW3 found (could be too far for colocation)
30+
lon (float): longitude of interest (eg SAR)
31+
lat (float): latitude of interest (eg SAR)
32+
time (datetime.datetime or tuple of int): Tuple of the form (year, month, day, hour, minute, second) (eg SAR)
33+
34+
35+
Returns:
36+
(xr.DataArray): time index of spatio-temporal closest point in the WW3 file
37+
"""
38+
coloc_ds = xr.Dataset()
39+
mytime = (
40+
np.datetime64(time)
41+
if type(time) == datetime.datetime
42+
else np.datetime64(datetime.datetime(*time))
43+
)
44+
# time_dist = np.abs(cartesianspWW3.time - mytime)
45+
46+
# isel = np.where(time_dist == time_dist.min())
47+
closest_dist_in_space = haversine(
48+
lon,
49+
lat,
50+
cartesianspWW3.attrs["longitude"],
51+
cartesianspWW3.attrs["latitude"],
52+
)
53+
# spatial_dist = haversine(
54+
# lon,
55+
# lat,
56+
# ww3spds[{"time": isel[0]}].longitude,
57+
# ww3spds[{"time": isel[0]}].latitude,
58+
# )
59+
# logging.debug(spatial_dist)
60+
# relative_index = np.argmin(spatial_dist.data)
61+
# absolute_index = isel[0][relative_index]
62+
# closest_dist_in_space = spatial_dist[relative_index].data
63+
# time_dist_minutes = (ww3spec[{"time": absolute_index}].time - mytime).data / (1e9 * 60)
64+
# time_dist_minutes = (
65+
# ww3spds[{"time": absolute_index}].time - mytime
66+
# ) / np.timedelta64(1, "m")
67+
time_dist_minutes = (cartesianspWW3.time - mytime) / np.timedelta64(1, "m")
68+
logging.debug(
69+
"Wave Watch III closest point @ {} km and {} minutes".format(
70+
closest_dist_in_space, time_dist_minutes
71+
)
72+
)
73+
if (
74+
abs(time_dist_minutes) > COLOC_LIMIT_TIME * 60
75+
or closest_dist_in_space > COLOC_LIMIT_SPACE
76+
):
77+
logging.debug(
78+
"closest in time then in space is beyond the limits -> No ww3 spectra will be associated "
79+
)
80+
absolute_index = INDEX_WW3_FILL_VALUE
81+
ww3_lon = np.nan
82+
ww3_lat = np.nan
83+
ww3_date = np.nan
84+
selection = xr.DataArray(
85+
absolute_index,
86+
attrs={
87+
"method": "closest in time then in space",
88+
"limit_of_selection_in_space": "%f km" % COLOC_LIMIT_SPACE,
89+
"limit_of_selection_in_time": "%f hours" % COLOC_LIMIT_TIME,
90+
"description": "index of the WW3 spectra selected (None=no WW3 spectra found)",
91+
},
92+
)
93+
else:
94+
# ww3_lon = ww3spds[{"time": isel[0]}].longitude.data[relative_index]
95+
# ww3_lat = ww3spds[{"time": isel[0]}].latitude.data[relative_index]
96+
# ww3_date = ww3spds[{"time": absolute_index}].time
97+
ww3_lon = cartesianspWW3.attrs["longitude"]
98+
ww3_lat = cartesianspWW3.attrs["latitude"]
99+
ww3_date = cartesianspWW3.time
100+
selection = xr.DataArray(
101+
cartesianspWW3.attrs["time_index"],
102+
attrs={
103+
"method": "closest in time then in space",
104+
"limit_of_selection_in_space": "%f km" % COLOC_LIMIT_SPACE,
105+
"limit_of_selection_in_time": "%f hours" % COLOC_LIMIT_TIME,
106+
"description": "index of the WW3 spectra selected (None=no WW3 spectra found)",
107+
"_FillValue": INDEX_WW3_FILL_VALUE,
108+
},
109+
)
110+
coloc_ds["WW3spectra_index"] = selection
111+
coloc_ds["WW3spectra_delta_time"] = xr.DataArray(
112+
time_dist_minutes,
113+
attrs={"description": "temporal distance (WW3-SAR) in minutes", "units": "minutes"},
114+
)
115+
coloc_ds["WW3spectra_delta_space"] = xr.DataArray(
116+
closest_dist_in_space,
117+
attrs={
118+
"description": "spatial distance between SAR tile center and WW3 spectra grid point (km)",
119+
"units": "kilometer",
120+
},
121+
)
122+
coloc_ds["WW3spectra_longitude"] = xr.DataArray(
123+
ww3_lon, attrs={"description": "longitude of colocated WW3 spectra"}
124+
)
125+
coloc_ds["WW3spectra_latitude"] = xr.DataArray(
126+
ww3_lat, attrs={"description": "latitude of colocated WW3 spectra"}
127+
)
128+
coloc_ds["WW3spectra_time"] = xr.DataArray(
129+
ww3_date, attrs={"description": "time associated to colocated WW3 spectra"}
130+
)
131+
coloc_ds["WW3spectra_path"] = xr.DataArray(
132+
cartesianspWW3.attrs["pathWW3"],
133+
attrs={"description": "file path used to colocate WW3 spectra"},
134+
)
135+
return coloc_ds
136+
18137

19138
def resampleWW3spectra_on_TOPS_SAR_cartesian_grid(dsar, xspeckind):
20139
"""
@@ -67,6 +186,7 @@ def resampleWW3spectra_on_TOPS_SAR_cartesian_grid(dsar, xspeckind):
67186

68187
list_ww3_cart_sp = []
69188
list_ww3_efth_sp = []
189+
list_ww3_coloc_sp_ds = []
70190
dk_az = np.diff(xsSAR["k_az"])
71191
dk_rg = np.diff(xsSAR["k_rg"])
72192
dk = (dk_rg.mean(), dk_az.mean())
@@ -80,10 +200,12 @@ def resampleWW3spectra_on_TOPS_SAR_cartesian_grid(dsar, xspeckind):
80200
if pathww3sp:
81201
if os.path.exists(pathww3sp):
82202
flag_ww3spectra_found = True
203+
83204
dsww3raw = xr.open_dataset(pathww3sp)
84205
for i in xndindex(gridsar):
85206
lonsar = dsar["longitude"][i].values
86207
latsar = dsar["latitude"][i].values
208+
sensing_time = dsar["sensing_time"][i].values
87209
heading = dsar["ground_heading"][i].values
88210
logging.debug("heading %s", heading)
89211
logging.debug("timesar :%s", start_date_dt)
@@ -100,37 +222,61 @@ def resampleWW3spectra_on_TOPS_SAR_cartesian_grid(dsar, xspeckind):
100222
clockwise_to_trigo=True,
101223
lon=lonsar,
102224
lat=latsar,
103-
time=start_date_dt,
104-
) # TODO use sensingTime
105-
ds_ww3_cartesian.attrs["source"] = "ww3"
106-
# TODO: check kx ky names to be same as the one from intra burst ds
107-
indiceww3spectra = find_closest_ww3(
108-
ww3_path=pathww3sp, lon=lonsar, lat=latsar, time=start_date_dt
109-
)
110-
# add the raw EFTH(f,dir) spectra from WW3
111-
rawspww3 = (
112-
dsww3raw["efth"]
113-
.isel(time=indiceww3spectra)
114-
.rename("ww3EFTHraw")
225+
time=sensing_time,
115226
)
116-
rawspww3.attrs["description"] = "raw EFTH(f,dir) spectra"
117-
ds_ww3_cartesian = ds_ww3_cartesian.swap_dims(
118-
{"kx": "k_rg", "ky": "k_az"}
119-
).T
120-
rawspww3 = rawspww3.assign_coords(i)
121-
rawspww3 = rawspww3.expand_dims(["tile_line", "tile_sample"])
122-
# rawspww3 = rawspww3.expand_dims(["burst", "tile_line", "tile_sample"])
123-
ds_ww3_cartesian = ds_ww3_cartesian.assign_coords(i)
124-
ds_ww3_cartesian = ds_ww3_cartesian.expand_dims(
125-
["tile_line", "tile_sample"]
227+
ds_ww3_cartesian.attrs["source"] = "ww3"
228+
ds_ww3_cartesian.attrs[
229+
"description"
230+
] = "WW3spectra_EFTHraw resampled on SAR cartesian grid"
231+
ds_ww3_cartesian = ds_ww3_cartesian.rename("WW3spectra_EFTHcart")
232+
colocww3sp_ds = check_colocation(
233+
cartesianspWW3=ds_ww3_cartesian,
234+
lon=lonsar,
235+
lat=latsar,
236+
time=start_date_dt,
126237
)
127-
# ds_ww3_cartesian = ds_ww3_cartesian.expand_dims(["burst", "tile_line", "tile_sample"])
128-
list_ww3_cart_sp.append(ds_ww3_cartesian)
129-
list_ww3_efth_sp.append(rawspww3)
238+
del ds_ww3_cartesian.attrs["longitude"]
239+
del ds_ww3_cartesian.attrs["latitude"]
240+
if colocww3sp_ds["WW3spectra_index"].data != INDEX_WW3_FILL_VALUE:
241+
# add the raw EFTH(f,dir) spectra from WW3
242+
rawspww3 = (
243+
dsww3raw["efth"]
244+
.isel(time=colocww3sp_ds["WW3spectra_index"].data)
245+
.rename("WW3spectra_EFTHraw")
246+
)
247+
rawspww3.attrs[
248+
"description"
249+
] = "colocated raw EFTH(f,dir) WAVEWATCH III wave height spectra"
250+
colocww3sp_ds = colocww3sp_ds.assign_coords(i)
251+
colocww3sp_ds = colocww3sp_ds.expand_dims(
252+
["tile_line", "tile_sample"]
253+
)
254+
list_ww3_coloc_sp_ds.append(colocww3sp_ds)
255+
256+
ds_ww3_cartesian = ds_ww3_cartesian.swap_dims(
257+
{"kx": "k_rg", "ky": "k_az"}
258+
).T
259+
rawspww3 = rawspww3.assign_coords(i)
260+
rawspww3 = rawspww3.expand_dims(["tile_line", "tile_sample"])
261+
ds_ww3_cartesian = ds_ww3_cartesian.assign_coords(i)
262+
ds_ww3_cartesian = ds_ww3_cartesian.expand_dims(
263+
["tile_line", "tile_sample"]
264+
)
265+
list_ww3_cart_sp.append(ds_ww3_cartesian)
266+
list_ww3_efth_sp.append(rawspww3)
267+
flag_ww3spectra_added = True
130268
ds_ww3_cartesian_merged = xr.merge(list_ww3_cart_sp)
131269
ds_ww3_efth_merged = xr.merge(list_ww3_efth_sp)
132-
dsar = xr.merge([dsar, ds_ww3_cartesian_merged, ds_ww3_efth_merged])
133-
flag_ww3spectra_added = True
270+
ds_ww3_coloc_merged = xr.merge(list_ww3_coloc_sp_ds)
271+
272+
dsar = xr.merge(
273+
[
274+
dsar,
275+
ds_ww3_cartesian_merged,
276+
ds_ww3_efth_merged,
277+
ds_ww3_coloc_merged,
278+
]
279+
)
134280
if xspeckind == "intra":
135281
dsar = dsar.drop_vars(["xspectra_2tau_Re", "xspectra_2tau_Im"])
136282
elif xspeckind == "inter":

‎slcl1butils/legacy_ocean/ocean/xPolarSpectrum.py

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -294,7 +294,9 @@ def from_ww3(ds, *, station_index=None, time_index=None, rotate=0., shortWavesFi
294294
myspec = myspec.assign_coords(dk=dk)
295295
# myspec.attrs.update({'dk':dk, 'dphi':dphi})
296296
myspec.attrs.update({'dphi':dphi})
297-
myspec.attrs.update({'longitude':specWW3[ww3_kwargs]['longitude'].data, 'latitude':specWW3[ww3_kwargs]['latitude'].data})
297+
file_path = specWW3.encoding.get('source', None)
298+
myspec.attrs.update({'pathWW3':file_path})
299+
myspec.attrs.update({'longitude':specWW3[ww3_kwargs]['longitude'].data, 'latitude':specWW3[ww3_kwargs]['latitude'].data,'time_index':time_index})
298300

299301
try:
300302
myspec.attrs.update({'wd':np.radians(specWW3[ww3_kwargs].wnddir.data+180), 'ws':specWW3[ww3_kwargs].wnd.data})# wnddir (in WW3 files) is the direction the wind is coming from
@@ -507,7 +509,12 @@ def find_closest_ww3(ww3_path, lon, lat, time):
507509
import numpy as np
508510
import xarray as xr
509511
ww3spec = xr.open_dataset(ww3_path)
510-
mytime = np.datetime64(time) if type(time)==datetime else np.datetime64(datetime(*time))
512+
if type(time)==datetime:
513+
mytime = np.datetime64(time)
514+
elif isinstance(time,np.datetime64):
515+
mytime = time
516+
else:
517+
mytime = np.datetime64(datetime(*time))
511518
time_dist = np.abs(ww3spec.time-mytime)
512519
isel = np.where(time_dist==time_dist.min())
513520
spatial_dist = haversine(lon, lat, ww3spec[{'time':isel[0]}].longitude, ww3spec[{'time':isel[0]}].latitude)

‎slcl1butils/legacy_ocean/ocean/xspectrum.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -311,7 +311,7 @@ def from_xPolarSpectrum(ds, **kwargs):
311311
myspec.attrs.update({'dkx':dkx, 'dky':dky})
312312
myspec = myspec.where(np.sqrt(myspec.kx**2+myspec.ky**2)>min(ds.k).data, other=0.)# spectrum at wavenumbers smaller than the minimum one of the polar spectrum are set to zero
313313
myspec.attrs.pop('dphi', None)
314-
return myspec.drop(('k','phi','dk'))
314+
return myspec.drop_vars(('k','phi','dk'))
315315

316316

317317
# ******************************************************************************************************************************************************************

‎slcl1butils/scripts/do_IW_L1C_SAFE_from_L1B_SAFE.py

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
import pdb
12
import argparse
23
import logging
34
import os
@@ -284,7 +285,7 @@ def append_ancillary_field(ancillary, ds_intra, ds_inter):
284285
raster_bb_ds = raster_cropping_in_polygon_bounding_box(
285286
polygons["swath"][0], raster_ds
286287
)
287-
288+
raster_bb_ds.attrs['name'] = ancillary['name']
288289
# Loop on the grid in the product
289290
burst_types = ["intra", "inter"]
290291
for burst_type in burst_types:
@@ -300,6 +301,7 @@ def append_ancillary_field(ancillary, ds_intra, ds_inter):
300301
# ds_intra_list.append(_ds_intra)
301302
# Merging the datasets
302303
ds_intra = xr.merge([ds_intra, _ds_intra])
304+
ds_intra.attrs[ancillary["name"] + "_pattern"] = filename
303305
else:
304306
# l1b_ds_inter = xr.open_dataset(_file,group=burst_type+'burst')
305307
# _ds = coloc_tiles_from_l1bgroup_with_raster(l1b_ds_inter, raster_bb_ds, apply_merging=False)
@@ -310,6 +312,7 @@ def append_ancillary_field(ancillary, ds_intra, ds_inter):
310312
# ds_inter_list.append(_ds_inter)
311313
# Merging the datasets
312314
ds_inter = xr.merge([ds_inter, _ds_inter])
315+
ds_inter.attrs[ancillary["name"] + "_pattern"] = filename
313316
logging.debug("ancillary fields added")
314317
return ds_intra, ds_inter, ancillary_fields_added
315318

‎slcl1butils/tests/test_find_closest_ww3sp.py

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
from slcl1butils.coloc.coloc_IW_WW3spectra import check_colocation
2+
from slcl1butils.legacy_ocean.ocean.xspectrum import from_ww3
3+
import pytest
4+
import os
5+
import slcl1butils
6+
import datetime
7+
from slcl1butils.get_config import get_conf
8+
conf = get_conf()
9+
ww3spectra_hindcast_file = os.path.abspath(os.path.join(os.path.dirname(slcl1butils.__file__),conf['data_dir'],'LOPS_WW3-GLOB-30M_202302_trck.nc'))
10+
cases = {
11+
'matching':{
12+
#matching both in time and space:
13+
'lon':177.5,
14+
'lat':-34,
15+
'datesar': datetime.datetime(2023,2,13,17,30),
16+
'expected': 56
17+
},
18+
'no_match_in_space':{
19+
#not matching in space
20+
'lon': -50,
21+
'lat': 40,
22+
'datesar': datetime.datetime(2023,2,13,17,30),
23+
'expected':-999
24+
},
25+
'no_match_in_time':{
26+
#not matching in time
27+
'lon':177.5,
28+
'lat':-34,
29+
'datesar': datetime.datetime(2023,2,13,14,29),
30+
'expected': -999,
31+
}
32+
}
33+
params_list = []
34+
for case in cases:
35+
params_list.append((cases[case]['lon'],cases[case]['lat'],cases[case]['datesar'],cases[case]['expected']))
36+
37+
@pytest.mark.parametrize(
38+
["lon", "lat","datesar","expected"],
39+
tuple(params_list),
40+
)
41+
def test_finder(lon,lat,datesar,expected):
42+
dk = (0.001795349, 0.0018014804)
43+
kmax = (0.8114978075027466, 0.045037008821964264)
44+
rotate = 79.63082299999999
45+
ds_ww3_cartesian = from_ww3(
46+
ww3spectra_hindcast_file,
47+
dk=dk,
48+
kmax=kmax,
49+
strict="dk",
50+
rotate=rotate,
51+
clockwise_to_trigo=True,
52+
lon=lon,
53+
lat=lat,
54+
time=datesar,
55+
) # TODO use sensingTime
56+
ds_ww3_cartesian.attrs["source"] = "ww3"
57+
ds_ww3_cartesian.attrs[
58+
"description"
59+
] = "WW3spectra_EFTHraw resampled on SAR cartesian grid"
60+
ds_ww3_cartesian = ds_ww3_cartesian.rename("WW3spectra_EFTHcart")
61+
da_index = check_colocation(cartesianspWW3=ds_ww3_cartesian,lon=lon,lat=lat,time=datesar)
62+
actual_index = da_index['WW3spectra_index'].data
63+
# print('da_index',da_index)
64+
assert actual_index==expected
65+
66+
if __name__ == '__main__':
67+
for case in cases:
68+
test_finder(lon=cases[case]['lon'],lat=cases[case]['lat'],
69+
datesar=cases[case]['datesar'],expected=cases[case]['expected'])
70+
print('OK successful test')

0 commit comments

Comments
 (0)