analysis/webservice/algorithms/Tomogram.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 os
import random
import re
import zipfile
from bisect import bisect_left
from io import BytesIO
from math import ceil, sqrt, isnan
from tempfile import TemporaryDirectory
from typing import Tuple

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from PIL import Image
from geopy.distance import geodesic
from shapely.geometry import Point, LineString

from nexustiles.exception import NexusTileServiceException
from webservice.NexusHandler import nexus_handler
from webservice.algorithms.NexusCalcHandler import NexusCalcHandler
from webservice.webmodel import NexusResults, NexusProcessingException, NoDataException

logger = logging.getLogger(__name__)
FLOAT_PATTERN = re.compile(r'^[+-]?(\d+(\.\d*)?|\.\d+)$')


class TomogramBaseClass(NexusCalcHandler):
    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)

        return ds, parameter_s

    def do_subset(self, ds, parameter, bounds, margin):
        tile_service = self._get_tile_service()

        # bounds = self.__slice_bounds

        if isinstance(bounds['lat'], slice):
            min_lat = bounds['lat'].start
            max_lat = bounds['lat'].stop
        else:
            min_lat = bounds['lat'] - margin
            max_lat = bounds['lat'] + margin

        if isinstance(bounds['lon'], slice):
            min_lon = bounds['lon'].start
            max_lon = bounds['lon'].stop
        else:
            min_lon = bounds['lon'] - margin
            max_lon = bounds['lon'] + margin

        if isinstance(bounds['elevation'], slice):
            min_elevation = bounds['elevation'].start
            max_elevation = bounds['elevation'].stop
        else:
            min_elevation = bounds['elevation'] - margin
            max_elevation = bounds['elevation'] + margin

        # In case slices are unbounded (ie, start/stop is None) replace the associated values with max/min for coord
        min_lat, min_lon, max_lat, max_lon = [v if v is not None else d for v, d in
                                              zip([min_lat, min_lon, max_lat, max_lon], [-90, -180, 90, 180])]

        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 tiles matched 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
                })

        if len(data) == 0:
            raise NoDataException(reason='No data fit the selected parameters')

        return data

    @staticmethod
    def data_subset_to_ds(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])

        vals = np.empty((len(lats), len(lons)))

        data_dict = {(d['latitude'], d['longitude']): d['data'] for d in data_in_bounds}

        for i, lat in enumerate(lats):
            for j, lon in enumerate(lons):
                vals[i, j] = data_dict.get((lat, lon), np.nan)

        ds = xr.Dataset(
            data_vars=dict(
                tomo=(('latitude', 'longitude'), vals)
            ),
            coords=dict(
                latitude=(['latitude'], lats),
                longitude=(['longitude'], lons)
            )
        )

        return lats, lons, vals, ds

    @staticmethod
    def bin_subset_elevations_to_range(data_in_bounds: list, elevation_range: np.ndarray, margin):
        binned_subset = []

        for d in data_in_bounds:
            elevation: float = d['elevation']
            d_copy = dict(**d)

            pos = bisect_left(elevation_range, elevation)

            if pos == 0:
                binned_elev = elevation_range[0]
            elif pos == len(elevation_range):
                binned_elev = elevation_range[-1]
            else:
                before = elevation_range[pos - 1]
                after = elevation_range[pos]
                if after - elevation < elevation - before:
                    binned_elev = after
                else:
                    binned_elev = before

            if abs(binned_elev - elevation) > margin:
                logger.warning(f'Could not bin point {d_copy}')
            else:
                d_copy['elevation'] = binned_elev
                binned_subset.append(d_copy)

        return binned_subset

    @staticmethod
    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])

        # Commented out is the old gridding method
        #
        # vals = np.empty((len(elevations), len(lats), len(lons)))
        #
        # time = elapsed('data_subset_to_ds_with_elevation: val init', time)
        #
        # data_dict = {(d['elevation'], d['latitude'], d['longitude']): d['data'] for d in data_in_bounds}
        #
        # time = elapsed('data_subset_to_ds_with_elevation: dict build', time)
        #
        # print(f'{len(elevations) * len(lats) * len(lons):,}')
        #
        # for i, elev in enumerate(elevations):
        #     for j, lat in enumerate(lats):
        #         for k, lon in enumerate(lons):
        #             vals[i, j, k] = data_dict.get((elev, lat, lon), np.nan)
        #
        # time = elapsed('data_subset_to_ds_with_elevation: populate', time)

        vals = np.full((len(elevations), len(lats), len(lons)), np.nan)

        elevation_idxs = {}
        lat_idxs = {}
        lon_idxs = {}

        for point in data_in_bounds:
            elevation = point['elevation']
            lat = point['latitude']
            lon = point['longitude']

            if elevation not in elevation_idxs:
                e_i = np.argwhere(elevations == elevation)
                elevation_idxs[elevation] = e_i
            else:
                e_i = elevation_idxs[elevation]

            if lat not in lat_idxs:
                lat_i = np.argwhere(lats == lat)
                lat_idxs[lat] = lat_i
            else:
                lat_i = lat_idxs[lat]

            if lon not in lon_idxs:
                lon_i = np.argwhere(lons == lon)
                lon_idxs[lon] = lon_i
            else:
                lon_i = lon_idxs[lon]

            vals[e_i, lat_i, lon_i] = point['data']

        ds = xr.Dataset(
            data_vars=dict(
                tomo=(('elevation', 'latitude', 'longitude'), vals)
            ),
            coords=dict(
                latitude=(['latitude'], lats),
                longitude=(['longitude'], lons),
                elevation=(['elevation'], elevations)
            )
        )

        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
                )

            tomo_ds[var].attrs.update(variables[var].attrs)

    @staticmethod
    def _get_elevation_step_size(data):
        tries = 5
        sample = [1]

        while tries >= 0:
            tries -= 1
            sample = data[random.randint(0, len(data) - 1)]

            sample = [d['elevation'] for d in data 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)}')

        return sample[0]


