# Licensed to the Apache Software Foundation (ASF) under one or more # contributor license agreements. See the NOTICE file distributed with # this work for additional information regarding copyright ownership. # The ASF licenses this file to You under the Apache License, Version 2.0 # (the "License"); you may not use this file except in compliance with # the License. You may obtain a copy of the License at # # http://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. import contextlib import logging import random from io import BytesIO from os.path import join from tempfile import TemporaryDirectory from urllib.parse import urlencode import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np import pandas as pd import requests from PIL import Image from matplotlib.colors import LinearSegmentedColormap, XKCD_COLORS from xarray import Dataset, apply_ufunc from webservice.NexusHandler import nexus_handler from webservice.algorithms.NexusCalcHandler import NexusCalcHandler from webservice.webmodel import NexusResults, NexusProcessingException, NoDataException logger = logging.getLogger(__name__) mpl.colormaps.register( LinearSegmentedColormap.from_list( 'forest', [ XKCD_COLORS['xkcd:dirt brown'], # (0.0, 1.0, 0.0), XKCD_COLORS['xkcd:dark forest green'], XKCD_COLORS['xkcd:forest green'], XKCD_COLORS['xkcd:forest'], XKCD_COLORS['xkcd:light forest green'], XKCD_COLORS['xkcd:leaf'], ]) ) @nexus_handler class Tomogram3D(NexusCalcHandler): name = '3D-Rendered Tomogram Subsetter' path = '/tomogram/3d' description = '3D visualization of tomogram subset, either as a static image or animated, orbiting GIF' params = { "ds": { "name": "Dataset", "type": "string", "description": "The Dataset shortname to use in calculation. Required" }, "parameter": { "name": "Parameter", "type": "string", "description": "The parameter of interest." }, "b": { "name": "Bounding box", "type": "comma-delimited float", "description": "Minimum (Western) Longitude, Minimum (Southern) Latitude, " "Maximum (Eastern) Longitude, Maximum (Northern) Latitude. Required." }, "minElevation": { "name": "Minimum Elevation", "type": "float", "description": "Minimum Elevation. Required." }, "maxElevation": { "name": "Maximum Elevation", "type": "float", "description": "Maximum Elevation. Required." }, "output": { "name": "Output format", "type": "string", "description": "Desired output format. Must be one of \"PNG\", \"GIF\", or \"CSV\". Required." }, "orbit": { "name": "Orbit settings", "type": "comma-delimited pair of ints", "description": "If output==GIF, specifies the orbit to be used in the animation. Format: " "elevation angle,orbit step. Default: 30, 10; Ranges: [-180,180],[1,90]" }, "viewAngle": { "name": "Static view angle", "type": "comma-delimited pair of ints", "description": "If output==PNG, specifies the angle to be used for the render. Format: " "azimuth,elevation angle. Default: 30,45; Ranges: [0,359],[-180,180]" }, "frameDuration": { "name": "Frame duration", "type": "int", "description": "If output==GIF, specifies the duration of each frame in the animation in milliseconds. " "Default: 100; Range: >=100" }, "filter": { "name": "Voxel filter", "type": "comma-delimited pair of numbers", "description": "Pair of numbers min,max. If defined, will filter out all points in the tomogram whose value" " does not satisfy min <= v <= max. Can be left unbounded, eg: 'filter=-15,' would filter " "everything but points >= -15. Ignored with \'elevPercentiles\'" }, "renderType": { "name": "Render Type", "type": "string", "description": "Type of render: Must be either \"scatter\" or \"peak\". Scatter will plot every voxel in " "the tomogram, colored by value, peak will plot a surface of all peaks for each lat/lon " "grid point, flat colored and shaded. Ignored with \'elevPercentiles\'. Default: scatter" }, "hideBasemap": { "name": "Hide Basemap", "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" }, "elevPercentiles": { "name": "Elevation percentiles", "type": "boolean", "description": "Display percentiles of cumulative returns over elevation. Requires both canopy_ds and " "ground_ds be set." }, "renderRH": { "name": "Render RH", "type": "int", "description": "With \'elevPercentiles\', this parameter is an integer between 0 and 100, which if set, " "will filter out all data whose cumulative return values is greater than the parameter " "value. It will also apply the color mapping to the voxel elevation rather than value, " "rendering a surface of that RH value." }, "cmap": { "name": "Color Map", "type": "string", "description": f"Color map to use. Will default to viridis if not provided or an unsupported map is " f"provided. Supported cmaps: {[k for k in sorted(mpl.colormaps.keys())]}" }, } singleton = True def __init__( self, tile_service_factory, **kwargs ): NexusCalcHandler.__init__(self, tile_service_factory) def parse_args(self, compute_options): try: ds = compute_options.get_dataset()[0] except: raise NexusProcessingException(reason="'ds' argument is required", code=400) parameter_s = compute_options.get_argument('parameter', None) try: bounding_poly = compute_options.get_bounding_polygon() except: raise NexusProcessingException(reason='Missing required parameter: b', code=400) def get_required_float(name): val = compute_options.get_float_arg(name, None) if val is None: raise NexusProcessingException(reason=f'Missing required parameter: {name}', code=400) return val min_elevation = get_required_float('minElevation') max_elevation = get_required_float('maxElevation') output = compute_options.get_argument('output', None) if output is not None: output = output.upper() if output not in ['PNG', 'GIF', 'CSV', 'NETCDF']: raise NexusProcessingException(reason=f'Missing or invalid required parameter: output = {output}', code=400) orbit_params = compute_options.get_argument('orbit', '30,10') view_params = compute_options.get_argument('viewAngle', '30,45') frame_duration = compute_options.get_int_arg('frameDuration', 100) if output == 'GIF': try: orbit_params = orbit_params.split(',') assert len(orbit_params) == 2 orbit_elev, orbit_step = tuple([int(p) for p in orbit_params]) assert -180 <= orbit_elev <= 180 assert 1 <= orbit_step <= 90 assert frame_duration >= 100 except: raise NexusProcessingException( reason=f'Invalid orbit parameters: {orbit_params} & {frame_duration}', code=400 ) view_azim, view_elev = None, None else: try: view_params = view_params.split(',') assert len(view_params) == 2 view_azim, view_elev = tuple([int(p) for p in view_params]) assert 0 <= view_azim <= 359 assert -180 <= view_elev <= 180 except: raise NexusProcessingException(reason=f'Invalid view angle string: {orbit_params}', code=400) orbit_elev, orbit_step = None, None ch_ds = compute_options.get_argument('canopy_ds', None) g_ds = compute_options.get_argument('ground_ds', None) percentiles = compute_options.get_boolean_arg('elevPercentiles') if percentiles and (ch_ds is None or g_ds is None): raise NexusProcessingException( code=400, reason='\'elevPercentiles\' argument requires \'canopy_ds\' and \'ground_ds\' be set' ) render_rh = compute_options.get_int_arg('renderRH', default=None) if render_rh is not None and not 0 <= render_rh <= 100: raise NexusProcessingException( code=400, reason="'renderRH' must be between 0 and 100" ) filter_arg = compute_options.get_argument('filter') if not percentiles else None if filter_arg is not None: try: filter_arg = filter_arg.split(',') assert len(filter_arg) == 2 vmin, vmax = filter_arg if vmin != '': vmin = float(vmin) else: vmin = None if vmax != '': vmax = float(vmax) else: vmax = None filter_arg = (vmin, vmax) except: raise NexusProcessingException( reason='Invalid filter arg, must be of format [number],[number]', code=400 ) else: filter_arg = (None, None) render_type = compute_options.get_argument('renderType', 'scatter').lower() if not percentiles else 'scatter' if render_type not in ['scatter', 'peak']: raise NexusProcessingException( reason='renderType must be either scatter or peak', code=400 ) hide_basemap = compute_options.get_boolean_arg('hideBasemap') cmap = mpl.colormaps.get( compute_options.get_argument('cmap', 'viridis'), mpl.colormaps['viridis'] ) 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, ch_ds, g_ds, percentiles, output, render_rh, cmap) def do_subset(self, ds, parameter, bounds): tile_service = self._get_tile_service() min_lat = bounds['lat'].start max_lat = bounds['lat'].stop min_lon = bounds['lon'].start max_lon = bounds['lon'].stop 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, ds=ds, min_elevation=min_elevation, max_elevation=max_elevation, fetch_data=False ) logger.info(f'Matched {len(tiles):,} tiles from Solr') if len(tiles) == 0: raise NoDataException(reason='No data was found within the selected parameters') data = [] for i in range(len(tiles) - 1, -1, -1): tile = tiles.pop(i) tile_id = tile.tile_id logger.info(f'Processing tile {tile_id} | {i=}') tile = tile_service.fetch_data_for_tiles(tile)[0] tile = tile_service.mask_tiles_to_bbox(min_lat, max_lat, min_lon, max_lon, [tile]) tile = tile_service.mask_tiles_to_elevation(min_elevation, max_elevation, tile) if len(tile) == 0: logger.info(f'Skipping empty tile {tile_id}') continue tile = tile[0] for nexus_point in tile.nexus_point_generator(): data_vals = nexus_point.data_vals if tile.is_multi else [nexus_point.data_vals] data_val = None for value, variable in zip(data_vals, tile.variables): if parameter is None or variable == parameter: data_val = value break if data_val is None: logger.warning( f'No variable {parameter} found at point {nexus_point.index} for tile {tile.tile_id}') data_val = np.nan data.append({ 'latitude': nexus_point.latitude, 'longitude': nexus_point.longitude, 'elevation': nexus_point.depth, 'data': data_val }) 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, percentiles, output, render_rh, cmap) = (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] 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 and output != 'NETCDF': 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, percentiles, output, render_rh is not None ) logger.info(f'Result repr:\n{df}') bounds = (bounds['lat'].start, bounds['lon'].start, bounds['lat'].stop, bounds['lon'].stop) if percentiles: if render_rh is None or output == 'NETCDF': style = 'percentile' else: df = df[df['tomo_value'] > 1 - (render_rh / 100)] style = 'percentile_elevation' else: style = 'db' return Tomogram3DResults( df, render_params, filter_arg, render_type, bounds=bounds, hide_basemap=hide_basemap, style=style, cmap=cmap ) def __mask_by_elev_map(self, data_points, ch_ds, g_ds, bounds, parameter, percentiles, output, invert): from .Tomogram import TomogramBaseClass as Tomo tries = 5 # Sample data_points to determine difference in elevation between slices while tries >= 0: # There's a bug where some samples yield non-uniform slices for what DEFINITELY was uniform data at ingest # Trying to get around it by just resampling for now, but this will definitely need further investigation tries -= 1 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): break elif tries < 0: raise NexusProcessingException( code=500, reason=f'There appears to be non-uniform elevation slices in the subset: ' f'{sample}, {min(sample)}, {max(sample)}, {len(sample)}' ) else: logger.warning(f'Possible bad elevation sampling. Trying again. ' f'{min(sample)}, {max(sample)}, {len(sample)}') 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['tomo'] = ds['tomo'].where(ds.tomo.elevation <= ds.ch) if 'gh' in ds.data_vars: ds['tomo'] = ds['tomo'].where(ds.tomo.elevation >= ds.gh) if percentiles: tomo = ds['tomo'].cumsum(dim='elevation', skipna=True).where(ds['tomo'].notnull()) ds['tomo'] = tomo / tomo.max(dim='elevation', skipna=True) if invert: ds['tomo'] = 1 - ds['tomo'] else: ds['tomo'] = apply_ufunc(lambda a: 10 * np.log10(a), ds.tomo) if output == 'NETCDF': return ds # 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 class Tomogram3DResults(NexusResults): def __init__(self, results=None, render_params=None, filter_args=None, render_type='scatter', style='db', bounds=None, meta=None, stats=None, computeOptions=None, status_code=200, **args): NexusResults.__init__(self, results, meta, stats, computeOptions, status_code, **args) self.render_params = render_params self.bounds = bounds if filter_args is None: self.filter = (None, None) else: self.filter = filter_args self.render_type = render_type self.style = style self.hide_basemap = 'hide_basemap' in args and args['hide_basemap'] self.cmap = args.get('cmap', mpl.colormaps['viridis']) def results(self): r: pd.DataFrame = NexusResults.results(self) r = Tomogram3DResults.min_max_filter(r, *self.filter) return r def result_dataset(self): ds = NexusResults.results(self) if not isinstance(ds, Dataset): raise NexusProcessingException( code=500, reason=f'Invalid result type {type(ds)}' ) return ds def __common(self): xyz = self.results()[['lon', 'lat', 'elevation']].values fig = plt.figure(figsize=(10, 7)) return xyz, (fig, fig.add_subplot(111, projection='3d')) @staticmethod def min_max_filter(df: pd.DataFrame, vmin=None, vmax=None): if vmin is not None or vmax is not None: n_points = len(df) if vmin is not None: df = df[df['tomo_value'] >= vmin] if vmax is not None: df = df[df['tomo_value'] <= vmax] logger.info(f'Filtered data from {n_points:,} to {len(df):,} points') return df def toImage(self): _, _, _, view_azim, view_elev = self.render_params results = self.results() xyz = results[['lon', 'lat', 'elevation']].values fig = plt.figure(figsize=(10, 7)) ax = fig.add_subplot(111, projection='3d') ax.view_init(elev=view_elev, azim=view_azim) if not self.hide_basemap: min_lat, min_lon, max_lat, max_lon = self.bounds # Originally used mpl_toolkits.basemap.Basemap.aspect for this, but basemap installs an LGPL3 data package # so we can't have that in a release. Thankfully, manual computation ends up being quite simple aspect = (max_lat - min_lat) / (max_lon - min_lon) basemap_size = 512 params = dict( bbox=f'{min_lon},{min_lat},{max_lon},{max_lat}', bboxSR=4326, imageSR=4326, size=f'{basemap_size},{int(aspect * basemap_size)}', dpi=2000, format='png32', transparent=True, f='image' ) url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/export' logger.info(f'Pulling basemap ({url}?{urlencode(params)})') try: elevations = results[['elevation']].values z_range = np.nanmax(elevations) - np.nanmin(elevations) z_coord = np.nanmin(elevations) - (0.1 * z_range) r = requests.get(url, params=params) r.raise_for_status() buf = BytesIO(r.content) img = Image.open(buf) img_data = np.flipud(np.array(img)) lats = np.linspace(min_lat, max_lat, num=img.height) lons = np.linspace(min_lon, max_lon, num=img.width) X, Y = np.meshgrid(lons, lats) Z = np.full(X.shape, z_coord) logger.info('Plotting basemap') ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=img_data / 255) except: logger.error('Failed to pull basemap, will not draw it') if self.render_type == 'scatter': logger.info('Plotting tomogram data') if self.style == 'db': v = dict(vmin=-30, vmax=-10) cbl = 'Tomogram (dB)' elif self.style == 'percentile': v = dict(vmin=0, vmax=1) cbl = 'Tomogram Cumulative value (%/100)' else: elevs = results[['elevation']].values v = dict(vmin=np.nanmin(elevs), vmax=np.nanmax(elevs)) cbl = 'Voxel elevation' s = ax.scatter( xyz[:, 0], xyz[:, 1], xyz[:, 2], marker=',', c=results[['tomo_value']].values if self.style != 'percentile_elevation' else results[['elevation']].values, zdir='z', depthshade=True, cmap=self.cmap, **v ) cbar = fig.colorbar(s, ax=ax) cbar.set_label(cbl) else: logger.info('Collecting data by lat/lon') lats = np.unique(results['lat'].values) lons = np.unique(results['lon'].values) data_dict = {} for r in results.itertuples(index=False): key = (r.lon, r.lat) if key not in data_dict: data_dict[key] = ([r.elevation], [r.tomo_value]) else: data_dict[key][0].append(r.elevation) data_dict[key][1].append(r.tomo_value) logger.info('Determining peaks') vals = np.empty((len(lats), len(lons))) for i, lat in enumerate(lats): for j, lon in enumerate(lons): if (lon, lat) in data_dict: elevs, tomo_vals = data_dict[(lon, lat)] i_max = np.argmax(np.array(tomo_vals)) vals[i, j] = elevs[i_max] else: vals[i, j] = np.nan X2, Y2 = np.meshgrid(lons, lats) logger.info('Plotting peak surface') s = ax.plot_surface( X2, Y2, vals, rstride=1, cstride=1, color='xkcd:leaf' ) ax.set_ylabel('Latitude') ax.set_xlabel('Longitude') ax.set_zlabel('Elevation w.r.t. dataset reference (m)') plt.tight_layout() buffer = BytesIO() logger.info('Writing plot to buffer') plt.savefig(buffer, format='png', facecolor='white') buffer.seek(0) plt.close(fig) return buffer.read() def toGif(self): orbit_elev, orbit_step, frame_duration, _, _ = self.render_params results = self.results() xyz = results[['lon', 'lat', 'elevation']].values fig = plt.figure(figsize=(10, 7)) ax = fig.add_subplot(111, projection='3d') ax.view_init(elev=orbit_elev, azim=0) if not self.hide_basemap: min_lat, min_lon, max_lat, max_lon = self.bounds # Originally used mpl_toolkits.basemap.Basemap.aspect for this, but basemap installs an LGPL3 data package # so we can't have that in a release. Thankfully, manual computation ends up being quite simple aspect = (max_lat - min_lat) / (max_lon - min_lon) basemap_size = 512 params = dict( bbox=f'{min_lon},{min_lat},{max_lon},{max_lat}', bboxSR=4326, imageSR=4326, size=f'{basemap_size},{int(aspect * basemap_size)}', dpi=2000, format='png32', transparent=True, f='image' ) url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/export' logger.info(f'Pulling basemap ({url}?{urlencode(params)})') try: elevations = results[['elevation']].values z_range = np.nanmax(elevations) - np.nanmin(elevations) z_coord = np.nanmin(elevations) - (0.1 * z_range) r = requests.get(url, params=params) r.raise_for_status() buf = BytesIO(r.content) img = Image.open(buf) img_data = np.flipud(np.array(img)) lats = np.linspace(min_lat, max_lat, num=img.height) lons = np.linspace(min_lon, max_lon, num=img.width) X, Y = np.meshgrid(lons, lats) Z = np.full(X.shape, z_coord) logger.info('Plotting basemap') ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=img_data / 255) except: logger.error('Failed to pull basemap, will not draw it') if self.render_type == 'scatter': logger.info('Plotting tomogram data') if self.style == 'db': v = dict(vmin=-30, vmax=-10) cbl = 'Tomogram (dB)' elif self.style == 'percentile': v = dict(vmin=0, vmax=1) cbl = 'Tomogram Cumulative value (%/100)' else: elevs = results[['elevation']].values v = dict(vmin=np.nanmin(elevs), vmax=np.nanmax(elevs)) cbl = 'Voxel elevation' s = ax.scatter( xyz[:, 0], xyz[:, 1], xyz[:, 2], marker=',', c=results[['tomo_value']].values if self.style != 'percentile_elevation' else results[['elevation']].values, zdir='z', depthshade=True, cmap=self.cmap, **v ) cbar = fig.colorbar(s, ax=ax) cbar.set_label(cbl) else: logger.info('Collecting data by lat/lon') lats = np.unique(results['lat'].values) lons = np.unique(results['lon'].values) data_dict = {} for r in results.itertuples(index=False): key = (r.lon, r.lat) if key not in data_dict: data_dict[key] = ([r.elevation], [r.tomo_value]) else: data_dict[key][0].append(r.elevation) data_dict[key][1].append(r.tomo_value) logger.info('Determining peaks') vals = np.empty((len(lats), len(lons))) for i, lat in enumerate(lats): for j, lon in enumerate(lons): if (lon, lat) in data_dict: elevs, tomo_vals = data_dict[(lon, lat)] i_max = np.argmax(np.array(tomo_vals)) vals[i, j] = elevs[i_max] else: vals[i, j] = np.nan X2, Y2 = np.meshgrid(lons, lats) logger.info('Plotting peak surface') ax.plot_surface( X2, Y2, vals, rstride=1, cstride=1, color='xkcd:leaf' ) ax.set_ylabel('Latitude') ax.set_xlabel('Longitude') ax.set_zlabel('Elevation w.r.t. dataset reference (m)') plt.tight_layout() buffer = BytesIO() with TemporaryDirectory() as td: for azim in range(0, 360, orbit_step): logger.info(f'Saving frame for azimuth = {azim}') ax.view_init(azim=azim) plt.savefig(join(td, f'fr_{azim}.png')) with contextlib.ExitStack() as stack: logger.info(f'Combining frames into final GIF') imgs = (stack.enter_context(Image.open(join(td, f'fr_{a}.png'))) for a in range(0, 360, orbit_step)) img = next(imgs) img.save(buffer, format='GIF', append_images=imgs, save_all=True, duration=frame_duration, loop=0) buffer.seek(0) plt.close(fig) return buffer.read() def toCSV(self): df = self.results()[['lon', 'lat', 'elevation', 'tomo_value']] buf = BytesIO() df.to_csv(buf, header=False, index=False) buf.seek(0) return buf.read() def toNetCDF(self): ds = self.result_dataset() comp = dict(zlib=True, complevel=9) encoding = {vname: comp for vname in ds.data_vars} with TemporaryDirectory() as td: ds.to_netcdf(join(td, 'ds.nc'), encoding=encoding) with open(join(td, 'ds.nc'), 'rb') as fp: buf = BytesIO(fp.read()) buf.seek(0) return buf.read()