"""
Functions for preprocessing products from the MAP Framework Team
"""

import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from gisutils.raster import get_values_at_points, write_raster
from mfsetup import load_modelgrid
from mfsetup.discretization import voxels_to_layers, fill_cells_vertically
from mfsetup.testing import point_is_on_nhg
from mfsetup.units import convert_length_units


def get_layer(botm_array, i, j, elev):
    """Return the botm_array for elevations at i, j locations.

    Parameters
    ----------
    botm_array : 3D numpy array
        layer bottom elevations
    i : scaler or sequence
        row index (zero-based)
    j : scaler or sequence
        column index
    elev : scaler or sequence
        elevation (in same units as model)

    Returns
    -------
    k : np.ndarray (1-D) or scalar
        zero-based layer index
    """

    def to_array(arg):
        if not isinstance(arg, np.ndarray):
            return np.array([arg])
        else:
            return arg

    i = to_array(i)
    j = to_array(j)
    nlay = botm_array.shape[0]
    elev = to_array(elev)
    botms = botm_array[:, i, j]  # .tolist()
    # identify layer botms that are above and below the elevations
    differences = np.round((botms - elev), 2)
    isabove = differences >= 0
    # layer is the number of botm_array that are above
    layers = np.sum(isabove, axis=0)
    # force elevations below model bottom into bottom layer
    layers[layers > nlay - 1] = nlay - 1
    layers = np.atleast_1d(np.squeeze(layers))
    if len(layers) == 1:
        layers = layers[0]
    return layers


def plot_slice(layer_elevations, property_data=None,
               row=0, column=slice(None),
               voxel_start_layer=0, voxel_zones=None, cmap='copper',
               voxel_cmap='viridis', unit_labels=None, add_surfaces=None):
    """Plot a single cross section slice

    Parameters
    ----------
    layer_elevations : 3D numpy array
        Array of layer elevations, starting with the model top. 
        (Length equal to the number of botm_array + 1)
    property_data : 3D numpy array 
        Array of zone numbers generated by setup_model_layers.
    row : int or slice instance
        If a cross section along a row is desired, row should be a integer,
        and column should be a slice instance indicating the range of columns to include.
        by default, 0.
    column : int or slice instance
        If a cross section along a column is desired, column should be a integer,
        and row should be a slice instance indicating the range of rows to include.
        by default, slice(None), which includes all columns.
    voxel_start_layer : int, optional
        First layer with voxel data, by default 0
    voxel_zones : sequence, optional
        Zone numbers within property_data that are voxel-based, 
        by default None
    cmap : str, optional
        Matplotlib colormap for non-voxel zone numbers, by default 'copper',
        to contrast with colormap for voxel-based zone numbers.
    voxel_cmap : str, optional
        Matplotlib colormap for voxel-based zone numbers, by default 'viridis'.
    unit_labels : dict, optional
        Dictionary mapping non-voxel zone numbers to hydrogeologic units, 
        by default None
        
    Returns
    -------
    ax : matplotlib AxesSubplot instance for figure
    """
    # cross section code
    nlay, nrow, ncol = layer_elevations.shape
    # create meshgrid for rows or columns
    # along a row
    if isinstance(column, slice):
        # x, z = np.meshgrid(range(ncol), np.array(z_edges))
        # x = grid.xcellcenters[row, column]
        ncells = ncol
        title = 'Row {}'.format(row)
        xlabel = 'Column in model'
    # along a column
    else:
        # x, z = np.meshgrid(range(nrow), np.array(z_edges))
        # x = grid.ycellcenters[row, column]
        ncells = nrow
        title = 'Column {}'.format(column)
        xlabel = 'Row in model'
    x = np.arange(ncells)
    z = layer_elevations[:, row, column].copy()
    # since z is used to define cell edges in the pcolormesh (below)
    # z cannot be masked or have nan values
    # set missing data values (outside of the model footprint) in z
    # to -9999
    # pcolormesh will still skip these cells, as they are defined
    # as no data by the mask for the property array
    z_nodata = -9999
    z[np.isnan(z)] = z_nodata
    # zero values will result in pcolormesh edges that dip to zero
    # on the edge of nodata areas
    # fill these with previous value in either direction
    # first drop any indices along the edges
    for side in -1, 1:
        k, j = np.where(z == z_nodata)
        interior_zeros = (j > 0) & (j < z.shape[1] - 1)
        j = j[interior_zeros]
        k = k[interior_zeros]
        # then reassign the zero elevations
        z[k, j] = z[k, j+side]

    #z = np.ma.masked_where(np.isnan(z), z)
    thicknesses = np.diff(z, axis=0) * -1
    thicknesses[thicknesses <= 0] = 0.

    fig, ax = plt.subplots(figsize=(11, 8.5))

    # optionally plot a property such as resistivity facies
    if property_data is not None:
        # drop na values
        # (areas with no voxel data at any depth)
        #loc = ~np.all(z.mask, axis=0)
        data = property_data[:, row, column].copy()
        vmin, vmax = property_data.min(), property_data.max()

        #x = np.squeeze(x[loc])  # + [x[-1] + 1]
        #x = np.ma.masked_array(x, mask=~loc)
        zstart = voxel_start_layer
        zend = voxel_start_layer + property_data.shape[0] + 1
        z = np.squeeze(z[zstart:zend, :])

        #if not np.any(z):
        #    return
        #data = np.squeeze(data[:, loc])
        #thicknesses = np.squeeze(thicknesses[:, loc])
        if np.any(z) and voxel_zones is not None:
            # get the min max values for the existing framework
            # and voxel-based property data
            is_voxel_3D = np.isin(property_data, voxel_zones)
            vmin = property_data[~is_voxel_3D].min()
            vmax = property_data[~is_voxel_3D].max()
            voxel_vmin = np.min(voxel_zones)
            voxel_vmax = np.max(voxel_zones)

            is_voxel = np.isin(data, voxel_zones)
            voxel_mask = (thicknesses <= 0) | ~is_voxel
            data_mask = (thicknesses <= 0) | is_voxel
            voxel_data = np.ma.masked_array(data, mask=voxel_mask)
            data = np.ma.masked_array(data, mask=data_mask)
            pcm2 = ax.pcolormesh(x, z, voxel_data, alpha=0.5, cmap=voxel_cmap,
                                 vmin=voxel_vmin, vmax=voxel_vmax)
            # add separate axis for colorbar to control position
            fig.colorbar(pcm2, label='Resistivity facies class (in order of increasing resistivity)',
                         pad=0.1, shrink=0.3,
                         orientation='horizontal',
                         ticks=voxel_zones[::2].astype(int),
                         )
        pcm = ax.pcolormesh(x, z, data, alpha=0.5, cmap='copper',
                            vmin=vmin, vmax=vmax)

    if add_surfaces is not None:
        for label, array2D in add_surfaces.items():
            ax.plot(array2D[row, column], c='k', label=label)
        surfaces_lg = ax.legend(loc='lower right', bbox_to_anchor=(1, 0))
        ax.add_artist(surfaces_lg)

    # plot the layer bottom elevations
    plot_layer = np.any(thicknesses > 0, axis=1)
    plot_layer = np.append(plot_layer, True)
    plot_ij = np.any(thicknesses > 0, axis=0)
    z[z == z_nodata] = np.nan
    if np.any(plot_ij):
        for layer_edge in z[plot_layer, :]:
            ax.plot(x, layer_edge, color='k', lw=0.1)
        xmin = x[plot_ij].min()
        xmax = x[plot_ij].max()
        ymin = np.nanmin(z[plot_layer, :])
        ymax = np.nanmax(z[plot_layer, :])
        ax.set_xlim(xmin, xmax)
        ax.set_ylim(ymin, ymax)
        ax.set_xlabel(xlabel)
        ax.set_ylabel('Elevation, in meters')
        ax.set_title(title)

        if property_data is not None:
            # make a discrete colors legend
            legend_title = "MERAS 2.0\nFramework Units"
            lg = add_discrete_colors_legend(ax, data, unit_labels=unit_labels,
                                            cmap=cmap, vmin=vmin, vmax=vmax, alpha=0.5,
                                            bbox_to_anchor=(1, -0.1), loc='upper right',
                                            title=legend_title, handleheight=2)
        plt.tight_layout()
    return ax