# @nexus_handler
class ElevationTomogramImpl(TomogramBaseClass):
    name = "Elevation-sliced Tomogram"
    path = "/tomogram/elevation"
    description = "Fetches a tomogram sliced by elevation"
    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."
        },
        "elevation": {
            "name": "Slice elevation",
            "type": "float",
            "description": "The desired elevation of the tomogram slice"
        },
        "margin": {
            "name": "Margin",
            "type": "float",
            "description": "Margin +/- desired elevation to include in output. Default: 0.5m"
        }
    }
    singleton = True

    def __init__(self, tile_service_factory, **kwargs):
        TomogramBaseClass.__init__(self, tile_service_factory, **kwargs)

    def parse_args(self, compute_options):
        try:
            bounding_poly = compute_options.get_bounding_polygon()
        except:
            raise NexusProcessingException(reason='Missing required parameter: b', code=400)

        elevation = compute_options.get_float_arg('elevation', None)

        if elevation is None:
            raise NexusProcessingException(reason='Missing required parameter: elevation', code=400)

        margin = compute_options.get_float_arg('margin', 0.5)

        ds, parameter = super().parse_args(compute_options)

        return ds, parameter, bounding_poly, elevation, margin

    def calc(self, compute_options, **args):
        dataset, parameter, bounding_poly, elevation, margin = self.parse_args(compute_options)

        min_lat = bounding_poly.bounds[1]
        max_lat = bounding_poly.bounds[3]
        min_lon = bounding_poly.bounds[0]
        max_lon = bounding_poly.bounds[2]

        slices = dict(
            lat=slice(min_lat, max_lat),
            lon=slice(min_lon, max_lon),
            elevation=float(elevation)
        )

        data_in_bounds = self.do_subset(dataset, parameter, slices, margin)
        logger.info(f'Fetched {len(data_in_bounds):,} data points, organizing them for plotting')

        lats, lons, vals, ds = TomogramBaseClass.data_subset_to_ds(data_in_bounds)

        result = ElevationTomoResults(
            (vals, lats, lons, dict(ds=dataset, elevation=elevation, margin=margin)),
        )

        return result


