analysis/webservice/algorithms/Tomogram3D.py

# 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()