def add_discrete_colors_legend(ax, data, unit_labels=None,
                               cmap='copper', vmin=None, vmax=None, alpha=1,
                               **kwargs):
    """Add a legend with a colored rectangle and label for each
    discrete value in data.

    Parameters
    ----------
    ax : matplotlib axes object
        An axes of the current figure.
    data : np.array or list-like
        Discrete data represented by legend.
    unit_labels : dict, optional
        Dictionary of unit labels keyed by integer value in data, 
        by default None, in which case all values in data
        are included with the integer value as the label.
    cmap : str, optional
        Matplotlib colormap, by default 'copper'
    vmin : float, optional
        Minimum value of colormap range, 
        by default None
    vmax : float, optional
        Maximum value of colormap range, , by default None
    **kwargs : keyword arguments to :func:`matplotlib.pyplot.legend`
    
    Return
    lg : :func:`matplotlib.pyplot.legend` handle
    
    """
    # make a discrete colors legend
    from matplotlib.colors import Normalize

    if isinstance(data, np.ma.masked_array):
        unique_zones = np.unique(data).compressed()
    else:
        unique_zones = np.unique(data)
    if not unit_labels:
        unit_labels = dict(list(zip(unique_zones, unique_zones)))
    else:
        unit_labels = {k: v for k, v in unit_labels.items() if k in unique_zones and len(v) > 0}
    cmap_obj = plt.get_cmap(cmap)

    # now get corresponding colormap value and label for each category and make legend
    norm = Normalize(vmin, vmax)  # normalize to range of 0-1
    handles = []
    labels_list = []
    for value, label in unit_labels.items():
        fc = cmap_obj(norm(value))
        handles.append(plt.Rectangle((0, 0), 1, 1, fc=fc, alpha=alpha))
        labels_list.append(label)
    if len(handles) > 0:
        lg = ax.legend(handles, labels_list, **kwargs)
        return lg


def plot_cross_sections(layers, out_pdf, property_data=None,
                        voxel_start_layer=0,  voxel_zones=None, cmap='copper',
                        voxel_cmap='viridis', unit_labels=None,
                        add_raster_surfaces=None,
                        modelgrid=None
                        ):
    """Generate a multi-page PDF of the layer cross sections.

    Parameters
    ----------
    layers : 3D numpy array 
        Array of layer elevations, starting with the model top. 
        (Length equal to the number of botm_array + 1)
    property_data : 3D numpy array 
        Array of zone numbers generated by setup_model_layers.
    out_pdf : str (filepath)
        Filename of multi-page PDF.
    voxel_start_layer : int, optional
        First layer with voxel data, by default 0
    voxel_zones : sequence, optional
        Zone numbers within property_data that are voxel-based, 
        by default None
    cmap : str, optional
        Matplotlib colormap for non-voxel zone numbers, by default 'copper',
        to contrast with colormap for voxel-based zone numbers.
    voxel_cmap : str, optional
        Matplotlib colormap for voxel-based zone numbers, by default 'viridis'.
    unit_labels : dict, optional
        Dictionary mapping non-voxel zone numbers to hydrogeologic units, 
        by default None
    """
    raster_arrays = None
    if add_raster_surfaces:
        if modelgrid is None:
            raise ValueError("add_raster_surfaces option requires a modelgrid")
        raster_arrays = {}
        _, nrow, ncol = modelgrid.shape
        x = modelgrid.xcellcenters.ravel()
        y = modelgrid.ycellcenters.ravel()
        for label, raster in add_raster_surfaces.items():
            values = get_values_at_points(raster, x, y, points_crs=modelgrid.crs)
            raster_arrays[label] = np.reshape(values, (nrow, ncol))

    with PdfPages(out_pdf) as pdf:
        nlay, nrow, ncol = layers.shape
        for row in range(0, nrow, 50):
            plot_slice(layers, property_data,
                       row=row, column=slice(None, None),
                       voxel_start_layer=voxel_start_layer,
                       voxel_zones=voxel_zones, cmap=cmap,
                       voxel_cmap=voxel_cmap, unit_labels=unit_labels,
                       add_surfaces=raster_arrays)
            pdf.savefig()
            plt.close()
        for column in range(0, ncol, 50):
            plot_slice(layers, property_data,
                       row=slice(None, None), column=column,
                       voxel_start_layer=voxel_start_layer,
                       voxel_zones=voxel_zones, cmap=cmap,
                       voxel_cmap=voxel_cmap, unit_labels=unit_labels,
                       add_surfaces=raster_arrays)
            pdf.savefig()
            plt.close()


def plot_zone_maps(layers, out_pdf, zones=None,
                   voxel_cmap='viridis', zones_cmap='copper',
                   voxel_zones=None, unit_labels=None):
    """Generate a multi-page PDF of the zones in each layer in plan view.

    Parameters
    ----------
    layers : 3D numpy array 
        Array of layer elevations, starting with the model top. 
        (Length equal to the number of botm_array + 1)
    zones : 3D numpy array 
        Array of zone numbers generated by setup_model_layers.
    out_pdf : str (filepath)
        Filename of multi-page PDF.
    voxel_zones : sequence, optional
        Zone numbers within property_data that are voxel-based, 
        by default None
    unit_labels : str, optional
        Matplotlib colormap for non-voxel zone numbers, by default 'copper',
        to contrast with colormap for voxel-based zone numbers.
    voxel_cmap : str, optional
        Matplotlib colormap for voxel-based zone numbers, by default 'viridis'.
    unit_labels : dict, optional
        Dictionary mapping non-voxel zone numbers to hydrogeologic units, 
        by default None
    """

    thickness = np.diff(layers, axis=0) * -1
    thickness[thickness <= 0] = 0.

    # split up arrays into voxel-based data (AEM)
    # and layer-based (existing framework)
    voxel_data = None
    if voxel_zones is not None:
        is_voxel = np.isin(zones, voxel_zones)
        voxel_data_mask = (thickness <= 0) | ~is_voxel

        # everything that isn't voxel is in the previous framework
        existing_framework_mask = (thickness <= 0) | is_voxel

        voxel_data = np.ma.masked_array(zones, mask=voxel_data_mask)
        voxel_vmin, voxel_vmax = voxel_data.min(), voxel_data.max()
        zones = np.ma.masked_array(zones, mask=existing_framework_mask)

    vmin, vmax = zones.min(), zones.max()

    if len(zones.shape) == 2:
        zones = [zones]

    with PdfPages(out_pdf) as pdf:
        for k, array in enumerate(zones):
            fig = plt.figure(figsize=(11, 8.5))
            ax = fig.add_axes([0.1, 0.1, 0.6, 0.8])
            ax.set_title('layer {}'.format(k))
            im = ax.imshow(array, cmap=zones_cmap, vmin=vmin, vmax=vmax)
            if voxel_zones is not None:
                im2 = ax.imshow(voxel_data[k], cmap=voxel_cmap, vmin=voxel_vmin, vmax=voxel_vmax)

            fig.colorbar(im2, label='Resistivity facies class (in order of increasing resistivity)',
                         pad=0.1, shrink=0.5, orientation='horizontal',
                         ticks=voxel_zones[::2].astype(int),
                         )

            # make a discrete colors legend
            title = "MERAS 2.0\nFramework Units"
            lg = add_discrete_colors_legend(ax, array, unit_labels=unit_labels,
                                            cmap=zones_cmap, vmin=vmin, vmax=vmax,
                                            bbox_to_anchor=(1, 0), loc='lower left',
                                            title=title, handleheight=1)
            plt.tight_layout()
            pdf.savefig()
            plt.close()


def rasters_to_grid(modelgrid, dem, rasters,
                    dem_elevation_units='meters', raster_elevation_units='meters',
                    dest_elevation_units='meters'):
    """Sample a sequence of rasters onto the i, j locations of a modelgrid,
    returning a 3D numpy array of the sampled elevations. Fill places with nodata
    using the next valid surface above.

    Parameters
    ----------
    modelgrid : Modflow-setup :class:`~mfsetup.grid.MFsetupGrid` instance
        Modflow-setup grid instance describing the model grid
    dem : str (filepath)
        Raster representing the land surface, at the highest resolution being contemplated for the model.
        Usually this is derived by sampling a higher resolution DEM using zonal statistics, taking
        the mean DEM value for each model cell.
    rasters : list of strings (filepaths)
        Raster surfaces describing hydrogelogic contacts surrounding the voxel data.
    dem_elevation_units : str, optional
        Elevation units of dem_means_raster, by default 'meters'
    framework_raster_elevation_units : str, optional
        Elevation units of the framework_rasters, by default 'meters'
    model_length_units : str, optional
        Length units used in the model, by default 'meters'

    References
    ----------
    See the documentation for the :func:`fill_cells_vertically <mfsetup.discretization.fill_cells_vertically>`
    function in Modflow-setup for an explanation of the filling process.

    """
    grid = modelgrid
    dem_elevations = get_values_at_points(dem, grid.xcellcenters, grid.ycellcenters,
                                          method='linear')
    # convert to model units
    dem_elevations *= convert_length_units(dem_elevation_units, dest_elevation_units)

    raster_elevations = []
    for raster in rasters:
        grid_cell_values = get_values_at_points(raster, grid.xcellcenters, grid.ycellcenters,
                                                method='linear')
        # convert to model units
        grid_cell_values *= convert_length_units(raster_elevation_units, dest_elevation_units)
        raster_elevations.append(grid_cell_values)
    raster_elevations = np.array(raster_elevations)
    
    # fill nans in the sampled original framework elevations
    # (nans are where a layer surface is absent)
    # fill the nans with the next surface above
    # see https://github.com/aleaf/modflow-setup/blob/develop/mfsetup/discretization.py
    model_top_filled, filled_raster_elevations = fill_cells_vertically(dem_elevations, raster_elevations)
    above_land_surface = filled_raster_elevations > dem_elevations
    # reset any values above land surface to land surface
    dem_means_3d = np.tile(dem_elevations, (filled_raster_elevations.shape[0], 1, 1))
    filled_raster_elevations[above_land_surface] = dem_means_3d[above_land_surface]
    del dem_means_3d
    filled_raster_elevations = np.vstack([np.reshape(dem_elevations, (1, *dem_elevations.shape)),
                                          filled_raster_elevations])
    return filled_raster_elevations
    

