"""
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
[docs]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
[docs]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
[docs]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
[docs]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()
[docs]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()
[docs]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
[docs]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
[docs]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 in range(voxel_start_layer, voxel_end_layer):
np.savetxt(os.path.join(zones_folder, 'res_fac{}.dat'.format(i)),
voxel_array[i], fmt='%.0f')
write_raster(os.path.join(zones_folder, 'rasters/res_fac{}.tif'.format(i)),
voxel_array[i], 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