@nexus_handler
class LongitudeTomogramImpl(TomogramBaseClass):
    name = "Tomogram Latitude profile view"
    path = "/tomogram/longitude"
    description = "Fetches a tomogram sliced by longitude"
    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."
        },
        "longitude": {
            "name": "Longitude",
            "type": "float, comma-delimited floats(>1), colon-delimited floats(<=3), literal['*']",
            "description": "Longitude(s) at which to slice the tomogram. Required. Valid formats: Single float - Slice "
                           "at a singular longitude; comma-delimited floats - slice at a list of longitudes; "
                           "colon-delimited floats - slice using Python slicing notation (both ends are inclusive): "
                           "[<start>]:[<stop>]:[<step>] where omitting <start> and/or <stop> produce unbounded slices "
                           "and omitting <step> equates to an all-inclusive step size. A special value of \"*\" "
                           "equates to ::1. Values other than singular floats are incompatible with PNG rendering "
                           "(use ZIP to a ZIP archive of .pngs)"
        },
        "minLat": {
            "name": "Minimum Latitude",
            "type": "float",
            "description": "Minimum (Southern) Latitude. Required."
        },
        "maxLat": {
            "name": "Maximum Latitude",
            "type": "float",
            "description": "Maximum (Northern) Latitude. Required."
        },
        "minElevation": {
            "name": "Minimum Elevation",
            "type": "int",
            "description": "Minimum Elevation. Required."
        },
        "maxElevation": {
            "name": "Maximum Elevation",
            "type": "int",
            "description": "Maximum Elevation. Required."
        },
        "horizontalMargin": {
            "name": "Horizontal Margin",
            "type": "float",
            "description": "Margin +/- desired lat/lon slice to include in output. Default: 0.001m"
        },
        "peaks": {
            "name": "Calculate peaks",
            "type": "boolean",
            "description": "Calculate peaks along tomogram slice (currently uses simplest approach). "
                           "Ignored with \'elevPercentiles\'"
        },
        "canopy_ds": {
            "name": "Canopy height dataset name",
            "type": "string, comma-separated pair of strings",
            "description": "Dataset(s) containing canopy heights (RH98). This is used to trim out tomogram voxels that "
                           "return over the measured or predicted canopy. If 2 are provided, the first is considered "
                           "more accurate and will be favored for masking."
        },
        "ground_ds": {
            "name": "Ground height dataset name",
            "type": "string, comma-separated pair of strings",
            "description": "Dataset containing ground height (DEM). This is used to trim out tomogram voxels that "
                           "return below the measured or predicted ground. If 2 are provided, the first is considered "
                           "more accurate and will be favored for masking."
        },
        "elevPercentiles": {
            "name": "Elevation percentiles",
            "type": "boolean",
            "description": "Display percentiles of cumulative returns over elevation. Requires both canopy_ds and "
                           "ground_ds be set."
        },
        "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())]}"
        },
        "cbarMin": {
            "name": "Color bar min",
            "type": "float",
            "description": "Minimum value of the color bar. Only applies if elevPercentiles is False. Value must be "
                           "less than cbarMax. Invalid configuration will be replaced by default values (-30)"
        },
        "cbarMax": {
            "name": "Color bar max",
            "type": "float",
            "description": "Maximum value of the color bar. Only applies if elevPercentiles is False. Value must be "
                           " greater than cbarMin. Invalid configuration will be replaced by default value (-10)"
        },
    }
    singleton = True

    def __init__(self, tile_service_factory, **kwargs):
        TomogramBaseClass.__init__(self, tile_service_factory, **kwargs)

    def parse_args(self, compute_options):
        ds, parameter = super().parse_args(compute_options)

        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

        def get_required_int(name):
            val = compute_options.get_int_arg(name, None)

            if val is None:
                raise NexusProcessingException(reason=f'Missing required parameter: {name}', code=400)

            return val

        longitude = compute_options.get_argument('longitude', None)
        output = compute_options.get_content_type().lower()

        if longitude is None:
            raise NexusProcessingException(reason='Missing required parameter: longitude', code=400)
        if re.match(FLOAT_PATTERN, longitude):
            longitude = float(longitude)

            if longitude < -180 or longitude > 180:
                raise NexusProcessingException(reason='Invalid longitude: out of range', code=400)

            longitude = [longitude]
        elif ',' in longitude:
            if output == 'png':
                raise NexusProcessingException(reason='Cannot use PNG output with multi-longitude slicing', code=400)

            try:
                longitude = list(sorted([float(l) for l in longitude.split(',')]))
            except ValueError:
                raise NexusProcessingException(reason='Invalid longitude: invalid format', code=400)

            if any([l < -180 or l > 180 for l in longitude]):
                raise NexusProcessingException(reason='Invalid longitude: out of range', code=400)
        elif ':' in longitude:
            if output == 'png':
                raise NexusProcessingException(reason='Cannot use PNG output with multi-longitude slicing', code=400)

            try:
                longitude = slice(*(float(p) if p else None for p in longitude.split(':')))
            except:
                raise NexusProcessingException(reason='Invalid longitude: invalid format', code=400)

            if longitude.start is not None and longitude.stop is not None and longitude.start >= longitude.stop:
                raise NexusProcessingException(reason='Invalid longitude: out of order', code=400)

            # Some helpers for a < b / a > b where a could be None and return false in that case
            def lt(a, b):
                return a is not None and a < b

            def gt(a, b):
                return a is not None and a > b

            if (lt(longitude.start, -180) or gt(longitude.start, 180) or
                    lt(longitude.stop, -180) or gt(longitude.stop, 180)):
                raise NexusProcessingException(reason='Invalid longitude: out of range', code=400)
        elif longitude == '*':
            if output == 'png':
                raise NexusProcessingException(reason='Cannot use PNG output with multi-longitude slicing', code=400)

            longitude = slice(None, None, None)
        else:
            raise NexusProcessingException(reason='Invalid longitude argument', code=400)

        min_lat = get_required_float('minLat')
        max_lat = get_required_float('maxLat')
        min_elevation = get_required_int('minElevation')
        max_elevation = get_required_int('maxElevation')

        horizontal_margin = compute_options.get_float_arg('horizontalMargin', 0.001)

        ch_ds = compute_options.get_argument('canopy_ds', None)
        g_ds = compute_options.get_argument('ground_ds', None)

        if ch_ds is not None:
            ch_ds = [v.strip() for v in ch_ds.split(',')]
        else:
            ch_ds = []

        if len(ch_ds) > 2:
            raise NexusProcessingException(
                reason='Too many datasets provided for canopy_ds', code=400
            )

        if g_ds is not None:
            g_ds = [v.strip() for v in g_ds.split(',')]
        else:
            g_ds = []

        if len(g_ds) > 2:
            raise NexusProcessingException(
                reason='Too many datasets provided for ground_ds', code=400
            )

        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'
            )

        peaks = compute_options.get_boolean_arg('peaks') and not percentiles

        cmap = mpl.colormaps.get(
            compute_options.get_argument('cmap', 'viridis'),
            mpl.colormaps['viridis']
        )

        return (ds, parameter, longitude, min_lat, max_lat, min_elevation, max_elevation, horizontal_margin, peaks,
                ch_ds, g_ds, percentiles, cmap)

    def calc(self, compute_options, **args):
        (dataset, parameter, longitude, min_lat, max_lat, min_elevation, max_elevation, horizontal_margin, calc_peaks,
         ch_ds, g_ds, percentiles, cmap) = self.parse_args(compute_options)

        slices = dict(
            lat=slice(min_lat, max_lat),
            elevation=slice(min_elevation, max_elevation)
        )

        if isinstance(longitude, slice):
            slices['lon'] = longitude
        elif isinstance(longitude, list):
            slices['lon'] = slice(longitude[0] - horizontal_margin, longitude[-1] + horizontal_margin)
        else:
            slices['lon'] = slice(longitude - horizontal_margin, longitude + horizontal_margin)

        data_in_bounds = self.do_subset(dataset, parameter, slices, 0)
        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]

        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 isinstance(longitude, list):
            lon_slices = [ds.sel(longitude=l, method='nearest') for l in longitude]
            ds = xr.concat(lon_slices, dim='longitude').drop_duplicates('longitude')
        else:
            ds = ds.sel(longitude=longitude)

        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)

        if calc_peaks:
            peaks = []

            for i in range(len(ds.tomo.longitude)):
                lon_peaks = []
                ds_sub = ds.isel(longitude=i)

                for j in range(len(ds.tomo.latitude)):
                    col = ds_sub.isel(latitude=j)

                    lat = col.latitude.item()

                    try:
                        idx = col.tomo.argmax(dim='elevation').item()

                        lon_peaks.append((lat, col.isel(elevation=idx).elevation.item()))

                    except ValueError:
                        lon_peaks.append((lat, np.nan))

                peaks.append((np.array([p[0] for p in lon_peaks]),
                              np.array([p[1] for p in lon_peaks])))
        else:
            peaks = None

        return ProfileTomoResults(
            results=(ds, peaks),
            s={'longitude': longitude},
            coords=dict(x=ds.latitude, y=ds.elevation),
            meta=dict(dataset=dataset),
            style='db' if not percentiles else 'percentile',
            cmap=cmap,
            computeOptions=compute_options
        )