def layers_to_zones(botm_array, model_cell_z_centers):
    """Map the layers from one grid to zones on another grid. In other words,
    given a 3D array of model cell bottom elevations (botm_array), and a second
    3D array of cell centers at the same i, j locations as botm_array,
    but with different vertical discretization (model_cell_z_centers),
    assign each cell in the model_cell_z_centers array a zone value based on its
    vertical position within botm array. For example, a cell that lies above all of the
    elevations in botm_array would be assigned a zone number of zero; a cell that lies below
    the first bottom elevation would be assigned a value of 1.

    Parameters
    ----------
    botm_array : 3D numpy array
        Layer bottom elevations
    model_cell_z_centers : 3D numpy array
        Cell centers for another grid that is aligned with
        botm_array along the i and j axis, but has a different vertical
        discretation.

    Returns
    -------
    zones : 3D numpy array
        Zone numbers indicating the layer in botm_array for each
        cell in model_cell_z_centers.
    """

    # convert from 2D to 3D
    is2D = False
    if len(model_cell_z_centers.shape) == 2:
        is2D = True
        model_cell_z_centers = np.array([model_cell_z_centers])
    assert model_cell_z_centers.shape[1:] == botm_array.shape[1:]

    # for each, i, j location and cell midpoint in the model_cell_z_centers array
    # get the corresponding layer within the original meras framework
    # get the i, j locations for each grid cell
    i, j = np.indices(model_cell_z_centers[0].shape)
    i = i.ravel()
    j = j.ravel()
    zones = []
    for k, layer_centers in enumerate(model_cell_z_centers):
        layer_in_botm_array = get_layer(botm_array, i, j, layer_centers.ravel())
        layer_in_botm_array = np.reshape(layer_in_botm_array, layer_centers.shape)
        zones.append(layer_in_botm_array)
    zones = np.array(zones, dtype=int)
    if is2D:
        zones = zones[0]
    return zones


def setup_model_layers(dem_means_raster, facies_classes_netcdf, framework_rasters, modelgrid,
                       facies_class_variable='fac_a',
                       dem_elevation_units='meters', framework_raster_elevation_units='meters',
                       model_length_units='meters', output_folder='.',
                       framework_unit_names=None):
    """Generate model layering and property zones from voxel-based zone numbers at uniform depths
    and raster surfaces that represent hydrogeologic contacts. The voxel-based zone numbers may
    represent hydrogeologically different sediments, as determined from an airborne electromagnetic survey, 
    and are assumed to have priority over the raster surfaces. The model top is set at the land surface (DEM). 
    Subsequent cell bottoms beneath the land surface are set based on the voxel depths subtracted from the DEM elevations, 
    as sampled at the grid cell centers. Voxel cells without valid data are filled with zone numbers based
    on their position within the raster surfaces of hydrogeologic contacts. For example, a voxel
    cell that lies above all of the raster surfaces would be assigned a zone number equal to the
    highest voxel-based zone number + 3 (an arbitrary gap in the zone numbering between the voxel data
    and the raster surfaces). Similarly, beneath the lowest voxel cell at each location, grid cells are 
    assigned zone values based on their vertical position relative to the raster surfaces. Below the voxel data,
    botm_array and their bottom elevations correspond to the framework_raster surfaces. The lowermost
    (last) raster surface forms the model bottom.
    

    Parameters
    ----------
    dem_means_raster : str (filepath)
        Raster representing the land surface, at the highest resolution being contemplated for the model.
        Usually this is derived by sampling a higher resolution DEM using zonal statistics, taking
        the mean DEM value for each model cell.
    facies_classes_netcdf : str (filepath)
        NetCDF file with the voxel zone data, with a facies_class_variable that contains the zone numbers. 
    framework_rasters : list of strings (filepaths)
        Raster surfaces describing hydrogelogic contacts surrounding the voxel data.
    modelgrid : mfsetup.MFsetupGrid instance
        Modflow-setup grid instance describing the model grid
    facies_class_variable : str, optional
        Variable in facies_classes_netcdf containing the zone information, by default 'fac_a'
    dem_elevation_units : str, optional
        Elevation units of dem_means_raster, by default 'meters'
    framework_raster_elevation_units : str, optional
        Elevation units of the framework_rasters, by default 'meters'
    model_length_units : str, optional
        Length units used in the model, by default 'meters'
    output_folder : str, optional
        Location where results will be saved, by default '.'
    framework_unit_names : list, optional
        Unit names for the framework_rasters, by default None

    Returns
    -------
    botm_array : 3D numpy array
        Array of layer elevations, starting with the model top. 
        (Length equal to the number of botm_array + 1)
    zone_array : 3D numpy array 
        Array of zone numbers, including the values from facies_classes_netcdf,
        followed by a gap of 3 and then the numbers for the framework_rasters.
        For example, if there are 5 zones in facies_classes_netcdf, and 3 raster surfaces
        that intersect cells outside of the voxel data, zones 1-5 would correspond
        to the voxel zones, and zones 8-10 would correspond to hydrogeologic units
        bounded by the framework_rasters.
        
    Notes
    -----
    In addition to the variables returned, the results are written
    to text-array and GeoTiff files within the output_folder.
    
    Two PDFs of figures are also generated- one showing the layering
    in cross sections at regular intervals throughout the model grid, and one
    showing the zones within each layer in plan view.
    """    

    # output folders
    figures_folder = os.path.join(output_folder, 'figures')
    layers_folder = os.path.join(output_folder, 'botm_array')
    zones_folder = os.path.join(output_folder, 'zones')
    # zone rasters are saved to zones_folder/rasters
    output_folders = [figures_folder,
                      layers_folder,
                      zones_folder,
                      os.path.join(zones_folder, 'rasters'),

                      ]
    for folder in output_folders:
        if not os.path.isdir(folder):
            os.makedirs(folder)

    # output PDF with cross sections
    out_xsection_pdf = os.path.join(figures_folder, 'zone_xsections.pdf')
    out_maps_pdf = os.path.join(figures_folder, 'zone_maps.pdf')

    # flopy model grid
    grid = modelgrid

    # read in the facies classes and resisitivity values
    # from framework team
    # and verify that the grid is aligned with the nhg
    ds = xr.load_dataset(facies_classes_netcdf)
    assert point_is_on_nhg(ds.x, ds.y, offset='center')

    # read in the dem and framework rasters and fill locations of no data
    # with the next valid elevation above
    filled_framework_layers = rasters_to_grid(modelgrid,
                                              dem_means_raster, framework_rasters,
                                              dem_elevation_units=dem_elevation_units,
                                              raster_elevation_units=framework_raster_elevation_units,
                                              dest_elevation_units=model_length_units)
    dem_means = filled_framework_layers[0]
    framework_layer_botms_filled = filled_framework_layers[1:]

    # make AEM framework z edge elevations by subtracting off depths from model top
    voxel_z_edges = []
    for depth in ds.zb.values:
        voxel_z_edges.append(dem_means - depth)
    # get the elevations of the voxel centers as well
    voxel_z_centers = []
    for depth in ds.z.values:
        voxel_z_centers.append(dem_means - depth)
    voxel_z_centers = np.array(voxel_z_centers)

    # get the resistivity facies for each model cell (nearest neighbor)
    voxel_array = ds[facies_class_variable].sel(x=grid.xcellcenters[0],
                                                y=grid.ycellcenters[:, 0], method='nearest')
    voxel_zones = np.unique(voxel_array.values[~np.isnan(voxel_array.values)])

    # fill areas of voxel_array that don't have AEM data with zone values
    # representing units in the original MERAS framework
    framework_layer = layers_to_zones(framework_layer_botms_filled, voxel_z_centers)

    # start framework zone numbers after resistivity facies, +3
    framework_zones = framework_layer + np.nanmax(voxel_array.values) + 3
    voxel_array = voxel_array.values
    voxel_array[np.isnan(voxel_array)] = framework_zones[np.isnan(voxel_array)]

    # put the botm_array together
    layers = voxels_to_layers(voxel_array, z_edges=voxel_z_edges, model_top=dem_means,
                              model_botm=framework_layer_botms_filled, no_data_value=0)

    # add framework zones below the voxel botm_array
    nlay, nrow, ncol = layers.shape
    nlay -= 1  # botm_array includes the model top
    framework_layers = list(range(voxel_array.shape[0], nlay))
    first_framework_zone = framework_zones.min()
    n_framwork_zones = framework_layer_botms_filled.shape[0]
    framework_zone_numbers = np.arange(first_framework_zone, first_framework_zone + n_framwork_zones)

    # assume that the last len(framework_layers) are sequential
    framework_zones_below = []
    for zone_number in framework_zone_numbers[-len(framework_layers):]:
        framework_zones_below.append(np.ones((nrow, ncol)) * zone_number)
    framework_zones_below = np.array(framework_zones_below)

    # make an array of all zones, including the resisitivity facies
    # and the zones representing the units in the original meras framework
    zone_array = np.vstack([voxel_array, framework_zones_below])

    # plot cross sections
    framework_unit_labels = dict(zip(framework_zone_numbers.astype(int), framework_unit_names))
    plot_cross_sections(layers, out_xsection_pdf, property_data=zone_array,
                        voxel_start_layer=0, voxel_zones=voxel_zones,
                        cmap='copper',
                        voxel_cmap='viridis', unit_labels=framework_unit_labels)

    # plot maps showing distribution of zones in each layer
    plot_zone_maps(layers, out_maps_pdf, zones=zone_array,
                   voxel_zones=voxel_zones, unit_labels=framework_unit_labels)

    # write out the layer elevation rasters
    write_raster(os.path.join(layers_folder, 'model_top.tif'), layers[0], xul=grid.xul, yul=grid.yul,
                 dx=grid.delc[0], dy=-grid.delr[0], rotation=0., proj_str='epsg:5070',
                 nodata=-9999, fieldname='elevation')
    for i, layer in enumerate(layers[1:]):
        write_raster(os.path.join(layers_folder, 'botm{}.tif'.format(i)),
                     layer, xul=grid.xul, yul=grid.yul,
                     dx=grid.delc[0], dy=-grid.delr[0],
                     rotation=0., crs=grid.crs,
                     nodata=-9999, fieldname='elevation')

    # write out zone arrays (in GeoTiff and text format)
    voxel_start_layer = 0
    voxel_end_layer = voxel_start_layer + voxel_array.shape[0]
    for i, zones2d in enumerate(zone_array):
        np.savetxt(os.path.join(zones_folder, 'res_fac{}.dat'.format(i)),
                   zones2d, fmt='%.0f')
        write_raster(os.path.join(zones_folder, 'rasters/res_fac{}.tif'.format(i)),
                     zones2d, xul=grid.xul, yul=grid.yul,
                     dx=grid.delc[0], dy=-grid.delr[0],
                     rotation=0., crs=grid.crs,
                     nodata=-9999, fieldname='elevation')
    return layers, zone_array
