This is an automated email from the ASF dual-hosted git repository. rkk pushed a commit to branch tmp-stv in repository https://gitbox.apache.org/repos/asf/sdap-nexus.git
commit aec77f5b2903b5bab05b025579a3c5dec288bea2 Author: rileykk <[email protected]> AuthorDate: Thu Jun 13 14:39:54 2024 -0700 Vertical truncation of tomograms by RH98/GND height maps --- analysis/webservice/algorithms/Tomogram.py | 143 +++++++++++++++++++++-- analysis/webservice/algorithms/Tomogram3D.py | 168 +++++++++++++++++++++------ 2 files changed, 267 insertions(+), 44 deletions(-) diff --git a/analysis/webservice/algorithms/Tomogram.py b/analysis/webservice/algorithms/Tomogram.py index 51ce091..e6220bd 100644 --- a/analysis/webservice/algorithms/Tomogram.py +++ b/analysis/webservice/algorithms/Tomogram.py @@ -173,7 +173,9 @@ class TomogramBaseClass(NexusCalcHandler): return binned_subset @staticmethod - def data_subset_to_ds_with_elevation(data_in_bounds: list) -> Tuple[np.ndarray, np.ndarray, np.ndarray, xr.Dataset]: + def data_subset_to_ds_with_elevation( + data_in_bounds: list + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, xr.Dataset]: lats = np.unique([d['latitude'] for d in data_in_bounds]) lons = np.unique([d['longitude'] for d in data_in_bounds]) elevations = np.unique([d['elevation'] for d in data_in_bounds]) @@ -200,6 +202,48 @@ class TomogramBaseClass(NexusCalcHandler): return lats, lons, vals, ds + @staticmethod + def add_variables_to_tomo(tomo_ds, **variables): + X, Y = np.meshgrid(tomo_ds.longitude.to_numpy(), tomo_ds.latitude.to_numpy()) + + for var in variables: + try: + xr.align(tomo_ds['tomo'], variables[var].tomo, join='exact') + tomo_ds[var] = variables[var].tomo + except: + from scipy.interpolate import griddata + + logger.info(f'Interpolating map to tomo grid') + + var_ds = variables[var] + data_points = [] + + for lon_i, lon in enumerate(var_ds.longitude): + for lat_i, lat in enumerate(var_ds.latitude): + val = var_ds.isel(latitude=lat_i, longitude=lon_i).tomo.item() + + data_points.append((lon, lat, val)) + + regridded = griddata( + list(zip( + [p[0] for p in data_points], [p[1] for p in data_points] + )), + np.array([p[2] for p in data_points]), + (X, Y), + method='nearest', + fill_value=np.nan + ) + + tomo_ds[var] = xr.DataArray( + data=regridded, + dims=['latitude', 'longitude'], + coords=dict( + latitude=tomo_ds.latitude, + longitude=tomo_ds.longitude + ), + name=var + ) + @nexus_handler class ElevationTomogramImpl(TomogramBaseClass): @@ -342,7 +386,19 @@ class LongitudeTomogramImpl(TomogramBaseClass): "name": "Calculate peaks", "type": "boolean", "description": "Calculate peaks along tomogram slice (currently uses simplest approach)" - } + }, + "canopy_ds": { + "name": "Canopy height dataset name", + "type": "string", + "description": "Dataset containing canopy heights (RH98). This is used to trim out tomogram voxels that " + "return over the measured or predicted canopy" + }, + "ground_ds": { + "name": "Ground height dataset name", + "type": "string", + "description": "Dataset containing ground height (DEM). This is used to trim out tomogram voxels that " + "return below the measured or predicted ground" + }, } singleton = True @@ -380,12 +436,15 @@ class LongitudeTomogramImpl(TomogramBaseClass): peaks = compute_options.get_boolean_arg('peaks') + ch_ds = compute_options.get_argument('canopy_ds', None) + g_ds = compute_options.get_argument('ground_ds', None) + return (ds, parameter, longitude, min_lat, max_lat, min_elevation, max_elevation, elevation_margin, - horizontal_margin, stride, peaks) + horizontal_margin, stride, peaks, ch_ds, g_ds) def calc(self, compute_options, **args): (dataset, parameter, longitude, min_lat, max_lat, min_elevation, max_elevation, elevation_margin, - horizontal_margin, stride, calc_peaks) = self.parse_args(compute_options) + horizontal_margin, stride, calc_peaks, ch_ds, g_ds) = self.parse_args(compute_options) slices = dict( lat=slice(min_lat, max_lat), @@ -393,6 +452,10 @@ class LongitudeTomogramImpl(TomogramBaseClass): elevation=slice(min_elevation, max_elevation) ) + # TODO: Do stride/elevation margin need to be provided or can they be computed? This is important because + # getting them wrong will mess up the result. (?): Stride === elevation diff between adjacent layers; + # margin === stride / 2 + r = np.arange(min_elevation, max_elevation + stride, stride) data_in_bounds = self.do_subset(dataset, parameter, slices, elevation_margin) @@ -403,8 +466,27 @@ class LongitudeTomogramImpl(TomogramBaseClass): TomogramBaseClass.bin_subset_elevations_to_range(data_in_bounds, r, elevation_margin) )[3] - if len(ds.longitude) > 1: - ds = ds.sel(longitude=longitude, method='nearest') + slices['elevation'] = slice(None, None) + elev_vars = {} + + if ch_ds is not None: + elev_vars['ch'] = TomogramBaseClass.data_subset_to_ds(self.do_subset(ch_ds, parameter, slices, 0))[3] + + if g_ds is not None: + elev_vars['gh'] = TomogramBaseClass.data_subset_to_ds(self.do_subset(g_ds, parameter, slices, 0))[3] + + # if len(ds.longitude) > 1: + # ds = ds.sel(longitude=longitude, method='nearest') + + if elev_vars: + TomogramBaseClass.add_variables_to_tomo(ds, **elev_vars) + + ds = ds.sel(longitude=longitude, method='nearest') + + if 'ch' in ds.data_vars: + ds = ds.where(ds.tomo.elevation <= ds.ch) + if 'gh' in ds.data_vars: + ds = ds.where(ds.tomo.elevation >= ds.gh) lats = ds.latitude.to_numpy() @@ -501,7 +583,19 @@ class LatitudeTomogramImpl(TomogramBaseClass): "name": "Calculate peaks", "type": "boolean", "description": "Calculate peaks along tomogram slice (currently uses simplest approach)" - } + }, + "canopy_ds": { + "name": "Canopy height dataset name", + "type": "string", + "description": "Dataset containing canopy heights (RH98). This is used to trim out tomogram voxels that " + "return over the measured or predicted canopy" + }, + "ground_ds": { + "name": "Ground height dataset name", + "type": "string", + "description": "Dataset containing ground height (DEM). This is used to trim out tomogram voxels that " + "return below the measured or predicted ground" + }, } singleton = True @@ -539,12 +633,15 @@ class LatitudeTomogramImpl(TomogramBaseClass): peaks = compute_options.get_boolean_arg('peaks') + ch_ds = compute_options.get_argument('canopy_ds', None) + g_ds = compute_options.get_argument('ground_ds', None)\ + return (ds, parameter, latitude, min_lon, max_lon, min_elevation, max_elevation, elevation_margin, - horizontal_margin, stride, peaks) + horizontal_margin, stride, peaks, ch_ds, g_ds) def calc(self, compute_options, **args): (dataset, parameter, latitude, min_lon, max_lon, min_elevation, max_elevation, elevation_margin, - horizontal_margin, stride, calc_peaks) = self.parse_args(compute_options) + horizontal_margin, stride, calc_peaks, ch_ds, g_ds) = self.parse_args(compute_options) slices = dict( lon=slice(min_lon, max_lon), @@ -562,8 +659,32 @@ class LatitudeTomogramImpl(TomogramBaseClass): TomogramBaseClass.bin_subset_elevations_to_range(data_in_bounds, r, elevation_margin) )[3] - if len(ds.latitude) > 1: - ds = ds.sel(latitude=latitude, method='nearest') + slices['elevation'] = slice(None, None) + elev_vars = {} + + if ch_ds is not None: + elev_vars['ch'] = TomogramBaseClass.data_subset_to_ds(self.do_subset(ch_ds, parameter, slices, 0))[3] + + if g_ds is not None: + elev_vars['gh'] = TomogramBaseClass.data_subset_to_ds(self.do_subset(g_ds, parameter, slices, 0))[3] + + try: + ds.to_netcdf('/tmp/tomo.nc') + except: + print('temp save failed :(') + + # if len(ds.latitude) > 1: + # ds = ds.sel(latitude=latitude, method='nearest') + + if elev_vars: + TomogramBaseClass.add_variables_to_tomo(ds, **elev_vars) + + ds = ds.sel(latitude=latitude, method='nearest') + + if 'ch' in ds.data_vars: + ds = ds.where(ds.tomo.elevation <= ds.ch) + if 'gh' in ds.data_vars: + ds = ds.where(ds.tomo.elevation >= ds.gh) lons = ds.longitude.to_numpy() diff --git a/analysis/webservice/algorithms/Tomogram3D.py b/analysis/webservice/algorithms/Tomogram3D.py index 3d15630..0214185 100644 --- a/analysis/webservice/algorithms/Tomogram3D.py +++ b/analysis/webservice/algorithms/Tomogram3D.py @@ -15,6 +15,7 @@ import contextlib import logging +import random from io import BytesIO from os.path import join from tempfile import TemporaryDirectory @@ -107,6 +108,18 @@ class Tomogram3D(NexusCalcHandler): "type": "boolean", "description": "If true, do not draw basemap beneath plot. This can be used to speed up render time a bit" }, + "canopy_ds": { + "name": "Canopy height dataset name", + "type": "string", + "description": "Dataset containing canopy heights (RH98). This is used to trim out tomogram voxels that " + "return over the measured or predicted canopy" + }, + "ground_ds": { + "name": "Ground height dataset name", + "type": "string", + "description": "Dataset containing ground height (DEM). This is used to trim out tomogram voxels that " + "return below the measured or predicted ground" + }, } singleton = True @@ -217,21 +230,24 @@ class Tomogram3D(NexusCalcHandler): hide_basemap = compute_options.get_boolean_arg('hideBasemap') + ch_ds = compute_options.get_argument('canopy_ds', None) + g_ds = compute_options.get_argument('ground_ds', None) + return (ds, parameter_s, bounding_poly, min_elevation, max_elevation, - (orbit_elev, orbit_step, frame_duration, view_azim, view_elev), filter_arg, render_type, hide_basemap) + (orbit_elev, orbit_step, frame_duration, view_azim, view_elev), filter_arg, render_type, hide_basemap, + ch_ds, g_ds) - def calc(self, computeOptions, **args): - (ds, parameter, bounding_poly, min_elevation, max_elevation, render_params, filter_arg, render_type, - hide_basemap) = (self.parse_args(computeOptions)) + def do_subset(self, ds, parameter, bounds): + tile_service = self._get_tile_service() - min_lat = bounding_poly.bounds[1] - max_lat = bounding_poly.bounds[3] - min_lon = bounding_poly.bounds[0] - max_lon = bounding_poly.bounds[2] + min_lat = bounds['lat'].start + max_lat = bounds['lat'].stop - bounds = (min_lat, min_lon, max_lat, max_lon) + min_lon = bounds['lon'].start + max_lon = bounds['lon'].stop - tile_service = self._get_tile_service() + min_elevation = bounds['elevation'].start + max_elevation = bounds['elevation'].stop tiles = tile_service.find_tiles_in_box( min_lat, max_lat, min_lon, max_lon, @@ -248,7 +264,7 @@ class Tomogram3D(NexusCalcHandler): data = [] - for i in range(len(tiles)-1, -1, -1): + for i in range(len(tiles) - 1, -1, -1): tile = tiles.pop(i) tile_id = tile.tile_id @@ -277,7 +293,8 @@ class Tomogram3D(NexusCalcHandler): break if data_val is None: - logger.warning(f'No variable {parameter} found at point {nexus_point.index} for tile {tile.tile_id}') + logger.warning( + f'No variable {parameter} found at point {nexus_point.index} for tile {tile.tile_id}') data_val = np.nan data.append({ @@ -287,32 +304,57 @@ class Tomogram3D(NexusCalcHandler): 'data': data_val }) - logger.info('Building dataframe') + return data + + def calc(self, computeOptions, **args): + (ds, parameter, bounding_poly, min_elevation, max_elevation, render_params, filter_arg, render_type, + hide_basemap, ch_ds, g_ds) = (self.parse_args(computeOptions)) + + min_lat = bounding_poly.bounds[1] + max_lat = bounding_poly.bounds[3] + min_lon = bounding_poly.bounds[0] + max_lon = bounding_poly.bounds[2] - lats = np.array([p['latitude'] for p in data]) - lons = np.array([p['longitude'] for p in data]) - elevs = np.array([p['elevation'] for p in data]) - - tomo = np.array([p['data'] for p in data]) - tomo = 10 * np.log10(tomo) - - df = pd.DataFrame( - np.hstack(( - lats[:, np.newaxis], - lons[:, np.newaxis], - elevs[:, np.newaxis], - tomo[:, np.newaxis] - )), - columns=[ - 'lat', - 'lon', - 'elevation', - 'tomo_value' - ] + bounds = dict( + lat=slice(min_lat, max_lat), + lon=slice(min_lon, max_lon), + elevation=slice(min_elevation, max_elevation) ) + data = self.do_subset(ds, parameter, bounds) + + logger.info('Building dataframe') + + if ch_ds is None and g_ds is None: + lats = np.array([p['latitude'] for p in data]) + lons = np.array([p['longitude'] for p in data]) + elevs = np.array([p['elevation'] for p in data]) + + tomo = np.array([p['data'] for p in data]) + tomo = 10 * np.log10(tomo) + + df = pd.DataFrame( + np.hstack(( + lats[:, np.newaxis], + lons[:, np.newaxis], + elevs[:, np.newaxis], + tomo[:, np.newaxis] + )), + columns=[ + 'lat', + 'lon', + 'elevation', + 'tomo_value' + ] + ) + else: + df = self.__mask_by_elev_map(data, ch_ds, g_ds, bounds.copy(), parameter) + df['tomo_value'] = 10 * np.log10(df['tomo_value']) + logger.info(f'DataFrame:\n{df}') + bounds = (bounds['lat'].start, bounds['lon'].start, bounds['lat'].stop, bounds['lon'].stop) + return Tomogram3DResults( df, render_params, @@ -322,6 +364,66 @@ class Tomogram3D(NexusCalcHandler): hide_basemap=hide_basemap, ) + def __mask_by_elev_map(self, data_points, ch_ds, g_ds, bounds, parameter): + from .Tomogram import TomogramBaseClass as Tomo + + # Sample data_points to determine difference in elevation between slices + + sample = data_points[random.randint(0, len(data_points) - 1)] + + sample = [d['elevation'] for d in data_points if + d['latitude'] == sample['latitude'] and d['longitude'] == sample['longitude']] + + sample.sort() + sample = [sample[i+1] - sample[i] for i in range(len(sample) - 1)] + + if min(sample) != max(sample): + raise NexusProcessingException( + code=500, + reason='There appears to be non-uniform elevation slices in the subset.' + ) + + stride = sample[0] + + # Bin and grid the subset data to an xarray Dataset + elevations = np.arange(bounds['elevation'].start, bounds['elevation'].stop + stride, stride/2) + bounds['elevation'] = slice(None, None) # Remove elevation from params as these are 2d datasets + + ds = Tomo.data_subset_to_ds_with_elevation( + Tomo.bin_subset_elevations_to_range( + data_points, elevations, stride/2 + ) + )[3] + + elev_vars = {} + + # Subset the elevation maps + if ch_ds is not None: + elev_vars['ch'] = Tomo.data_subset_to_ds(self.do_subset(ch_ds, parameter, bounds))[3] + + if g_ds is not None: + elev_vars['gh'] = Tomo.data_subset_to_ds(self.do_subset(g_ds, parameter, bounds))[3] + + # Add them to the Dataset, regridding if necessary + if elev_vars: + Tomo.add_variables_to_tomo(ds, **elev_vars) + + # Trim the data + if 'ch' in ds.data_vars: + ds = ds.where(ds.tomo.elevation <= ds.ch) + if 'gh' in ds.data_vars: + ds = ds.where(ds.tomo.elevation >= ds.gh) + + # Convert to DataFrame, move coords to columns, rename to be consistent with other approach and select cols + filtered_df = ds.to_dataframe().reset_index(['elevation', 'latitude', 'longitude']).rename( + columns=dict( + latitude='lat', longitude='lon', tomo='tomo_value' + ) + )[['lat', 'lon', 'elevation', 'tomo_value']] + + return filtered_df + # return filtered_df[pd.notna(filtered_df['tomo_value'])] + class Tomogram3DResults(NexusResults): def __init__(self, results=None, render_params=None, filter_args=None, render_type='scatter',