@nexus_handler
class LatitudeTomogramImpl(TomogramBaseClass):
    name = "Tomogram Latitude profile view"
    path = "/tomogram/latitude"
    description = "Fetches a tomogram sliced by latitude"
    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."
        },
        "latitude": {
            "name": "Latitude",
            "type": "float, comma-delimited floats(>1), colon-delimited floats(<=3), literal['*']",
            "description": "Latitude(s) at which to slice the tomogram. Required. Valid formats: Single float - Slice "
                           "at a singular latitude; comma-delimited floats - slice at a list of latitudes; "
                           "colon-delimited floats - slice using Python slicing notation (both ends are inclusive): "
                           "[<start>]:[<stop>]:[<step>] where omitting <start> and/or <stop> produce unbounded slices "
                           "and omitting <step> equates to an all-inclusive step size. A special value of \"*\" "
                           "equates to ::1. Values other than singular floats are incompatible with PNG rendering "
                           "(use ZIP to a ZIP archive of .pngs)"
        },
        "minLon": {
            "name": "Minimum Longitude",
            "type": "float",
            "description": "Minimum (Eastern) Longitude. Required."
        },
        "maxLon": {
            "name": "Maximum Longitude",
            "type": "float",
            "description": "Maximum (Western) Longitude. Required."
        },
        "minElevation": {
            "name": "Minimum Elevation",
            "type": "int",
            "description": "Minimum Elevation. Required."
        },
        "maxElevation": {
            "name": "Maximum Elevation",
            "type": "int",
            "description": "Maximum Elevation. Required."
        },
        "horizontalMargin": {
            "name": "Horizontal Margin",
            "type": "float",
            "description": "Margin +/- desired lat/lon slice to include in output. Default: 0.001m"
        },
        "peaks": {
            "name": "Calculate peaks",
            "type": "boolean",
            "description": "Calculate peaks along tomogram slice (currently uses simplest approach). "
                           "Ignored with \'elevPercentiles\'"
        },
        "canopy_ds": {
            "name": "Canopy height dataset name",
            "type": "string, comma-separated pair of strings",
            "description": "Dataset(s) containing canopy heights (RH98). This is used to trim out tomogram voxels that "
                           "return over the measured or predicted canopy. If 2 are provided, the first is considered "
                           "more accurate and will be favored for masking."
        },
        "ground_ds": {
            "name": "Ground height dataset name",
            "type": "string, comma-separated pair of strings",
            "description": "Dataset containing ground height (DEM). This is used to trim out tomogram voxels that "
                           "return below the measured or predicted ground. If 2 are provided, the first is considered "
                           "more accurate and will be favored for masking."
        },
        "elevPercentiles": {
            "name": "Elevation percentiles",
            "type": "boolean",
            "description": "Display percentiles of cumulative returns over elevation. Requires both canopy_ds and "
                           "ground_ds be set."
        },
        "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())]}"
        },
        "cbarMin": {
            "name": "Color bar min",
            "type": "float",
            "description": "Minimum value of the color bar. Only applies if elevPercentiles is False. Must be set with"
                           " and less than cbarMax. Invalid configuration will be replaced by default value (-30)"
        },
        "cbarMax": {
            "name": "Color bar max",
            "type": "float",
            "description": "Maximum value of the color bar. Only applies if elevPercentiles is False. Must be set with"
                           " and greater than cbarMin. Invalid configuration will be replaced by default value (-10)"
        },
    }
    singleton = True

    def __init__(self, tile_service_factory, **kwargs):
        TomogramBaseClass.__init__(self, tile_service_factory, **kwargs)

    def parse_args(self, compute_options):
        ds, parameter = super().parse_args(compute_options)

        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

        def get_required_int(name):
            val = compute_options.get_int_arg(name, None)

            if val is None:
                raise NexusProcessingException(reason=f'Missing required parameter: {name}', code=400)

            return val

        latitude = compute_options.get_argument('latitude', None)
        output = compute_options.get_content_type().lower()

        if latitude is None:
            raise NexusProcessingException(reason='Missing required parameter: latitude', code=400)
        if re.match(FLOAT_PATTERN, latitude):
            latitude = float(latitude)

            if latitude < -180 or latitude > 180:
                raise NexusProcessingException(reason='Invalid latitude: out of range', code=400)

            latitude = [latitude]
        elif ',' in latitude:
            if output == 'png':
                raise NexusProcessingException(reason='Cannot use PNG output with multi-latitude slicing', code=400)

            try:
                latitude = list(sorted([float(l) for l in latitude.split(',')]))
            except ValueError:
                raise NexusProcessingException(reason='Invalid latitude: invalid format', code=400)

            if any([l < -180 or l > 180] for l in latitude):
                raise NexusProcessingException(reason='Invalid latitude: out of range', code=400)
        elif ':' in latitude:
            if output == 'png':
                raise NexusProcessingException(reason='Cannot use PNG output with multi-longitude slicing', code=400)

            try:
                latitude = slice(*(float(p) if p else None for p in latitude.split(':')))
            except:
                raise NexusProcessingException(reason='Invalid latitude: invalid format', code=400)

            if latitude.start is not None and latitude.stop is not None and latitude.start >= latitude.stop:
                raise NexusProcessingException(reason='Invalid latitude: out of order', code=400)

            # Some helpers for a < b / a > b where a could be None and return false in that case
            def lt(a, b):
                return a is not None and a < b

            def gt(a, b):
                return a is not None and a > b

            if (lt(latitude.start, -180) or gt(latitude.start, 180) or
                    lt(latitude.stop, -180) or gt(latitude.stop, 180)):
                raise NexusProcessingException(reason='Invalid latitude: out of range', code=400)
        elif latitude == '*':
            if output == 'png':
                raise NexusProcessingException(reason='Cannot use PNG output with multi-latitude slicing', code=400)

            latitude = slice(None, None, None)
        else:
            raise NexusProcessingException(reason='Invalid latitude argument', code=400)

        min_lon = get_required_float('minLon')
        max_lon = get_required_float('maxLon')
        min_elevation = get_required_int('minElevation')
        max_elevation = get_required_int('maxElevation')

        horizontal_margin = compute_options.get_float_arg('horizontalMargin', 0.001)

        ch_ds = compute_options.get_argument('canopy_ds', None)
        g_ds = compute_options.get_argument('ground_ds', None)

        if ch_ds is not None:
            ch_ds = [v.strip() for v in ch_ds.split(',')]
        else:
            ch_ds = []

        if len(ch_ds) > 2:
            raise NexusProcessingException(
                reason='Too many datasets provided for canopy_ds', code=400
            )

        if g_ds is not None:
            g_ds = [v.strip() for v in g_ds.split(',')]
        else:
            g_ds = []

        if len(g_ds) > 2:
            raise NexusProcessingException(
                reason='Too many datasets provided for ground_ds', code=400
            )

        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'
            )

        peaks = compute_options.get_boolean_arg('peaks') and not percentiles

        cmap = mpl.colormaps.get(
            compute_options.get_argument('cmap', 'viridis'),
            mpl.colormaps['viridis']
        )

        return (ds, parameter, latitude, min_lon, max_lon, min_elevation, max_elevation,horizontal_margin, peaks, ch_ds,
                g_ds, percentiles, cmap)

    def calc(self, compute_options, **args):
        (dataset, parameter, latitude, min_lon, max_lon, min_elevation, max_elevation, horizontal_margin, calc_peaks,
         ch_ds, g_ds, percentiles, cmap) = self.parse_args(compute_options)

        slices = dict(
            lon=slice(min_lon, max_lon),
            elevation=slice(min_elevation, max_elevation)
        )

        if isinstance(latitude, slice):
            slices['lat'] = latitude
        elif isinstance(latitude, list):
            slices['lat'] = slice(latitude[0] - horizontal_margin, latitude[-1] + horizontal_margin)
        else:
            slices['lat'] = slice(latitude - horizontal_margin, latitude + horizontal_margin)

        data_in_bounds = self.do_subset(dataset, parameter, slices, 0)
        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]

        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 isinstance(latitude, list):
            lat_slices = [ds.sel(latitude=l, method='nearest') for l in latitude]
            ds = xr.concat(lat_slices, dim='latitude').drop_duplicates('latitude')
        else:
            ds = ds.sel(latitude=latitude)

        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)

        if calc_peaks:
            peaks = []

            for i in range(len(ds.tomo.latitude)):
                lat_peaks = []
                ds_sub = ds.isel(latitude=i)

                for j in range(len(ds.tomo.longitude)):
                    col = ds_sub.isel(longitude=j)

                    lon = col.longitude.item()

                    try:
                        idx = col.tomo.argmax(dim='elevation').item()

                        lat_peaks.append((lon, col.isel(elevation=idx).elevation.item()))

                    except ValueError:
                        lat_peaks.append((lon, np.nan))

                peaks.append((np.array([p[0] for p in lat_peaks]),
                              np.array([p[1] for p in lat_peaks])))
        else:
            peaks = None

        return ProfileTomoResults(
            results=(ds, peaks),
            s={'latitude': latitude},
            coords=dict(x=ds.longitude, y=ds.elevation),
            meta=dict(dataset=dataset),
            style='db' if not percentiles else 'percentile',
            cmap=cmap,
            computeOptions=compute_options
        )