"""
Functions for preprocessing products from the MAP Framework Team
"""

import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from gisutils.raster import get_values_at_points, write_raster
from mfsetup import load_modelgrid
from mfsetup.discretization import voxels_to_layers, fill_cells_vertically
from mfsetup.testing import point_is_on_nhg
from mfsetup.units import convert_length_units


def get_layer(botm_array, i, j, elev):
    """Return the botm_array for elevations at i, j locations.

    Parameters
    ----------
    botm_array : 3D numpy array
        layer bottom elevations
    i : scaler or sequence
        row index (zero-based)
    j : scaler or sequence
        column index
    elev : scaler or sequence
        elevation (in same units as model)

    Returns
    -------
    k : np.ndarray (1-D) or scalar
        zero-based layer index
    """

    def to_array(arg):
        if not isinstance(arg, np.ndarray):
            return np.array([arg])
        else:
            return arg

    i = to_array(i)
    j = to_array(j)
    nlay = botm_array.shape[0]
    elev = to_array(elev)
    botms = botm_array[:, i, j]  # .tolist()
    # identify layer botms that are above and below the elevations
    differences = np.round((botms - elev), 2)
    isabove = differences >= 0
    # layer is the number of botm_array that are above
    layers = np.sum(isabove, axis=0)
    # force elevations below model bottom into bottom layer
    layers[layers > nlay - 1] = nlay - 1
    layers = np.atleast_1d(np.squeeze(layers))
    if len(layers) == 1:
        layers = layers[0]
    return layers


def plot_slice(layer_elevations, property_data=None,
               row=0, column=slice(None),
               voxel_start_layer=0, voxel_zones=None, cmap='copper',
               voxel_cmap='viridis', unit_labels=None, add_surfaces=None):
    """Plot a single cross section slice

    Parameters
    ----------
    layer_elevations : 3D numpy array
        Array of layer elevations, starting with the model top. 
        (Length equal to the number of botm_array + 1)
    property_data : 3D numpy array 
        Array of zone numbers generated by setup_model_layers.
    row : int or slice instance
        If a cross section along a row is desired, row should be a integer,
        and column should be a slice instance indicating the range of columns to include.
        by default, 0.
    column : int or slice instance
        If a cross section along a column is desired, column should be a integer,
        and row should be a slice instance indicating the range of rows to include.
        by default, slice(None), which includes all columns.
    voxel_start_layer : int, optional
        First layer with voxel data, by default 0
    voxel_zones : sequence, optional
        Zone numbers within property_data that are voxel-based, 
        by default None
    cmap : str, optional
        Matplotlib colormap for non-voxel zone numbers, by default 'copper',
        to contrast with colormap for voxel-based zone numbers.
    voxel_cmap : str, optional
        Matplotlib colormap for voxel-based zone numbers, by default 'viridis'.
    unit_labels : dict, optional
        Dictionary mapping non-voxel zone numbers to hydrogeologic units, 
        by default None
        
    Returns
    -------
    ax : matplotlib AxesSubplot instance for figure
    """
    # cross section code
    nlay, nrow, ncol = layer_elevations.shape
    # create meshgrid for rows or columns
    # along a row
    if isinstance(column, slice):
        # x, z = np.meshgrid(range(ncol), np.array(z_edges))
        # x = grid.xcellcenters[row, column]
        ncells = ncol
        title = 'Row {}'.format(row)
        xlabel = 'Column in model'
    # along a column
    else:
        # x, z = np.meshgrid(range(nrow), np.array(z_edges))
        # x = grid.ycellcenters[row, column]
        ncells = nrow
        title = 'Column {}'.format(column)
        xlabel = 'Row in model'
    x = np.arange(ncells)
    z = layer_elevations[:, row, column].copy()
    # since z is used to define cell edges in the pcolormesh (below)
    # z cannot be masked or have nan values
    # set missing data values (outside of the model footprint) in z
    # to -9999
    # pcolormesh will still skip these cells, as they are defined
    # as no data by the mask for the property array
    z_nodata = -9999
    z[np.isnan(z)] = z_nodata
    # zero values will result in pcolormesh edges that dip to zero
    # on the edge of nodata areas
    # fill these with previous value in either direction
    # first drop any indices along the edges
    for side in -1, 1:
        k, j = np.where(z == z_nodata)
        interior_zeros = (j > 0) & (j < z.shape[1] - 1)
        j = j[interior_zeros]
        k = k[interior_zeros]
        # then reassign the zero elevations
        z[k, j] = z[k, j+side]

    #z = np.ma.masked_where(np.isnan(z), z)
    thicknesses = np.diff(z, axis=0) * -1
    thicknesses[thicknesses <= 0] = 0.

    fig, ax = plt.subplots(figsize=(11, 8.5))

    # optionally plot a property such as resistivity facies
    if property_data is not None:
        # drop na values
        # (areas with no voxel data at any depth)
        #loc = ~np.all(z.mask, axis=0)
        data = property_data[:, row, column].copy()
        vmin, vmax = property_data.min(), property_data.max()

        #x = np.squeeze(x[loc])  # + [x[-1] + 1]
        #x = np.ma.masked_array(x, mask=~loc)
        zstart = voxel_start_layer
        zend = voxel_start_layer + property_data.shape[0] + 1
        z = np.squeeze(z[zstart:zend, :])

        #if not np.any(z):
        #    return
        #data = np.squeeze(data[:, loc])
        #thicknesses = np.squeeze(thicknesses[:, loc])
        if np.any(z) and voxel_zones is not None:
            # get the min max values for the existing framework
            # and voxel-based property data
            is_voxel_3D = np.isin(property_data, voxel_zones)
            vmin = property_data[~is_voxel_3D].min()
            vmax = property_data[~is_voxel_3D].max()
            voxel_vmin = np.min(voxel_zones)
            voxel_vmax = np.max(voxel_zones)

            is_voxel = np.isin(data, voxel_zones)
            voxel_mask = (thicknesses <= 0) | ~is_voxel
            data_mask = (thicknesses <= 0) | is_voxel
            voxel_data = np.ma.masked_array(data, mask=voxel_mask)
            data = np.ma.masked_array(data, mask=data_mask)
            pcm2 = ax.pcolormesh(x, z, voxel_data, alpha=0.5, cmap=voxel_cmap,
                                 vmin=voxel_vmin, vmax=voxel_vmax)
            # add separate axis for colorbar to control position
            fig.colorbar(pcm2, label='Resistivity facies class (in order of increasing resistivity)',
                         pad=0.1, shrink=0.3,
                         orientation='horizontal',
                         ticks=voxel_zones[::2].astype(int),
                         )
        pcm = ax.pcolormesh(x, z, data, alpha=0.5, cmap='copper',
                            vmin=vmin, vmax=vmax)

    if add_surfaces is not None:
        for label, array2D in add_surfaces.items():
            ax.plot(array2D[row, column], c='k', label=label)
        surfaces_lg = ax.legend(loc='lower right', bbox_to_anchor=(1, 0))
        ax.add_artist(surfaces_lg)

    # plot the layer bottom elevations
    plot_layer = np.any(thicknesses > 0, axis=1)
    plot_layer = np.append(plot_layer, True)
    plot_ij = np.any(thicknesses > 0, axis=0)
    z[z == z_nodata] = np.nan
    if np.any(plot_ij):
        for layer_edge in z[plot_layer, :]:
            ax.plot(x, layer_edge, color='k', lw=0.1)
        xmin = x[plot_ij].min()
        xmax = x[plot_ij].max()
        ymin = np.nanmin(z[plot_layer, :])
        ymax = np.nanmax(z[plot_layer, :])
        ax.set_xlim(xmin, xmax)
        ax.set_ylim(ymin, ymax)
        ax.set_xlabel(xlabel)
        ax.set_ylabel('Elevation, in meters')
        ax.set_title(title)

        if property_data is not None:
            # make a discrete colors legend
            legend_title = "MERAS 2.0\nFramework Units"
            lg = add_discrete_colors_legend(ax, data, unit_labels=unit_labels,
                                            cmap=cmap, vmin=vmin, vmax=vmax, alpha=0.5,
                                            bbox_to_anchor=(1, -0.1), loc='upper right',
                                            title=legend_title, handleheight=2)
        plt.tight_layout()
    return ax


