in analysis/webservice/algorithms/Tomogram.py [0:0]
def calc(self, compute_options, **args):
(dataset, parameter, start_point, end_point, min_elevation, max_elevation, horizontal_margin, calc_peaks,
ch_ds, g_ds, percentiles, cmap) = self.parse_args(compute_options)
bbox_miny = min(start_point.y, end_point.y) - horizontal_margin
bbox_maxy = max(start_point.y, end_point.y) - horizontal_margin
bbox_minx = min(start_point.x, end_point.x) + horizontal_margin
bbox_maxx = max(start_point.x, end_point.x) + horizontal_margin
slices = dict(
lat=slice(bbox_miny, bbox_maxy),
lon=slice(bbox_minx, bbox_maxx),
elevation=slice(min_elevation, max_elevation)
)
try:
data_in_bounds = self.__subset_on_line(dataset, start_point, end_point, min_elevation, max_elevation,
parameter, horizontal_margin)
except NexusTileServiceException as e:
logger.warning('Selected dataset does not support subsetting along a line, will fall back to box subset')
data_in_bounds = self.do_subset(dataset, parameter, slices, 0)
except Exception as e:
logger.error('An unexpected exception occurred', exc_info=e)
raise NexusProcessingException(
code=500,
reason=str(e)
)
z_step = TomogramBaseClass._get_elevation_step_size(data_in_bounds)
r = np.arange(min_elevation, max_elevation + z_step, z_step)
logger.info(f'Fetched {len(data_in_bounds):,} data points at this elevation')
ds = TomogramBaseClass.data_subset_to_ds_with_elevation(
TomogramBaseClass.bin_subset_elevations_to_range(data_in_bounds, r, z_step / 2)
)[3]
slices['elevation'] = slice(None, None)
elev_vars = {}
if len(ch_ds) > 0:
elev_vars['ch'] = TomogramBaseClass.data_subset_to_ds(
self.do_subset(ch_ds[0], parameter, slices, 0)
)[3]
elev_vars['ch'].attrs['_source'] = ch_ds[0]
if len(ch_ds) == 2:
elev_vars['ch_secondary'] = TomogramBaseClass.data_subset_to_ds(
self.do_subset(ch_ds[1], parameter, slices, 0)
)[3]
elev_vars['ch_secondary'].attrs['_source'] = ch_ds[1]
# TODO: Special checks for NoDataException around these (primary map gets in all Impls) to reraise with better error message
if len(g_ds) > 0:
elev_vars['gh'] = TomogramBaseClass.data_subset_to_ds(
self.do_subset(g_ds[0], parameter, slices, 0)
)[3]
elev_vars['gh'].attrs['_source'] = g_ds[0]
if len(g_ds) == 2:
elev_vars['gh_secondary'] = TomogramBaseClass.data_subset_to_ds(
self.do_subset(g_ds[1], parameter, slices, 0)
)[3]
elev_vars['gh_secondary'].attrs['_source'] = g_ds[1]
if elev_vars:
TomogramBaseClass.add_variables_to_tomo(ds, **elev_vars)
if 'ch' in ds.data_vars:
if 'ch_secondary' in ds.data_vars:
ds['ch_filter'] = xr.where(np.isnan(ds.ch), ds.ch_secondary, ds.ch)
else:
ds['ch_filter'] = ds.ch
ds['tomo'] = ds['tomo'].where(ds.tomo.elevation <= ds.ch_filter)
if 'gh' in ds.data_vars:
if 'gh_secondary' in ds.data_vars:
ds['gh_filter'] = xr.where(np.isnan(ds.gh), ds.gh_secondary, ds.gh)
else:
ds['gh_filter'] = ds.gh
ds['tomo'] = ds['tomo'].where(ds.tomo.elevation >= ds.gh_filter)
if percentiles:
tomo = ds['tomo'].cumsum(dim='elevation', skipna=True).where(ds['tomo'].notnull())
ds['tomo'] = tomo / tomo.max(dim='elevation', skipna=True)
else:
ds['tomo'] = xr.apply_ufunc(lambda a: 10 * np.log10(a), ds.tomo)
line = LineString([start_point, end_point])
i_lat_1 = list(ds.latitude.values).index(ds.sel(latitude=start_point.y, method='nearest').latitude)
i_lat_2 = list(ds.latitude.values).index(ds.sel(latitude=end_point.y, method='nearest').latitude)
i_lon_1 = list(ds.longitude.values).index(ds.sel(longitude=start_point.x, method='nearest').longitude)
i_lon_2 = list(ds.longitude.values).index(ds.sel(longitude=end_point.x, method='nearest').longitude)
steps = ceil(sqrt((i_lat_2 - i_lat_1) ** 2 + (i_lon_2 - i_lon_1) ** 2))
point_coords = []
for p in np.linspace(0, 1, steps):
point_coords.append(line.interpolate(p, True).coords[0])
point_coords = [dict(longitude=p[0], latitude=p[1]) for p in point_coords]
points = [ds.sel(p, method='nearest') for p in point_coords]
ds = xr.concat(points, dim='distance')
def dist(lat, lon):
return abs(geodesic((start_point.y, start_point.x), (lat, lon)).m)
dists = []
for lat, lon in zip(ds.latitude.values, ds.longitude.values):
dists.append(dist(lat, lon))
ds = (ds.assign_coords(dict(distance=(['distance'], dists))).set_index(distance='distance').
expand_dims(dim='slice').drop_duplicates('distance'))
if calc_peaks:
peaks = []
for i in range(len(ds.tomo.slice)):
slice_peaks = []
ds_sub = ds.isel(slice=i)
for j in range(len(ds.tomo.distance)):
col = ds_sub.isel(distance=j)
distance = col.distance.item()
try:
idx = col.tomo.argmax(dim='elevation').item()
slice_peaks.append((distance, col.isel(elevation=idx).elevation.item()))
except ValueError:
slice_peaks.append((distance, np.nan))
peaks.append((np.array([p[0] for p in slice_peaks]),
np.array([p[1] for p in slice_peaks])))
else:
peaks = None
return ProfileTomoResults(
results=(ds, peaks),
s={'start_point': start_point, 'end_point': end_point},
coords=dict(x=ds.distance, y=ds.elevation),
meta=dict(dataset=dataset),
style='db' if not percentiles else 'percentile',
cmap=cmap,
computeOptions=compute_options
)