@nexus_handler
class CustomProfileTomogramImpl(TomogramBaseClass):
    name = "Tomogram profile view between two points"
    path = "/tomogram/profile"
    description = "Fetches a tomogram sliced by a line between two points"
    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."
        },
        "p1": {
            "name": "Start Point",
            "type": "comma-delimited pair of floats",
            "description": "Point to start the slice line at: longitude,latitude, REQUIRED"
        },
        "p2": {
            "name": "End Point",
            "type": "comma-delimited pair of floats",
            "description": "Point to end the slice line at: longitude,latitude, REQUIRED"
        },
        "minElevation": {
            "name": "Minimum Elevation",
            "type": "int",
            "description": "Minimum Elevation. Required."
        },
        "maxElevation": {
            "name": "Maximum Elevation",
            "type": "int",
            "description": "Maximum Elevation. Required."
        },
        "horizontalMargin": {
            "name": "Horizontal Margin",
            "type": "float",
            "description": "Margin +/- desired lat/lon slice to include in output. Default: 0.001m"
        },
        "peaks": {
            "name": "Calculate peaks",
            "type": "boolean",
            "description": "Calculate peaks along tomogram slice (currently uses simplest approach). "
                           "Ignored with \'elevPercentiles\'"
        },
        "canopy_ds": {
            "name": "Canopy height dataset name",
            "type": "string, comma-separated pair of strings",
            "description": "Dataset(s) containing canopy heights (RH98). This is used to trim out tomogram voxels that "
                           "return over the measured or predicted canopy. If 2 are provided, the first is considered "
                           "more accurate and will be favored for masking."
        },
        "ground_ds": {
            "name": "Ground height dataset name",
            "type": "string, comma-separated pair of strings",
            "description": "Dataset containing ground height (DEM). This is used to trim out tomogram voxels that "
                           "return below the measured or predicted ground. If 2 are provided, the first is considered "
                           "more accurate and will be favored for masking."
        },
        "elevPercentiles": {
            "name": "Elevation percentiles",
            "type": "boolean",
            "description": "Display percentiles of cumulative returns over elevation. Requires both canopy_ds and "
                           "ground_ds be set."
        },
        "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())]}"
        },
        "cbarMin": {
            "name": "Color bar min",
            "type": "float",
            "description": "Minimum value of the color bar. Only applies if elevPercentiles is False. Must be set with"
                           " and less than cbarMax. Invalid configuration will be replaced by default value (-30)"
        },
        "cbarMax": {
            "name": "Color bar max",
            "type": "float",
            "description": "Maximum value of the color bar. Only applies if elevPercentiles is False. Must be set with"
                           " and greater than cbarMin. Invalid configuration will be replaced by default value (-10)"
        },
    }
    singleton = True

    def __init__(self, tile_service_factory, **kwargs):
        TomogramBaseClass.__init__(self, tile_service_factory, **kwargs)

    def parse_args(self, compute_options):
        ds, parameter = super().parse_args(compute_options)

        def get_required_int(name):
            val = compute_options.get_int_arg(name, None)

            if val is None:
                raise NexusProcessingException(reason=f'Missing required parameter: {name}', code=400)

            return val

        def get_required_point(name):
            val = compute_options.get_argument(name, None)

            if val is None:
                raise NexusProcessingException(reason=f'Missing required parameter: {name}', code=400)

            val = val.split(',')

            if len(val) != 2:
                raise NexusProcessingException(reason=f'Invalid parameter: {name}={val}. Incorrect number of fields, '
                                                      f'need 2', code=400)

            try:
                longitude = float(val[0])
            except ValueError:
                raise NexusProcessingException(reason='Longitude must be numeric', code=400)

            try:
                latitude = float(val[1])
            except ValueError:
                raise NexusProcessingException(reason='Latitude must be numeric', code=400)

            if longitude < -180 or longitude > 180 or isnan(longitude):
                raise NexusProcessingException(
                    reason='Longitude must be non-nan and between -180 and 180', code=400
                )

            if latitude < -90 or latitude > 90 or isnan(latitude):
                raise NexusProcessingException(
                    reason='Latitude must be non-nan and between -90 and 90', code=400
                )

            return Point(longitude, latitude)

        start_point = get_required_point('p1')
        end_point = get_required_point('p2')

        min_elevation = get_required_int('minElevation')
        max_elevation = get_required_int('maxElevation')

        if min_elevation >= max_elevation:
            raise NexusProcessingException(
                reason='Min elevation must be less than max elevation', code=400
            )

        horizontal_margin = compute_options.get_float_arg('horizontalMargin', 0.001)

        ch_ds = compute_options.get_argument('canopy_ds', None)
        g_ds = compute_options.get_argument('ground_ds', None)

        if ch_ds is not None:
            ch_ds = [v.strip() for v in ch_ds.split(',')]
        else:
            ch_ds = []

        if len(ch_ds) > 2:
            raise NexusProcessingException(
                reason='Too many datasets provided for canopy_ds', code=400
            )

        if g_ds is not None:
            g_ds = [v.strip() for v in g_ds.split(',')]
        else:
            g_ds = []

        if len(g_ds) > 2:
            raise NexusProcessingException(
                reason='Too many datasets provided for ground_ds', code=400
            )

        output = compute_options.get_content_type().lower()

        if output in ['csv', 'zip']:
            raise NexusProcessingException(reason=f'Output {output} unsupported', code=400)

        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'
            )

        peaks = compute_options.get_boolean_arg('peaks') and not percentiles

        cmap = mpl.colormaps.get(
            compute_options.get_argument('cmap', 'viridis'),
            mpl.colormaps['viridis']
        )

        return (ds, parameter, start_point, end_point, min_elevation, max_elevation,horizontal_margin, peaks, ch_ds,
                g_ds, percentiles, cmap)

    def __subset_on_line(self, ds, start_point, end_point, min_elevation, max_elevation, parameter, horizontal_margin):
        tile_service = self._get_tile_service()

        tiles = tile_service.find_tiles_along_line(
            start_point, end_point,
            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 tiles matched the selected paramters")

        data = []

        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

        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(bbox_miny, bbox_maxy, bbox_minx, bbox_maxx, [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
                })

        if len(data) == 0:
            raise NoDataException(reason='No data fit the selected parameters')

        return data

    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
        )