def add_discrete_colors_legend(ax, data, unit_labels=None,
                               cmap='copper', vmin=None, vmax=None, alpha=1,
                               **kwargs):
    """Add a legend with a colored rectangle and label for each
    discrete value in data.

    Parameters
    ----------
    ax : matplotlib axes object
        An axes of the current figure.
    data : np.array or list-like
        Discrete data represented by legend.
    unit_labels : dict, optional
        Dictionary of unit labels keyed by integer value in data, 
        by default None, in which case all values in data
        are included with the integer value as the label.
    cmap : str, optional
        Matplotlib colormap, by default 'copper'
    vmin : float, optional
        Minimum value of colormap range, 
        by default None
    vmax : float, optional
        Maximum value of colormap range, , by default None
    **kwargs : keyword arguments to :func:`matplotlib.pyplot.legend`
    
    Return
    lg : :func:`matplotlib.pyplot.legend` handle
    
    """
    # make a discrete colors legend
    from matplotlib.colors import Normalize

    if isinstance(data, np.ma.masked_array):
        unique_zones = np.unique(data).compressed()
    else:
        unique_zones = np.unique(data)
    if not unit_labels:
        unit_labels = dict(list(zip(unique_zones, unique_zones)))
    else:
        unit_labels = {k: v for k, v in unit_labels.items() if k in unique_zones and len(v) > 0}
    cmap_obj = plt.get_cmap(cmap)

    # now get corresponding colormap value and label for each category and make legend
    norm = Normalize(vmin, vmax)  # normalize to range of 0-1
    handles = []
    labels_list = []
    for value, label in unit_labels.items():
        fc = cmap_obj(norm(value))
        handles.append(plt.Rectangle((0, 0), 1, 1, fc=fc, alpha=alpha))
        labels_list.append(label)
    if len(handles) > 0:
        lg = ax.legend(handles, labels_list, **kwargs)
        return lg


def plot_cross_sections(layers, out_pdf, property_data=None,
                        voxel_start_layer=0,  voxel_zones=None, cmap='copper',
                        voxel_cmap='viridis', unit_labels=None,
                        add_raster_surfaces=None,
                        modelgrid=None
                        ):
    """Generate a multi-page PDF of the layer cross sections.

    Parameters
    ----------
    layers : 3D numpy array 
        Array of layer elevations, starting with the model top. 
        (Length equal to the number of botm_array + 1)
    property_data : 3D numpy array 
        Array of zone numbers generated by setup_model_layers.
    out_pdf : str (filepath)
        Filename of multi-page PDF.
    voxel_start_layer : int, optional
        First layer with voxel data, by default 0
    voxel_zones : sequence, optional
        Zone numbers within property_data that are voxel-based, 
        by default None
    cmap : str, optional
        Matplotlib colormap for non-voxel zone numbers, by default 'copper',
        to contrast with colormap for voxel-based zone numbers.
    voxel_cmap : str, optional
        Matplotlib colormap for voxel-based zone numbers, by default 'viridis'.
    unit_labels : dict, optional
        Dictionary mapping non-voxel zone numbers to hydrogeologic units, 
        by default None
    """
    raster_arrays = None
    if add_raster_surfaces:
        if modelgrid is None:
            raise ValueError("add_raster_surfaces option requires a modelgrid")
        raster_arrays = {}
        _, nrow, ncol = modelgrid.shape
        x = modelgrid.xcellcenters.ravel()
        y = modelgrid.ycellcenters.ravel()
        for label, raster in add_raster_surfaces.items():
            values = get_values_at_points(raster, x, y, points_crs=modelgrid.crs)
            raster_arrays[label] = np.reshape(values, (nrow, ncol))

    with PdfPages(out_pdf) as pdf:
        nlay, nrow, ncol = layers.shape
        for row in range(0, nrow, 50):
            plot_slice(layers, property_data,
                       row=row, column=slice(None, None),
                       voxel_start_layer=voxel_start_layer,
                       voxel_zones=voxel_zones, cmap=cmap,
                       voxel_cmap=voxel_cmap, unit_labels=unit_labels,
                       add_surfaces=raster_arrays)
            pdf.savefig()
            plt.close()
        for column in range(0, ncol, 50):
            plot_slice(layers, property_data,
                       row=slice(None, None), column=column,
                       voxel_start_layer=voxel_start_layer,
                       voxel_zones=voxel_zones, cmap=cmap,
                       voxel_cmap=voxel_cmap, unit_labels=unit_labels,
                       add_surfaces=raster_arrays)
            pdf.savefig()
            plt.close()


def plot_zone_maps(layers, out_pdf, zones=None,
                   voxel_cmap='viridis', zones_cmap='copper',
                   voxel_zones=None, unit_labels=None):
    """Generate a multi-page PDF of the zones in each layer in plan view.

    Parameters
    ----------
    layers : 3D numpy array 
        Array of layer elevations, starting with the model top. 
        (Length equal to the number of botm_array + 1)
    zones : 3D numpy array 
        Array of zone numbers generated by setup_model_layers.
    out_pdf : str (filepath)
        Filename of multi-page PDF.
    voxel_zones : sequence, optional
        Zone numbers within property_data that are voxel-based, 
        by default None
    unit_labels : str, optional
        Matplotlib colormap for non-voxel zone numbers, by default 'copper',
        to contrast with colormap for voxel-based zone numbers.
    voxel_cmap : str, optional
        Matplotlib colormap for voxel-based zone numbers, by default 'viridis'.
    unit_labels : dict, optional
        Dictionary mapping non-voxel zone numbers to hydrogeologic units, 
        by default None
    """

    thickness = np.diff(layers, axis=0) * -1
    thickness[thickness <= 0] = 0.

    # split up arrays into voxel-based data (AEM)
    # and layer-based (existing framework)
    voxel_data = None
    if voxel_zones is not None:
        is_voxel = np.isin(zones, voxel_zones)
        voxel_data_mask = (thickness <= 0) | ~is_voxel

        # everything that isn't voxel is in the previous framework
        existing_framework_mask = (thickness <= 0) | is_voxel

        voxel_data = np.ma.masked_array(zones, mask=voxel_data_mask)
        voxel_vmin, voxel_vmax = voxel_data.min(), voxel_data.max()
        zones = np.ma.masked_array(zones, mask=existing_framework_mask)

    vmin, vmax = zones.min(), zones.max()

    if len(zones.shape) == 2:
        zones = [zones]

    with PdfPages(out_pdf) as pdf:
        for k, array in enumerate(zones):
            fig = plt.figure(figsize=(11, 8.5))
            ax = fig.add_axes([0.1, 0.1, 0.6, 0.8])
            ax.set_title('layer {}'.format(k))
            im = ax.imshow(array, cmap=zones_cmap, vmin=vmin, vmax=vmax)
            if voxel_zones is not None:
                im2 = ax.imshow(voxel_data[k], cmap=voxel_cmap, vmin=voxel_vmin, vmax=voxel_vmax)

            fig.colorbar(im2, label='Resistivity facies class (in order of increasing resistivity)',
                         pad=0.1, shrink=0.5, orientation='horizontal',
                         ticks=voxel_zones[::2].astype(int),
                         )

            # make a discrete colors legend
            title = "MERAS 2.0\nFramework Units"
            lg = add_discrete_colors_legend(ax, array, unit_labels=unit_labels,
                                            cmap=zones_cmap, vmin=vmin, vmax=vmax,
                                            bbox_to_anchor=(1, 0), loc='lower left',
                                            title=title, handleheight=1)
            plt.tight_layout()
            pdf.savefig()
            plt.close()


def rasters_to_grid(modelgrid, dem, rasters,
                    dem_elevation_units='meters', raster_elevation_units='meters',
                    dest_elevation_units='meters'):
    """Sample a sequence of rasters onto the i, j locations of a modelgrid,
    returning a 3D numpy array of the sampled elevations. Fill places with nodata
    using the next valid surface above.

    Parameters
    ----------
    modelgrid : Modflow-setup :class:`~mfsetup.grid.MFsetupGrid` instance
        Modflow-setup grid instance describing the model grid
    dem : str (filepath)
        Raster representing the land surface, at the highest resolution being contemplated for the model.
        Usually this is derived by sampling a higher resolution DEM using zonal statistics, taking
        the mean DEM value for each model cell.
    rasters : list of strings (filepaths)
        Raster surfaces describing hydrogelogic contacts surrounding the voxel data.
    dem_elevation_units : str, optional
        Elevation units of dem_means_raster, by default 'meters'
    framework_raster_elevation_units : str, optional
        Elevation units of the framework_rasters, by default 'meters'
    model_length_units : str, optional
        Length units used in the model, by default 'meters'

    References
    ----------
    See the documentation for the :func:`fill_cells_vertically <mfsetup.discretization.fill_cells_vertically>`
    function in Modflow-setup for an explanation of the filling process.

    """
    grid = modelgrid
    dem_elevations = get_values_at_points(dem, grid.xcellcenters, grid.ycellcenters,
                                          method='linear')
    # convert to model units
    dem_elevations *= convert_length_units(dem_elevation_units, dest_elevation_units)

    raster_elevations = []
    for raster in rasters:
        grid_cell_values = get_values_at_points(raster, grid.xcellcenters, grid.ycellcenters,
                                                method='linear')
        # convert to model units
        grid_cell_values *= convert_length_units(raster_elevation_units, dest_elevation_units)
        raster_elevations.append(grid_cell_values)
    raster_elevations = np.array(raster_elevations)
    
    # fill nans in the sampled original framework elevations
    # (nans are where a layer surface is absent)
    # fill the nans with the next surface above
    # see https://github.com/aleaf/modflow-setup/blob/develop/mfsetup/discretization.py
    model_top_filled, filled_raster_elevations = fill_cells_vertically(dem_elevations, raster_elevations)
    above_land_surface = filled_raster_elevations > dem_elevations
    # reset any values above land surface to land surface
    dem_means_3d = np.tile(dem_elevations, (filled_raster_elevations.shape[0], 1, 1))
    filled_raster_elevations[above_land_surface] = dem_means_3d[above_land_surface]
    del dem_means_3d
    filled_raster_elevations = np.vstack([np.reshape(dem_elevations, (1, *dem_elevations.shape)),
                                          filled_raster_elevations])
    return filled_raster_elevations
    

def layers_to_zones(botm_array, model_cell_z_centers):
    """Map the layers from one grid to zones on another grid. In other words,
    given a 3D array of model cell bottom elevations (botm_array), and a second
    3D array of cell centers at the same i, j locations as botm_array,
    but with different vertical discretization (model_cell_z_centers),
    assign each cell in the model_cell_z_centers array a zone value based on its
    vertical position within botm array. For example, a cell that lies above all of the
    elevations in botm_array would be assigned a zone number of zero; a cell that lies below
    the first bottom elevation would be assigned a value of 1.

    Parameters
    ----------
    botm_array : 3D numpy array
        Layer bottom elevations
    model_cell_z_centers : 3D numpy array
        Cell centers for another grid that is aligned with
        botm_array along the i and j axis, but has a different vertical
        discretation.

    Returns
    -------
    zones : 3D numpy array
        Zone numbers indicating the layer in botm_array for each
        cell in model_cell_z_centers.
    """

    # convert from 2D to 3D
    is2D = False
    if len(model_cell_z_centers.shape) == 2:
        is2D = True
        model_cell_z_centers = np.array([model_cell_z_centers])
    assert model_cell_z_centers.shape[1:] == botm_array.shape[1:]

    # for each, i, j location and cell midpoint in the model_cell_z_centers array
    # get the corresponding layer within the original meras framework
    # get the i, j locations for each grid cell
    i, j = np.indices(model_cell_z_centers[0].shape)
    i = i.ravel()
    j = j.ravel()
    zones = []
    for k, layer_centers in enumerate(model_cell_z_centers):
        layer_in_botm_array = get_layer(botm_array, i, j, layer_centers.ravel())
        layer_in_botm_array = np.reshape(layer_in_botm_array, layer_centers.shape)
        zones.append(layer_in_botm_array)
    zones = np.array(zones, dtype=int)
    if is2D:
        zones = zones[0]
    return zones