class ElevationTomoResults(NexusResults):
    def __init__(self, results=None, meta=None, stats=None, computeOptions=None, status_code=200, **args):
        NexusResults.__init__(self, results, meta, stats, computeOptions, status_code, **args)

    def toImage(self):
        data, lats, lons, attrs = self.results()

        lats = lats.tolist()
        lons = lons.tolist()

        logger.info('Building plot')

        plt.figure(figsize=(15,11))
        plt.imshow(np.squeeze(10*np.log10(data)), vmax=-10, vmin=-30)
        plt.colorbar(label=f'Tomogram, z={attrs["elevation"]} ± {attrs["margin"]} m (dB)')
        plt.title(f'{attrs["ds"]} tomogram @ {attrs["elevation"]} ± {attrs["margin"]} m elevation')
        plt.ylabel('Latitude')
        plt.xlabel('Longitude')

        xticks, xlabels = plt.xticks()
        xlabels = [f'{lons[int(t)]:.4f}' if int(t) in range(len(lons)) else '' for t in xticks]
        plt.xticks(xticks, xlabels, rotation=-90)

        yticks, ylabels = plt.yticks()
        ylabels = [f'{lats[int(t)]:.4f}' if int(t) in range(len(lats)) else '' for t in yticks]
        plt.yticks(yticks, ylabels, )

        buffer = BytesIO()

        logger.info('Writing plot to buffer')
        plt.savefig(buffer, format='png', facecolor='white')

        buffer.seek(0)
        plt.close()
        return buffer.read()