def setup_model_layers(dem_means_raster, facies_classes_netcdf, framework_rasters, modelgrid,
                       facies_class_variable='fac_a',
                       dem_elevation_units='meters', framework_raster_elevation_units='meters',
                       model_length_units='meters', output_folder='.',
                       framework_unit_names=None):
    """Generate model layering and property zones from voxel-based zone numbers at uniform depths
    and raster surfaces that represent hydrogeologic contacts. The voxel-based zone numbers may
    represent hydrogeologically different sediments, as determined from an airborne electromagnetic survey, 
    and are assumed to have priority over the raster surfaces. The model top is set at the land surface (DEM). 
    Subsequent cell bottoms beneath the land surface are set based on the voxel depths subtracted from the DEM elevations, 
    as sampled at the grid cell centers. Voxel cells without valid data are filled with zone numbers based
    on their position within the raster surfaces of hydrogeologic contacts. For example, a voxel
    cell that lies above all of the raster surfaces would be assigned a zone number equal to the
    highest voxel-based zone number + 3 (an arbitrary gap in the zone numbering between the voxel data
    and the raster surfaces). Similarly, beneath the lowest voxel cell at each location, grid cells are 
    assigned zone values based on their vertical position relative to the raster surfaces. Below the voxel data,
    botm_array and their bottom elevations correspond to the framework_raster surfaces. The lowermost
    (last) raster surface forms the model bottom.
    

    Parameters
    ----------
    dem_means_raster : str (filepath)
        Raster representing the land surface, at the highest resolution being contemplated for the model.
        Usually this is derived by sampling a higher resolution DEM using zonal statistics, taking
        the mean DEM value for each model cell.
    facies_classes_netcdf : str (filepath)
        NetCDF file with the voxel zone data, with a facies_class_variable that contains the zone numbers. 
    framework_rasters : list of strings (filepaths)
        Raster surfaces describing hydrogelogic contacts surrounding the voxel data.
    modelgrid : mfsetup.MFsetupGrid instance
        Modflow-setup grid instance describing the model grid
    facies_class_variable : str, optional
        Variable in facies_classes_netcdf containing the zone information, by default 'fac_a'
    dem_elevation_units : str, optional
        Elevation units of dem_means_raster, by default 'meters'
    framework_raster_elevation_units : str, optional
        Elevation units of the framework_rasters, by default 'meters'
    model_length_units : str, optional
        Length units used in the model, by default 'meters'
    output_folder : str, optional
        Location where results will be saved, by default '.'
    framework_unit_names : list, optional
        Unit names for the framework_rasters, by default None

    Returns
    -------
    botm_array : 3D numpy array
        Array of layer elevations, starting with the model top. 
        (Length equal to the number of botm_array + 1)
    zone_array : 3D numpy array 
        Array of zone numbers, including the values from facies_classes_netcdf,
        followed by a gap of 3 and then the numbers for the framework_rasters.
        For example, if there are 5 zones in facies_classes_netcdf, and 3 raster surfaces
        that intersect cells outside of the voxel data, zones 1-5 would correspond
        to the voxel zones, and zones 8-10 would correspond to hydrogeologic units
        bounded by the framework_rasters.
        
    Notes
    -----
    In addition to the variables returned, the results are written
    to text-array and GeoTiff files within the output_folder.
    
    Two PDFs of figures are also generated- one showing the layering
    in cross sections at regular intervals throughout the model grid, and one
    showing the zones within each layer in plan view.
    """    

    # output folders
    figures_folder = os.path.join(output_folder, 'figures')
    layers_folder = os.path.join(output_folder, 'botm_array')
    zones_folder = os.path.join(output_folder, 'zones')
    # zone rasters are saved to zones_folder/rasters
    output_folders = [figures_folder,
                      layers_folder,
                      zones_folder,
                      os.path.join(zones_folder, 'rasters'),

                      ]
    for folder in output_folders:
        if not os.path.isdir(folder):
            os.makedirs(folder)

    # output PDF with cross sections
    out_xsection_pdf = os.path.join(figures_folder, 'zone_xsections.pdf')
    out_maps_pdf = os.path.join(figures_folder, 'zone_maps.pdf')

    # flopy model grid
    grid = modelgrid

    # read in the facies classes and resisitivity values
    # from framework team
    # and verify that the grid is aligned with the nhg
    ds = xr.load_dataset(facies_classes_netcdf)
    assert point_is_on_nhg(ds.x, ds.y, offset='center')

    # read in the dem and framework rasters and fill locations of no data
    # with the next valid elevation above
    filled_framework_layers = rasters_to_grid(modelgrid,
                                              dem_means_raster, framework_rasters,
                                              dem_elevation_units=dem_elevation_units,
                                              raster_elevation_units=framework_raster_elevation_units,
                                              dest_elevation_units=model_length_units)
    dem_means = filled_framework_layers[0]
    framework_layer_botms_filled = filled_framework_layers[1:]

    # make AEM framework z edge elevations by subtracting off depths from model top
    voxel_z_edges = []
    for depth in ds.zb.values:
        voxel_z_edges.append(dem_means - depth)
    # get the elevations of the voxel centers as well
    voxel_z_centers = []
    for depth in ds.z.values:
        voxel_z_centers.append(dem_means - depth)
    voxel_z_centers = np.array(voxel_z_centers)

    # get the resistivity facies for each model cell (nearest neighbor)
    voxel_array = ds[facies_class_variable].sel(x=grid.xcellcenters[0],
                                                y=grid.ycellcenters[:, 0], method='nearest')
    voxel_zones = np.unique(voxel_array.values[~np.isnan(voxel_array.values)])

    # fill areas of voxel_array that don't have AEM data with zone values
    # representing units in the original MERAS framework
    framework_layer = layers_to_zones(framework_layer_botms_filled, voxel_z_centers)

    # start framework zone numbers after resistivity facies, +3
    framework_zones = framework_layer + np.nanmax(voxel_array.values) + 3
    voxel_array = voxel_array.values
    voxel_array[np.isnan(voxel_array)] = framework_zones[np.isnan(voxel_array)]

    # put the botm_array together
    layers = voxels_to_layers(voxel_array, z_edges=voxel_z_edges, model_top=dem_means,
                              model_botm=framework_layer_botms_filled, no_data_value=0)

    # add framework zones below the voxel botm_array
    nlay, nrow, ncol = layers.shape
    nlay -= 1  # botm_array includes the model top
    framework_layers = list(range(voxel_array.shape[0], nlay))
    first_framework_zone = framework_zones.min()
    n_framwork_zones = framework_layer_botms_filled.shape[0]
    framework_zone_numbers = np.arange(first_framework_zone, first_framework_zone + n_framwork_zones)

    # assume that the last len(framework_layers) are sequential
    framework_zones_below = []
    for zone_number in framework_zone_numbers[-len(framework_layers):]:
        framework_zones_below.append(np.ones((nrow, ncol)) * zone_number)
    framework_zones_below = np.array(framework_zones_below)

    # make an array of all zones, including the resisitivity facies
    # and the zones representing the units in the original meras framework
    zone_array = np.vstack([voxel_array, framework_zones_below])

    # plot cross sections
    framework_unit_labels = dict(zip(framework_zone_numbers.astype(int), framework_unit_names))
    plot_cross_sections(layers, out_xsection_pdf, property_data=zone_array,
                        voxel_start_layer=0, voxel_zones=voxel_zones,
                        cmap='copper',
                        voxel_cmap='viridis', unit_labels=framework_unit_labels)

    # plot maps showing distribution of zones in each layer
    plot_zone_maps(layers, out_maps_pdf, zones=zone_array,
                   voxel_zones=voxel_zones, unit_labels=framework_unit_labels)

    # write out the layer elevation rasters
    write_raster(os.path.join(layers_folder, 'model_top.tif'), layers[0], xul=grid.xul, yul=grid.yul,
                 dx=grid.delc[0], dy=-grid.delr[0], rotation=0., proj_str='epsg:5070',
                 nodata=-9999, fieldname='elevation')
    for i, layer in enumerate(layers[1:]):
        write_raster(os.path.join(layers_folder, 'botm{}.tif'.format(i)),
                     layer, xul=grid.xul, yul=grid.yul,
                     dx=grid.delc[0], dy=-grid.delr[0],
                     rotation=0., crs=grid.crs,
                     nodata=-9999, fieldname='elevation')

    # write out zone arrays (in GeoTiff and text format)
    voxel_start_layer = 0
    voxel_end_layer = voxel_start_layer + voxel_array.shape[0]
    for i, zones2d in enumerate(zone_array):
        np.savetxt(os.path.join(zones_folder, 'res_fac{}.dat'.format(i)),
                   zones2d, fmt='%.0f')
        write_raster(os.path.join(zones_folder, 'rasters/res_fac{}.tif'.format(i)),
                     zones2d, xul=grid.xul, yul=grid.yul,
                     dx=grid.delc[0], dy=-grid.delr[0],
                     rotation=0., crs=grid.crs,
                     nodata=-9999, fieldname='elevation')
    return layers, zone_array