class ProfileTomoResults(NexusResults):
    def __init__(
            self,
            results=None,
            s=None,
            coords=None,
            meta=None,
            style='db',
            stats=None,
            computeOptions=None,
            status_code=200,
            **args
    ):
        NexusResults.__init__(self, results, meta, stats, None, status_code, **args)
        # Workaround to avoid error from non-standard parameters
        self.__computeOptions = computeOptions
        self.__coords = coords
        self.__slice = s
        self.__style = style
        self.__cmap = args.get('cmap', mpl.colormaps['viridis'])

    def toImage(self, i=0):
        ds, peaks = self.results()

        options = self.__computeOptions

        if 'longitude' in self.__slice:
            lon = self.__slice['longitude'][0]
            title_row = f'Longitude: {lon}'
            xlabel = 'Latitude'
            coord = 'longitude'
        elif 'latitude' in self.__slice:
            lat = self.__slice['latitude'][0]
            title_row = f'Latitude: {lat}'
            xlabel = 'Longitude'
            coord = 'latitude'
        else:
            start = self.__slice['start_point']
            end = self.__slice['end_point']

            start = (start.x, start.y)
            end = (end.x, end.y)

            title_row = f'Slice from: {start} to {end}'
            xlabel = 'Distance along slice [m]'
            coord = 'slice'

        ds = ds.isel({coord: i}).transpose(self.__coords['y'].name, self.__coords['x'].name)

        rows = ds.tomo.to_numpy()
        if peaks is not None:
            peaks = peaks[i]

        logger.info('Building plot')

        plt.figure(figsize=(11, 7))
        if self.__style == 'db':
            opt_min = options.get_float_arg('cbarMin', None)
            opt_max = options.get_float_arg('cbarMax', None)

            v = dict(vmin=-30, vmax=-10)

            if opt_min is not None:
                v['vmin'] = opt_min
            if opt_max is not None:
                v['vmax'] = opt_max

            if v['vmin'] >= v['vmax']:
                logger.warning('cbarMin >= cbarMax. Using defaults')
                v = dict(vmin=-30, vmax=-10)
        else:
            v = dict(vmin=0, vmax=1)
        plt.pcolormesh(self.__coords['x'], self.__coords['y'], rows, cmap=self.__cmap, **v)
        plt.colorbar()
        plt.title(f'{self.meta()["dataset"]} tomogram slice\n{title_row}')
        plt.xlabel(xlabel)
        plt.ylabel(f'Elevation w.r.t. {self.meta()["dataset"]} reference (m)')

        legend_handles = []

        if 'ch' in ds.data_vars and 'ch_secondary' in ds.data_vars:
            x_dim = ds['ch'].dims[0]

            # Plot secondary first so primary covers it
            second_ch = plt.plot(
                ds[x_dim].to_numpy(),
                ds.ch_secondary.to_numpy(),
                label=f'Secondary canopy height: {ds.ch_secondary.attrs["_source"]}'
            )

            prim_ch = plt.plot(
                ds[x_dim].to_numpy(),
                ds.ch.to_numpy(),
                label=f'Primary canopy height: {ds.ch.attrs["_source"]}'
            )

            legend_handles.extend(prim_ch)
            legend_handles.extend(second_ch)

        if 'gh' in ds.data_vars and 'gh_secondary' in ds.data_vars:
            x_dim = ds['gh'].dims[0]

            # Plot secondary first so primary covers it
            second_ch = plt.plot(
                ds[x_dim].to_numpy(),
                ds.gh_secondary.to_numpy(),
                label=f'Secondary ground height: {ds.gh_secondary.attrs["_source"]}'
            )

            prim_ch = plt.plot(
                ds[x_dim].to_numpy(),
                ds.gh.to_numpy(),
                label=f'Primary ground height: {ds.gh.attrs["_source"]}'
            )

            legend_handles.extend(prim_ch)
            legend_handles.extend(second_ch)

        if peaks is not None:
            peak_plot = plt.plot(
                peaks[0], peaks[1], color='red', alpha=0.75, label='Peak value', linestyle='--'
            )

            legend_handles.extend(peak_plot)

        plt.legend(handles=legend_handles)

        plt.ticklabel_format(useOffset=False)

        buffer = BytesIO()

        logger.info('Writing plot to buffer')
        plt.savefig(buffer, format='png', facecolor='white')

        buffer.seek(0)
        plt.close()
        return buffer.read()

    def toNetCDF(self):
        ds, peaks = self.results()

        if 'start_point' in self.__slice:
            coords = dict(slice=ds['slice'], distance=ds['distance'])
            dims = ['slice', 'distance']
        else:
            coords = dict(latitude=ds['latitude'], longitude=ds['longitude'])
            dims = ['longitude', 'latitude']

        if peaks is not None:
            ds['peaks'] = xr.DataArray(
                np.array([p[1] for p in peaks]),
                coords=coords,
                dims=dims,
                name='peaks'
            )

        comp = dict(zlib=True, complevel=9)
        encoding = {vname: comp for vname in ds.data_vars}

        with TemporaryDirectory() as td:
            ds.to_netcdf(os.path.join(td, 'ds.nc'), encoding=encoding)

            with open(os.path.join(td, 'ds.nc'), 'rb') as fp:
                buf = BytesIO(fp.read())

        buf.seek(0)
        return buf.read()

    def __plot_to_images(self, td):
        ds, peaks = self.results()

        options = self.__computeOptions

        if 'latitude' in self.__slice:
            along = 'latitude'
            across = 'longitude'
            title_row = 'Latitude: {:3.4f}'
            xlabel = 'Longitude'
        else:
            along = 'longitude'
            across = 'latitude'
            title_row = 'Longitude: {:3.4f}'
            xlabel = 'Latitude'

        min_across = ds[across].min().item()
        max_across = ds[across].max().item()

        files = []

        for al_i in range(len(ds[along])):
            ds_sub = ds.isel({along: al_i})
            al = ds_sub[along].item()

            rows = ds_sub.tomo.to_numpy()

            peaks_sub = None

            if peaks is not None:
                peaks_sub = peaks[al_i]

            logger.info(f'Building plot for along slice {al_i} of {len(ds[along])}')

            plt.figure(figsize=(11, 7))
            if self.__style == 'db':
                opt_min = options.get_float_arg('cbarMin', None)
                opt_max = options.get_float_arg('cbarMax', None)

                v = dict(vmin=-30, vmax=-10)

                if opt_min is not None:
                    v['vmin'] = opt_min
                if opt_max is not None:
                    v['vmax'] = opt_max

                if v['vmin'] >= v['vmax']:
                    logger.warning('cbarMin >= cbarMax. Using defaults')
                    v = dict(vmin=-30, vmax=-10)
            else:
                v = dict(vmin=0, vmax=1)

            plt.pcolormesh(self.__coords['x'], self.__coords['y'], rows, cmap=self.__cmap, **v)
            plt.colorbar()
            plt.title(f'{self.meta()["dataset"]} tomogram slice\n{title_row.format(al)}')
            plt.xlabel(xlabel)
            plt.ylabel(f'Elevation w.r.t. {self.meta()["dataset"]} reference (m)')

            legend_handles = []

            if 'ch' in ds_sub.data_vars and 'ch_secondary' in ds_sub.data_vars:
                x_dim = ds_sub['ch'].dims[0]

                # Plot secondary first so primary covers it
                second_ch = plt.plot(
                    ds_sub[x_dim].to_numpy(),
                    ds_sub.ch_secondary.to_numpy(),
                    label=f'Secondary canopy height: {ds_sub.ch_secondary.attrs["_source"]}'
                )

                prim_ch = plt.plot(
                    ds_sub[x_dim].to_numpy(),
                    ds_sub.ch.to_numpy(),
                    label=f'Primary canopy height: {ds_sub.ch.attrs["_source"]}'
                )

                legend_handles.extend(prim_ch)
                legend_handles.extend(second_ch)

            if 'gh' in ds_sub.data_vars and 'gh_secondary' in ds_sub.data_vars:
                x_dim = ds_sub['gh'].dims[0]

                # Plot secondary first so primary covers it
                second_ch = plt.plot(
                    ds_sub[x_dim].to_numpy(),
                    ds_sub.gh_secondary.to_numpy(),
                    label=f'Secondary ground height: {ds_sub.gh_secondary.attrs["_source"]}'
                )

                prim_ch = plt.plot(
                    ds_sub[x_dim].to_numpy(),
                    ds_sub.gh.to_numpy(),
                    label=f'Primary ground height: {ds_sub.gh.attrs["_source"]}'
                )

                legend_handles.extend(prim_ch)
                legend_handles.extend(second_ch)

            if peaks_sub is not None:
                peak_plot = plt.plot(
                    peaks_sub[0], peaks_sub[1], color='red', alpha=0.75, linestyle='--'
                )

                legend_handles.extend(peak_plot)

            plt.legend(handles=legend_handles, loc='lower right')

            plt.ticklabel_format(useOffset=False)

            filename = f'tomogram_{self.meta()["dataset"]}_{along}_{al:3.4f}_{across}_{min_across:3.4f}_to_{max_across:3.4f}.png'

            plt.savefig(
                os.path.join(
                    td, filename
                ), format='png', facecolor='white')

            files.append((filename, os.path.join(td, filename)))
            plt.close()

        return files

    def toZip(self):
        buffer = BytesIO()

        with TemporaryDirectory() as td:
            files = self.__plot_to_images(td)

            with zipfile.ZipFile(buffer, 'a', zipfile.ZIP_DEFLATED) as zf:
                for filename, file_path in files:
                    zf.write(file_path, filename)

        buffer.seek(0)
        return buffer.read()

    def toGif(self):
        buffer = BytesIO()

        with TemporaryDirectory() as td:
            files = self.__plot_to_images(td)

            with contextlib.ExitStack() as stack:
                logger.info(f'Combining frames into final GIF')

                imgs = (stack.enter_context(Image.open(f[1])) for f in files)
                img = next(imgs)
                img.save(
                    buffer,
                    format='GIF',
                    save_all=True,
                    append_images=imgs,
                    duration=([100] * (len(files) - 1)) + [800],
                    loop=0
                )

        buffer.seek(0)
        return buffer.read()