Functions for preprocessing products from the MAP Framework Team Return the botm_array for elevations at i, j locations.      Parameters     ----------     botm_array : 3D numpy array         layer bottom elevations     i : scaler or sequence         row index (zero-based)     j : scaler or sequence         column index     elev : scaler or sequence         elevation (in same units as model)      Returns     -------     k : np.ndarray (1-D) or scalar         zero-based layer index # .tolist() # identify layer botms that are above and below the elevations # layer is the number of botm_array that are above # force elevations below model bottom into bottom layer Plot a single cross section slice      Parameters     ----------     layer_elevations : 3D numpy array         Array of layer elevations, starting with the model top.          (Length equal to the number of botm_array + 1)     property_data : 3D numpy array          Array of zone numbers generated by setup_model_layers.     row : int or slice instance         If a cross section along a row is desired, row should be a integer,         and column should be a slice instance indicating the range of columns to include.         by default, 0.     column : int or slice instance         If a cross section along a column is desired, column should be a integer,         and row should be a slice instance indicating the range of rows to include.         by default, slice(None), which includes all columns.     voxel_start_layer : int, optional         First layer with voxel data, by default 0     voxel_zones : sequence, optional         Zone numbers within property_data that are voxel-based,          by default None     cmap : str, optional         Matplotlib colormap for non-voxel zone numbers, by default 'copper',         to contrast with colormap for voxel-based zone numbers.     voxel_cmap : str, optional         Matplotlib colormap for voxel-based zone numbers, by default 'viridis'.     unit_labels : dict, optional         Dictionary mapping non-voxel zone numbers to hydrogeologic units,          by default None              Returns     -------     ax : matplotlib AxesSubplot instance for figure # cross section code # create meshgrid for rows or columns # along a row # x, z = np.meshgrid(range(ncol), np.array(z_edges)) # x = grid.xcellcenters[row, column] # along a column # x, z = np.meshgrid(range(nrow), np.array(z_edges)) # x = grid.ycellcenters[row, column] # since z is used to define cell edges in the pcolormesh (below) # z cannot be masked or have nan values # set missing data values (outside of the model footprint) in z # to -9999 # pcolormesh will still skip these cells, as they are defined # as no data by the mask for the property array # zero values will result in pcolormesh edges that dip to zero # on the edge of nodata areas # fill these with previous value in either direction # first drop any indices along the edges # then reassign the zero elevations #z = np.ma.masked_where(np.isnan(z), z) # optionally plot a property such as resistivity facies # drop na values # (areas with no voxel data at any depth) #loc = ~np.all(z.mask, axis=0) #x = np.squeeze(x[loc])  # + [x[-1] + 1] #x = np.ma.masked_array(x, mask=~loc) #if not np.any(z): #    return #data = np.squeeze(data[:, loc]) #thicknesses = np.squeeze(thicknesses[:, loc]) # get the min max values for the existing framework # and voxel-based property data # add separate axis for colorbar to control position # plot the layer bottom elevations # make a discrete colors legend Add a legend with a colored rectangle and label for each     discrete value in data.      Parameters     ----------     ax : matplotlib axes object         An axes of the current figure.     data : np.array or list-like         Discrete data represented by legend.     unit_labels : dict, optional         Dictionary of unit labels keyed by integer value in data,          by default None, in which case all values in data         are included with the integer value as the label.     cmap : str, optional         Matplotlib colormap, by default 'copper'     vmin : float, optional         Minimum value of colormap range,          by default None     vmax : float, optional         Maximum value of colormap range, , by default None     **kwargs : keyword arguments to :func:`matplotlib.pyplot.legend`          Return     lg : :func:`matplotlib.pyplot.legend` handle # make a discrete colors legend # now get corresponding colormap value and label for each category and make legend # normalize to range of 0-1 Generate a multi-page PDF of the layer cross sections.      Parameters     ----------     layers : 3D numpy array          Array of layer elevations, starting with the model top.          (Length equal to the number of botm_array + 1)     property_data : 3D numpy array          Array of zone numbers generated by setup_model_layers.     out_pdf : str (filepath)         Filename of multi-page PDF.     voxel_start_layer : int, optional         First layer with voxel data, by default 0     voxel_zones : sequence, optional         Zone numbers within property_data that are voxel-based,          by default None     cmap : str, optional         Matplotlib colormap for non-voxel zone numbers, by default 'copper',         to contrast with colormap for voxel-based zone numbers.     voxel_cmap : str, optional         Matplotlib colormap for voxel-based zone numbers, by default 'viridis'.     unit_labels : dict, optional         Dictionary mapping non-voxel zone numbers to hydrogeologic units,          by default None Generate a multi-page PDF of the zones in each layer in plan view.      Parameters     ----------     layers : 3D numpy array          Array of layer elevations, starting with the model top.          (Length equal to the number of botm_array + 1)     zones : 3D numpy array          Array of zone numbers generated by setup_model_layers.     out_pdf : str (filepath)         Filename of multi-page PDF.     voxel_zones : sequence, optional         Zone numbers within property_data that are voxel-based,          by default None     unit_labels : str, optional         Matplotlib colormap for non-voxel zone numbers, by default 'copper',         to contrast with colormap for voxel-based zone numbers.     voxel_cmap : str, optional         Matplotlib colormap for voxel-based zone numbers, by default 'viridis'.     unit_labels : dict, optional         Dictionary mapping non-voxel zone numbers to hydrogeologic units,          by default None # split up arrays into voxel-based data (AEM) # and layer-based (existing framework) # everything that isn't voxel is in the previous framework # make a discrete colors legend Sample a sequence of rasters onto the i, j locations of a modelgrid,     returning a 3D numpy array of the sampled elevations. Fill places with nodata     using the next valid surface above.      Parameters     ----------     modelgrid : Modflow-setup :class:`~mfsetup.grid.MFsetupGrid` instance         Modflow-setup grid instance describing the model grid     dem : str (filepath)         Raster representing the land surface, at the highest resolution being contemplated for the model.         Usually this is derived by sampling a higher resolution DEM using zonal statistics, taking         the mean DEM value for each model cell.     rasters : list of strings (filepaths)         Raster surfaces describing hydrogelogic contacts surrounding the voxel data.     dem_elevation_units : str, optional         Elevation units of dem_means_raster, by default 'meters'     framework_raster_elevation_units : str, optional         Elevation units of the framework_rasters, by default 'meters'     model_length_units : str, optional         Length units used in the model, by default 'meters'      References     ----------     See the documentation for the :func:`fill_cells_vertically <mfsetup.discretization.fill_cells_vertically>`     function in Modflow-setup for an explanation of the filling process. # convert to model units # convert to model units # fill nans in the sampled original framework elevations # (nans are where a layer surface is absent) # fill the nans with the next surface above # see https://github.com/aleaf/modflow-setup/blob/develop/mfsetup/discretization.py # reset any values above land surface to land surface Map the layers from one grid to zones on another grid. In other words,     given a 3D array of model cell bottom elevations (botm_array), and a second     3D array of cell centers at the same i, j locations as botm_array,     but with different vertical discretization (model_cell_z_centers),     assign each cell in the model_cell_z_centers array a zone value based on its     vertical position within botm array. For example, a cell that lies above all of the     elevations in botm_array would be assigned a zone number of zero; a cell that lies below     the first bottom elevation would be assigned a value of 1.      Parameters     ----------     botm_array : 3D numpy array         Layer bottom elevations     model_cell_z_centers : 3D numpy array         Cell centers for another grid that is aligned with         botm_array along the i and j axis, but has a different vertical         discretation.      Returns     -------     zones : 3D numpy array         Zone numbers indicating the layer in botm_array for each         cell in model_cell_z_centers. # convert from 2D to 3D # for each, i, j location and cell midpoint in the model_cell_z_centers array # get the corresponding layer within the original meras framework # get the i, j locations for each grid cell Generate model layering and property zones from voxel-based zone numbers at uniform depths     and raster surfaces that represent hydrogeologic contacts. The voxel-based zone numbers may     represent hydrogeologically different sediments, as determined from an airborne electromagnetic survey,      and are assumed to have priority over the raster surfaces. The model top is set at the land surface (DEM).      Subsequent cell bottoms beneath the land surface are set based on the voxel depths subtracted from the DEM elevations,      as sampled at the grid cell centers. Voxel cells without valid data are filled with zone numbers based     on their position within the raster surfaces of hydrogeologic contacts. For example, a voxel     cell that lies above all of the raster surfaces would be assigned a zone number equal to the     highest voxel-based zone number + 3 (an arbitrary gap in the zone numbering between the voxel data     and the raster surfaces). Similarly, beneath the lowest voxel cell at each location, grid cells are      assigned zone values based on their vertical position relative to the raster surfaces. Below the voxel data,     botm_array and their bottom elevations correspond to the framework_raster surfaces. The lowermost     (last) raster surface forms the model bottom.           Parameters     ----------     dem_means_raster : str (filepath)         Raster representing the land surface, at the highest resolution being contemplated for the model.         Usually this is derived by sampling a higher resolution DEM using zonal statistics, taking         the mean DEM value for each model cell.     facies_classes_netcdf : str (filepath)         NetCDF file with the voxel zone data, with a facies_class_variable that contains the zone numbers.      framework_rasters : list of strings (filepaths)         Raster surfaces describing hydrogelogic contacts surrounding the voxel data.     modelgrid : mfsetup.MFsetupGrid instance         Modflow-setup grid instance describing the model grid     facies_class_variable : str, optional         Variable in facies_classes_netcdf containing the zone information, by default 'fac_a'     dem_elevation_units : str, optional         Elevation units of dem_means_raster, by default 'meters'     framework_raster_elevation_units : str, optional         Elevation units of the framework_rasters, by default 'meters'     model_length_units : str, optional         Length units used in the model, by default 'meters'     output_folder : str, optional         Location where results will be saved, by default '.'     framework_unit_names : list, optional         Unit names for the framework_rasters, by default None      Returns     -------     botm_array : 3D numpy array         Array of layer elevations, starting with the model top.          (Length equal to the number of botm_array + 1)     zone_array : 3D numpy array          Array of zone numbers, including the values from facies_classes_netcdf,         followed by a gap of 3 and then the numbers for the framework_rasters.         For example, if there are 5 zones in facies_classes_netcdf, and 3 raster surfaces         that intersect cells outside of the voxel data, zones 1-5 would correspond         to the voxel zones, and zones 8-10 would correspond to hydrogeologic units         bounded by the framework_rasters.              Notes     -----     In addition to the variables returned, the results are written     to text-array and GeoTiff files within the output_folder.          Two PDFs of figures are also generated- one showing the layering     in cross sections at regular intervals throughout the model grid, and one     showing the zones within each layer in plan view. # output folders # zone rasters are saved to zones_folder/rasters # output PDF with cross sections # flopy model grid # read in the facies classes and resisitivity values # from framework team # and verify that the grid is aligned with the nhg # read in the dem and framework rasters and fill locations of no data # with the next valid elevation above # make AEM framework z edge elevations by subtracting off depths from model top # get the elevations of the voxel centers as well # get the resistivity facies for each model cell (nearest neighbor) # fill areas of voxel_array that don't have AEM data with zone values # representing units in the original MERAS framework # start framework zone numbers after resistivity facies, +3 # put the botm_array together # add framework zones below the voxel botm_array # botm_array includes the model top # assume that the last len(framework_layers) are sequential # make an array of all zones, including the resisitivity facies # and the zones representing the units in the original meras framework # plot cross sections # plot maps showing distribution of zones in each layer # write out the layer elevation rasters # write out zone arrays (in GeoTiff and text format)