From 27f80f6b98a1e0d17f96c7a95d7d6255c165a769 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Wed, 20 Mar 2024 17:50:16 -0400 Subject: [PATCH] refactor(modpathfile): step toward generic particle tracking api * introduce particletracking module and ParticleTrackFile base class * rename _ModpathSeries -> ModpathFile * deduplicate shared logic in ModpathFile * add dtypes as public class attributes * clarify canonical (minimal) fields * prep to support PRT output files --- autotest/test_modpathfile.py | 2 +- autotest/test_mp5.py | 12 +- flopy/utils/flopy_io.py | 1 - flopy/utils/modpathfile.py | 1628 ++++++++----------------------- flopy/utils/particletracking.py | 333 +++++++ 5 files changed, 758 insertions(+), 1218 deletions(-) create mode 100644 flopy/utils/particletracking.py diff --git a/autotest/test_modpathfile.py b/autotest/test_modpathfile.py index 32a0531eff..83795759c5 100644 --- a/autotest/test_modpathfile.py +++ b/autotest/test_modpathfile.py @@ -347,7 +347,7 @@ def test_write_shapefile(function_tmpdir, mp7_small, longfieldname): # write the pathline recarray to shapefile pathline_file.write_shapefile( - pathline_data=pathlines, + data=pathlines, shpname=shp_file, one_per_particle=False, mg=grid, diff --git a/autotest/test_mp5.py b/autotest/test_mp5.py index 7c4515444b..a013e2da34 100644 --- a/autotest/test_mp5.py +++ b/autotest/test_mp5.py @@ -1,10 +1,8 @@ import os import numpy as np -import pandas as pd from autotest.test_mp6 import eval_timeseries from matplotlib import pyplot as plt -from modflow_devtools.markers import requires_pkg from flopy.modflow import Modflow from flopy.plot import PlotMapView @@ -50,14 +48,8 @@ def test_mp5_load(function_tmpdir, example_data_path): for n in pthobj.nid: p = pthobj.get_data(partid=n) e = endobj.get_data(partid=n) - try: - mm.plot_pathline(p, colors=colors[n], layer="all") - except: - assert False, f'could not plot pathline {n + 1} with layer="all"' - try: - mm.plot_endpoint(e) - except: - assert False, f'could not plot endpoint {n + 1} with layer="all"' + mm.plot_pathline(p, colors=colors[n], layer="all") + mm.plot_endpoint(e) # plot the grid and ibound array mm.plot_grid(lw=0.5) diff --git a/flopy/utils/flopy_io.py b/flopy/utils/flopy_io.py index 126d47e904..ffd09c5dc7 100644 --- a/flopy/utils/flopy_io.py +++ b/flopy/utils/flopy_io.py @@ -352,7 +352,6 @@ def loadtxt( ra : np.recarray Numpy record array of file contents. """ - from ..utils import import_optional_dependency if use_pandas: if delimiter.isspace(): diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py index 76b8c6ed18..efc53c68eb 100644 --- a/flopy/utils/modpathfile.py +++ b/flopy/utils/modpathfile.py @@ -1,228 +1,106 @@ """ -Module to read MODPATH output files. The module contains two -important classes that can be accessed by the user. - -* EndpointFile (ascii endpoint file) -* PathlineFile (ascii pathline file) - +Support for MODPATH output files. """ import itertools import os -from pathlib import Path -from typing import Union +from typing import List, Union, Tuple, Optional import numpy as np -from numpy.lib.recfunctions import append_fields, repack_fields, stack_arrays - -from ..utils.flopy_io import loadtxt - - -class _ModpathSeries: - """ - Base class for PathlineFile and TimeseriesFile objects. - - This class served to reduce the amount of duplicated code and - increase maintainability of the modpath output methods - - Parameters - ---------- - filename : str or PathLike - Path of pathline or modpath file - verbose : bool - Write information to the screen. Default is False - output_type : str - pathline or timeseries file type - - """ - - def __init__(self, filename, verbose=False, output_type="pathline"): - self.fname = Path(filename).expanduser().absolute() - self.verbose = verbose - self.output_type = output_type.upper() - - self._build_index() - - # set output type - self.outdtype = self._get_outdtype() - - def _build_index(self): - """ - Set position of the start of the pathline data. - """ - compact = False - self.skiprows = 0 - self.file = open(self.fname, "r") - while True: - line = self.file.readline() - if isinstance(line, bytes): - line = line.decode() - if self.skiprows < 1: - if f"MODPATH_{self.output_type}_FILE 6" in line.upper(): - self.version = 6 - elif ( - f"MODPATH_{self.output_type}_FILE 7" - in line.upper() - ): - self.version = 7 - elif "MODPATH 5.0" in line.upper(): - self.version = 5 - if "COMPACT" in line.upper(): - compact = True - elif "MODPATH Version 3.00" in line.upper(): - self.version = 3 - else: - self.version = None - if self.version is None: - errmsg = f"{self.fname} is not a valid {self.output_type.lower()} file" - raise Exception(errmsg) - self.skiprows += 1 - if self.version == 6 or self.version == 7: - if "end header" in line.lower(): - break - elif self.version == 3 or self.version == 5: - break - - # set compact - self.compact = compact - - # return to start of file - self.file.seek(0) - - def _get_outdtype(self): - outdtype = np.dtype( - [ - ("x", np.float32), - ("y", np.float32), - ("z", np.float32), - ("time", np.float32), - ("k", np.int32), - ("particleid", np.int32), - ] - ) - return outdtype - - def get_maxid(self): - """ - Get the maximum timeseries number in the file timeseries file - - Returns - ---------- - out : int - Maximum pathline number. - - """ - return self._data["particleid"].max() - - def get_maxtime(self): - """ - Get the maximum time in pathline file - - Returns - ---------- - out : float - Maximum pathline time. - - """ - return self._data["time"].max() - - def get_data(self, partid=0, totim=None, ge=True): - """ - get pathline data from the pathline file for a single pathline. - - Parameters - ---------- - partid : int - The zero-based particle id - totim : float - The simulation time. - ge : bool - Boolean that determines if pathline times greater than or equal - to or less than or equal to totim. +from numpy.lib.recfunctions import append_fields, repack_fields - Returns - ---------- - ra : np.recarray - Recarray with the x, y, z, time, k, and particleid. +from flopy.utils.particletracking import ParticleTrackFile - """ - ra = self._data - if totim is not None: - if ge: - idx = np.where( - (ra["time"] >= totim) & (ra["particleid"] == partid) - )[0] - else: - idx = np.where( - (ra["time"] <= totim) & (ra["particleid"] == partid) - )[0] - else: - idx = np.where(ra["particleid"] == partid)[0] - ra = ra[idx] - return ra[["x", "y", "z", "time", "k", "particleid"]] - - def get_alldata(self, totim=None, ge=True): - """ - get data from the output file for all particles and all times. - - Parameters - ---------- - totim : float - The simulation time. - ge : bool - Boolean that determines if pathline times greater than or equal - to or less than or equal to totim. +from ..utils.flopy_io import loadtxt - Returns - ---------- - plist : list of numpy record arrays - A list of numpy recarrays - """ - ra = self._data - if totim is not None: - if ge: - idx = np.where(ra["time"] >= totim)[0] - else: - idx = np.where(ra["time"] <= totim)[0] - if len(idx) > 0: - ra = ra[idx] - ra = ra[["x", "y", "z", "time", "k", "particleid"]] - return [ra[ra["particleid"] == i] for i in range(self.nid.size)] +class ModpathFile(ParticleTrackFile): + """Provides MODPATH output file support.""" - def get_destination_data(self, dest_cells, to_recarray=True): - """ - Get data for set of destination cells. + def __init__( + self, filename: Union[str, os.PathLike], verbose: bool = False + ): + super().__init__(filename, verbose) + self.output_type = self.__class__.__name__.lower().replace("file", "") + ( + self.modpath, + self.compact, + self.skiprows, + self.version, + self.direction, + ) = self.parse(filename, self.output_type) + + @staticmethod + def parse( + file_path: Union[str, os.PathLike], file_type: str + ) -> Tuple[bool, int, int, Optional[int]]: + """ + Extract preliminary information from a MODPATH output file: + - whether in compact format + - how many rows to skip + - the MODPATH version Parameters ---------- - dest_cells : list or array of tuples - (k, i, j) of each destination cell for MODPATH versions less than - MODPATH 7 or node number of each destination cell. (zero based) - to_recarray : bool - Boolean that controls returned series. If to_recarray is True, - a single recarray with all of the pathlines that intersect - dest_cells are returned. If to_recarray is False, a list of - recarrays (the same form as returned by get_alldata method) - that intersect dest_cells are returned (default is False). + file_path : str or PathLike + The output file path + file_type : str + The output file type: pathline, endpoint, or timeseries Returns ------- - series : np.recarray - Slice of data array (e.g. PathlineFile._data, TimeseriesFile._data) - containing endpoint, pathline, or timeseries data that intersect - (k,i,j) or (node) dest_cells. - + out : bool, int, int + A tuple (compact, skiprows, version) """ - # create local copy of _data - ra = np.array(self._data) + modpath = True + compact = False + idx = 0 + skiprows = 0 + version = None + direction = None + with open(file_path, "r") as f: + while True: + line = f.readline() + if isinstance(line, bytes): + line = line.decode() + if skiprows < 1: + if f"MODPATH_{file_type.upper()}_FILE 6" in line.upper(): + version = 6 + elif ( + f"MODPATH_{file_type.upper()}_FILE 7" + in line.upper() + ): + version = 7 + elif "MODPATH 5.0" in line.upper(): + version = 5 + if "COMPACT" in line.upper(): + compact = True + elif "MODPATH Version 3.00" in line.upper(): + version = 3 + if version is None: + modpath = False + skiprows += 1 + if version in [6, 7]: + if file_type.lower() == "endpoint": + if idx == 1: + direction = 1 + if int(line.strip()[0]) == 2: + direction = -1 + idx += 1 + if "end header" in line.lower(): + break + else: + break - # find the intersection of pathlines and dest_cells - # convert dest_cells to same dtype for comparison + return modpath, compact, skiprows, version, direction + + def intersect( + self, cells, to_recarray + ) -> Union[List[np.recarray], np.recarray]: if self.version < 7: try: - raslice = ra[["k", "i", "j"]] + raslice = self._data[["k", "i", "j"]] except (KeyError, ValueError): raise KeyError( "could not extract 'k', 'i', and 'j' keys " @@ -230,258 +108,63 @@ def get_destination_data(self, dest_cells, to_recarray=True): ) else: try: - raslice = ra[["node"]] + raslice = self._data[["node"]] except (KeyError, ValueError): msg = "could not extract 'node' key from {} data".format( self.output_type.lower() ) raise KeyError(msg) - if isinstance(dest_cells, (list, tuple)): - allint = all(isinstance(el, int) for el in dest_cells) + if isinstance(cells, (list, tuple)): + allint = all(isinstance(el, int) for el in cells) # convert to a list of tuples if allint: t = [] - for el in dest_cells: + for el in cells: t.append((el,)) - dest_cells = t + cells = t - dest_cells = np.array(dest_cells, dtype=raslice.dtype) - inds = np.in1d(raslice, dest_cells) - epdest = ra[inds].copy().view(np.recarray) + cells = np.array(cells, dtype=raslice.dtype) + inds = np.in1d(raslice, cells) + epdest = self._data[inds].copy().view(np.recarray) if to_recarray: # use particle ids to get the rest of the paths - inds = np.in1d(ra["particleid"], epdest.particleid) - series = ra[inds].copy() + inds = np.in1d(self._data["particleid"], epdest.particleid) + series = self._data[inds].copy() series.sort(order=["particleid", "time"]) series = series.view(np.recarray) else: - # get list of unique particleids in selection + # collect unique particleids in selection partids = np.unique(epdest["particleid"]) - - # build list of unique particleids in selection series = [self.get_data(partid) for partid in partids] return series - def write_shapefile( - self, - data=None, - one_per_particle=True, - direction="ending", - shpname="endpoints.shp", - mg=None, - crs=None, - **kwargs, - ): - """ - Write pathlines or timeseries to a shapefile - - data : np.recarray - Record array of same form as that returned by - get_alldata(). (if none, get_alldata() is exported). - one_per_particle : boolean (default True) - True writes a single LineString with a single set of attribute - data for each particle. False writes a record/geometry for each - pathline segment (each row in the PathLine file). This option can - be used to visualize attribute information (time, model layer, - etc.) across a pathline in a GIS. - direction : str - String defining if starting or ending particle locations should be - included in shapefile attribute information. Only used if - one_per_particle=False. (default is 'ending') - shpname : str - File path for shapefile - mg : flopy.discretization.grid instance - Used to scale and rotate Global x,y,z values. - crs : pyproj.CRS, int, str, optional - Coordinate reference system (CRS) for the model grid - (must be projected; geographic CRS are not supported). - The value can be anything accepted by - :meth:`pyproj.CRS.from_user_input() `, - such as an authority string (eg "EPSG:26916") or a WKT string. - kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp - - .. deprecated:: 3.5 - The following keyword options will be removed for FloPy 3.6: - - - ``epsg`` (int): use ``crs`` instead. - """ - from ..discretization import StructuredGrid - from ..export.shapefile_utils import recarray2shp - from ..utils import geometry - from ..utils.geometry import LineString - - series = data - if series is None: - series = self._data.view(np.recarray) - else: - # convert pathline list to a single recarray - if isinstance(series, list): - s = series[0] - print(s.dtype) - for n in range(1, len(series)): - s = stack_arrays((s, series[n])) - series = s.view(np.recarray) - - series = series.copy() - series.sort(order=["particleid", "time"]) - - if mg is None: - raise ValueError("A modelgrid object was not provided.") - - particles = np.unique(series.particleid) - geoms = [] - - # create dtype with select attributes in pth - names = series.dtype.names - dtype = [] - atts = ["particleid", "particlegroup", "time", "k", "i", "j", "node"] - for att in atts: - if att in names: - t = np.int32 - if att == "time": - t = np.float32 - dtype.append((att, t)) - dtype = np.dtype(dtype) - - # reset names to the selected names in the created dtype - names = dtype.names - - # 1 geometry for each path - if one_per_particle: - loc_inds = 0 - if direction == "ending": - loc_inds = -1 - - sdata = [] - for pid in particles: - ra = series[series.particleid == pid] - - x, y = geometry.transform( - ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians - ) - z = ra.z - geoms.append(LineString(list(zip(x, y, z)))) - - t = [pid] - if "particlegroup" in names: - t.append(ra.particlegroup[0]) - t.append(ra.time.max()) - if "node" in names: - t.append(ra.node[loc_inds]) - else: - if "k" in names: - t.append(ra.k[loc_inds]) - if "i" in names: - t.append(ra.i[loc_inds]) - if "j" in names: - t.append(ra.j[loc_inds]) - sdata.append(tuple(t)) - - sdata = np.array(sdata, dtype=dtype).view(np.recarray) - - # geometry for each row in PathLine file - else: - dtype = series.dtype - sdata = [] - for pid in particles: - ra = series[series.particleid == pid] - if mg is not None: - x, y = geometry.transform( - ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians - ) - else: - x, y = geometry.transform(ra.x, ra.y, 0, 0, 0) - z = ra.z - geoms += [ - LineString( - [(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])] - ) - for i in np.arange(1, (len(ra))) - ] - sdata += ra[1:].tolist() - sdata = np.array(sdata, dtype=dtype).view(np.recarray) - - # convert back to one-based - for n in set(self.kijnames).intersection(set(sdata.dtype.names)): - sdata[n] += 1 - - # write the final recarray to a shapefile - recarray2shp(sdata, geoms, shpname=shpname, crs=crs, **kwargs) - - -class PathlineFile(_ModpathSeries): +class PathlineFile(ModpathFile): """ - PathlineFile Class. + Particle pathline file. Parameters ---------- filename : str or PathLike Path of the pathline file verbose : bool - Write information to the screen. Default is False. + Show verbose output. Default is False. Examples -------- >>> import flopy - >>> pthobj = flopy.utils.PathlineFile('model.mppth') - >>> p1 = pthobj.get_data(partid=1) + >>> pl_file = flopy.utils.PathlineFile('model.mppth') + >>> pl1 = pthobj.get_data(partid=1) """ - kijnames = [ - "k", - "i", - "j", - "node", - "particleid", - "particlegroup", - "linesegmentindex", - "particleidloc", - "sequencenumber", - ] - - def __init__(self, filename: Union[str, os.PathLike], verbose=False): - """ - Class constructor. - - """ - - super().__init__(filename, verbose=verbose, output_type="pathline") - - # set data dtype and read pathline data - if self.version == 7: - self.dtype, self._data = self._get_mp7data() - else: - self.dtype = self._get_dtypes() - self._data = loadtxt( - self.file, dtype=self.dtype, skiprows=self.skiprows - ) - - # convert layer, row, and column indices; particle id and group; and - # line segment indices to zero-based - for n in self.kijnames: - if n in self._data.dtype.names: - self._data[n] -= 1 - - # set number of particle ids - self.nid = np.unique(self._data["particleid"]) - - # sort data - self._data.sort(order=["particleid", "time"]) - - # close the input file - self.file.close() - - def _get_dtypes(self): - """ - Build numpy dtype for the MODPATH 6 pathline file. - """ - if self.version == 3 or self.version == 5: - dtype = np.dtype( + dtypes = { + **dict.fromkeys( + [3, 5], + np.dtype( [ ("particleid", np.int32), ("x", np.float32), @@ -494,52 +177,29 @@ def _get_dtypes(self): ("k", np.int32), ("cumulativetimestep", np.int32), ] - ) - elif self.version == 6: - dtype = np.dtype( - [ - ("particleid", np.int32), - ("particlegroup", np.int32), - ("timepointindex", np.int32), - ("cumulativetimestep", np.int32), - ("time", np.float32), - ("x", np.float32), - ("y", np.float32), - ("z", np.float32), - ("k", np.int32), - ("i", np.int32), - ("j", np.int32), - ("grid", np.int32), - ("xloc", np.float32), - ("yloc", np.float32), - ("zloc", np.float32), - ("linesegmentindex", np.int32), - ] - ) - elif self.version == 7: - raise TypeError( - "_get_dtypes() should not be called for " - "MODPATH 7 pathline files" - ) - return dtype - - def _get_mp7data(self): - dtyper = np.dtype( + ), + ), + 6: np.dtype( [ - ("node", np.int32), + ("particleid", np.int32), + ("particlegroup", np.int32), + ("timepointindex", np.int32), + ("cumulativetimestep", np.int32), + ("time", np.float32), ("x", np.float32), ("y", np.float32), ("z", np.float32), - ("time", np.float32), + ("k", np.int32), + ("i", np.int32), + ("j", np.int32), + ("grid", np.int32), ("xloc", np.float32), ("yloc", np.float32), ("zloc", np.float32), - ("k", np.int32), - ("stressperiod", np.int32), - ("timestep", np.int32), + ("linesegmentindex", np.int32), ] - ) - dtype = np.dtype( + ), + 7: np.dtype( [ ("particleid", np.int32), ("particlegroup", np.int32), @@ -557,152 +217,115 @@ def _get_mp7data(self): ("stressperiod", np.int32), ("timestep", np.int32), ] - ) - idx = 0 - part_dict = {} - ndata = 0 - while True: - if idx == 0: - for n in range(self.skiprows): - line = self.file.readline() - # read header line - try: - line = self.file.readline().strip() - if self.verbose: - print(line) - if len(line) < 1: - break - except: - break - t = [int(s) for j, s in enumerate(line.split()) if j < 4] - sequencenumber, group, particleid, pathlinecount = t[0:4] - ndata += pathlinecount - # read in the particle data - d = np.loadtxt( - itertools.islice(self.file, 0, pathlinecount), dtype=dtyper - ) - key = (idx, sequencenumber, group, particleid, pathlinecount) - part_dict[key] = d.copy() - idx += 1 - - # create data array - data = np.zeros(ndata, dtype=dtype) - - # fill data - ipos0 = 0 - for key, value in part_dict.items(): - idx, sequencenumber, group, particleid, pathlinecount = key[0:5] - ipos1 = ipos0 + pathlinecount - # fill constant items for particle - # particleid is not necessarily unique for all pathlines - use - # sequencenumber which is unique - data["particleid"][ipos0:ipos1] = sequencenumber - # set particlegroup and sequence number - data["particlegroup"][ipos0:ipos1] = group - data["sequencenumber"][ipos0:ipos1] = sequencenumber - # save particleidloc to particleid - data["particleidloc"][ipos0:ipos1] = particleid - # fill particle data - for name in value.dtype.names: - data[name][ipos0:ipos1] = value[name] - ipos0 = ipos1 + ), + } - return dtype, data - - def get_maxid(self): - """ - Get the maximum pathline number in the file pathline file - - Returns - ---------- - out : int - Maximum pathline number. - - """ - return super().get_maxid() - - def get_maxtime(self): - """ - Get the maximum time in pathline file - - Returns - ---------- - out : float - Maximum pathline time. - - """ - return super().get_maxtime() - - def get_data(self, partid=0, totim=None, ge=True): - """ - get pathline data from the pathline file for a single pathline. - - Parameters - ---------- - partid : int - The zero-based particle id. The first record is record 0. - totim : float - The simulation time. All pathline points for particle partid - that are greater than or equal to (ge=True) or less than or - equal to (ge=False) totim will be returned. Default is None - ge : bool - Boolean that determines if pathline times greater than or equal - to or less than or equal to totim is used to create a subset - of pathlines. Default is True. - - Returns - ---------- - ra : numpy record array - A numpy recarray with the x, y, z, time, k, and particleid for - pathline partid. - - - See Also - -------- - - Notes - ----- - - Examples - -------- - - >>> import flopy.utils.modpathfile as mpf - >>> pthobj = flopy.utils.PathlineFile('model.mppth') - >>> p1 = pthobj.get_data(partid=1) + kijnames = [ + "k", + "i", + "j", + "node", + "particleid", + "particlegroup", + "linesegmentindex", + "particleidloc", + "sequencenumber", + ] - """ - return super().get_data(partid=partid, totim=totim, ge=ge) + def __init__( + self, filename: Union[str, os.PathLike], verbose: bool = False + ): + super().__init__(filename, verbose=verbose) + self.dtype, self._data = self._load() + self.nid = np.unique(self._data["particleid"]) - def get_alldata(self, totim=None, ge=True): - """ - get pathline data from the pathline file for all pathlines and all times. + def _load(self) -> Tuple[np.dtype, np.ndarray]: + dtype = self.dtypes[self.version] + if self.version == 7: + dtyper = np.dtype( + [ + ("node", np.int32), + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("time", np.float32), + ("xloc", np.float32), + ("yloc", np.float32), + ("zloc", np.float32), + ("k", np.int32), + ("stressperiod", np.int32), + ("timestep", np.int32), + ] + ) + idx = 0 + particle_pathlines = {} + nrows = 0 + with open(self.fname) as f: + while True: + if idx == 0: + for n in range(self.skiprows): + line = f.readline() + # read header line + try: + line = f.readline().strip() + if self.verbose: + print(line) + if len(line) < 1: + break + except: + break + t = [int(s) for j, s in enumerate(line.split()) if j < 4] + sequencenumber, group, particleid, pathlinecount = t[0:4] + nrows += pathlinecount + # read in the particle data + d = np.loadtxt( + itertools.islice(f, 0, pathlinecount), dtype=dtyper + ) + key = ( + idx, + sequencenumber, + group, + particleid, + pathlinecount, + ) + particle_pathlines[key] = d.copy() + idx += 1 - Parameters - ---------- - totim : float - The simulation time. All pathline points for particle partid - that are greater than or equal to (ge=True) or less than or - equal to (ge=False) totim will be returned. Default is None - ge : bool - Boolean that determines if pathline times greater than or equal - to or less than or equal to totim is used to create a subset - of pathlines. Default is True. + # create data array + data = np.zeros(nrows, dtype=dtype) - Returns - ---------- - plist : a list of numpy record array - A list of numpy recarrays with the x, y, z, time, k, and particleid - for all pathlines. + # fill data + ipos0 = 0 + for key, value in particle_pathlines.items(): + idx, sequencenumber, group, particleid, pathlinecount = key[ + 0:5 + ] + ipos1 = ipos0 + pathlinecount + # fill constant items for particle + # particleid is not necessarily unique for all pathlines - use + # sequencenumber which is unique + data["particleid"][ipos0:ipos1] = sequencenumber + # set particlegroup and sequence number + data["particlegroup"][ipos0:ipos1] = group + data["sequencenumber"][ipos0:ipos1] = sequencenumber + # save particleidloc to particleid + data["particleidloc"][ipos0:ipos1] = particleid + # fill particle data + for name in value.dtype.names: + data[name][ipos0:ipos1] = value[name] + ipos0 = ipos1 + else: + data = loadtxt(self.fname, dtype=dtype, skiprows=self.skiprows) - Examples - -------- + # convert indices to zero-based + for n in self.kijnames: + if n in data.dtype.names: + data[n] -= 1 - >>> import flopy.utils.modpathfile as mpf - >>> pthobj = flopy.utils.PathlineFile('model.mppth') - >>> p = pthobj.get_alldata() + # sort by particle ID and time + data.sort(order=["particleid", "time"]) - """ - return super().get_alldata(totim=totim, ge=ge) + return dtype, data def get_destination_pathline_data(self, dest_cells, to_recarray=False): """ @@ -722,7 +345,7 @@ def get_destination_pathline_data(self, dest_cells, to_recarray=False): Returns ------- - pthldest : np.recarray + data : np.recarray Slice of pathline data array (e.g. PathlineFile._data) containing only pathlines that pass through (k,i,j) or (node) dest_cells. @@ -740,84 +363,120 @@ def get_destination_pathline_data(self, dest_cells, to_recarray=False): dest_cells=dest_cells, to_recarray=to_recarray ) - def write_shapefile( - self, - pathline_data=None, - one_per_particle=True, - direction="ending", - shpname="pathlines.shp", - mg=None, - crs=None, - **kwargs, - ): - """ - Write pathlines to a shapefile - - pathline_data : np.recarray - Record array of same form as that returned by - PathlineFile.get_alldata(). (if none, PathlineFile.get_alldata() - is exported). - one_per_particle : boolean (default True) - True writes a single LineString with a single set of attribute - data for each particle. False writes a record/geometry for each - pathline segment (each row in the PathLine file). This option can - be used to visualize attribute information (time, model layer, - etc.) across a pathline in a GIS. - direction : str - String defining if starting or ending particle locations should be - included in shapefile attribute information. Only used if - one_per_particle=False. (default is 'ending') - shpname : str - File path for shapefile - mg : flopy.discretization.grid instance - Used to scale and rotate Global x,y,z values in MODPATH Pathline - file. - crs : pyproj.CRS, int, str, optional - Coordinate reference system (CRS) for the model grid - (must be projected; geographic CRS are not supported). - The value can be anything accepted by - :meth:`pyproj.CRS.from_user_input() `, - such as an authority string (eg "EPSG:26916") or a WKT string. - kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp - .. deprecated:: 3.5 - The following keyword options will be removed for FloPy 3.6: - - - ``epsg`` (int): use ``crs`` instead. - - """ - super().write_shapefile( - data=pathline_data, - one_per_particle=one_per_particle, - direction=direction, - shpname=shpname, - mg=mg, - crs=crs, - **kwargs, - ) - - -class EndpointFile: +class EndpointFile(ModpathFile): """ - EndpointFile Class. + Particle endpoint file. Parameters ---------- - filename : string - Name of the endpoint file + filename : str or PathLike + Path of the endpoint file verbose : bool - Write information to the screen. Default is False. + Show verbose output. Default is False. Examples -------- >>> import flopy - >>> endobj = flopy.utils.EndpointFile('model.mpend') - >>> e1 = endobj.get_data(partid=1) - + >>> ep_file = flopy.utils.EndpointFile('model.mpend') + >>> ep1 = endobj.get_data(partid=1) """ + dtypes = { + **dict.fromkeys( + [3, 5], + np.dtype( + [ + ("zone", np.int32), + ("j", np.int32), + ("i", np.int32), + ("k", np.int32), + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("zloc", np.float32), + ("time", np.float32), + ("x0", np.float32), + ("y0", np.float32), + ("zloc0", np.float32), + ("j0", np.int32), + ("i0", np.int32), + ("k0", np.int32), + ("zone0", np.int32), + ("cumulativetimestep", np.int32), + ("ipcode", np.int32), + ("time0", np.float32), + ] + ), + ), + 6: np.dtype( + [ + ("particleid", np.int32), + ("particlegroup", np.int32), + ("status", np.int32), + ("time0", np.float32), + ("time", np.float32), + ("initialgrid", np.int32), + ("k0", np.int32), + ("i0", np.int32), + ("j0", np.int32), + ("cellface0", np.int32), + ("zone0", np.int32), + ("xloc0", np.float32), + ("yloc0", np.float32), + ("zloc0", np.float32), + ("x0", np.float32), + ("y0", np.float32), + ("z0", np.float32), + ("finalgrid", np.int32), + ("k", np.int32), + ("i", np.int32), + ("j", np.int32), + ("cellface", np.int32), + ("zone", np.int32), + ("xloc", np.float32), + ("yloc", np.float32), + ("zloc", np.float32), + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("label", "|S40"), + ] + ), + 7: np.dtype( + [ + ("particleid", np.int32), + ("particlegroup", np.int32), + ("particleidloc", np.int32), + ("status", np.int32), + ("time0", np.float32), + ("time", np.float32), + ("node0", np.int32), + ("k0", np.int32), + ("xloc0", np.float32), + ("yloc0", np.float32), + ("zloc0", np.float32), + ("x0", np.float32), + ("y0", np.float32), + ("z0", np.float32), + ("zone0", np.int32), + ("initialcellface", np.int32), + ("node", np.int32), + ("k", np.int32), + ("xloc", np.float32), + ("yloc", np.float32), + ("zloc", np.float32), + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("zone", np.int32), + ("cellface", np.int32), + ] + ), + } + kijnames = [ "k0", "i0", @@ -834,262 +493,42 @@ class EndpointFile: "zone", ] - def __init__(self, filename, verbose=False): - """ - Class constructor. + def __init__( + self, filename: Union[str, os.PathLike], verbose: bool = False + ): + super().__init__(filename, verbose) + self.dtype, self._data = self._load() + self.nid = np.unique(self._data["particleid"]) - """ - self.fname = filename - self.verbose = verbose - self._build_index() - self.dtype = self._get_dtypes() - self._data = loadtxt( - self.file, dtype=self.dtype, skiprows=self.skiprows - ) - # add particleid if required - self._add_particleid() + def _load(self) -> Tuple[np.dtype, np.ndarray]: + dtype = self.dtypes[self.version] + data = loadtxt(self.fname, dtype=dtype, skiprows=self.skiprows) - # convert layer, row, and column indices; particle id and group; and - # line segment indices to zero-based + # convert indices to zero-based for n in self.kijnames: - if n in self._data.dtype.names: - self._data[n] -= 1 - - # set number of particle ids - self.nid = np.unique(self._data["particleid"]).shape[0] - - # close the input file - self.file.close() - return - - def _build_index(self): - """ - Set position of the start of the pathline data. - """ - self.skiprows = 0 - self.file = open(self.fname, "r") - idx = 0 - while True: - line = self.file.readline() - if isinstance(line, bytes): - line = line.decode() - if self.skiprows < 1: - if "MODPATH_ENDPOINT_FILE 6" in line.upper(): - self.version = 6 - elif "MODPATH_ENDPOINT_FILE 7" in line.upper(): - self.version = 7 - elif "MODPATH 5.0" in line.upper(): - self.version = 5 - elif "MODPATH Version 3.00" in line.upper(): - self.version = 3 - else: - self.version = None - if self.version is None: - errmsg = f"{self.fname} is not a valid endpoint file" - raise Exception(errmsg) - self.skiprows += 1 - if self.version == 6 or self.version == 7: - if idx == 1: - t = line.strip() - self.direction = 1 - if int(t[0]) == 2: - self.direction = -1 - idx += 1 - if "end header" in line.lower(): - break - else: - break - self.file.seek(0) + if n in data.dtype.names: + data[n] -= 1 - if self.verbose: - print(f"MODPATH version {self.version} endpoint file") - - def _get_dtypes(self): - """ - Build numpy dtype for the MODPATH 6 endpoint file. - """ - if self.version == 3 or self.version == 5: - dtype = self._get_mp35_dtype() - elif self.version == 6: - dtype = self._get_mp6_dtype() - elif self.version == 7: - dtype = self._get_mp7_dtype() - return dtype - - def _get_mp35_dtype(self, add_id=False): - dtype = [ - ("zone", np.int32), - ("j", np.int32), - ("i", np.int32), - ("k", np.int32), - ("x", np.float32), - ("y", np.float32), - ("z", np.float32), - ("zloc", np.float32), - ("time", np.float32), - ("x0", np.float32), - ("y0", np.float32), - ("zloc0", np.float32), - ("j0", np.int32), - ("i0", np.int32), - ("k0", np.int32), - ("zone0", np.int32), - ("cumulativetimestep", np.int32), - ("ipcode", np.int32), - ("time0", np.float32), - ] - if add_id: - dtype.insert(0, ("particleid", np.int32)) - return np.dtype(dtype) - - def _get_mp6_dtype(self): - dtype = [ - ("particleid", np.int32), - ("particlegroup", np.int32), - ("status", np.int32), - ("time0", np.float32), - ("time", np.float32), - ("initialgrid", np.int32), - ("k0", np.int32), - ("i0", np.int32), - ("j0", np.int32), - ("cellface0", np.int32), - ("zone0", np.int32), - ("xloc0", np.float32), - ("yloc0", np.float32), - ("zloc0", np.float32), - ("x0", np.float32), - ("y0", np.float32), - ("z0", np.float32), - ("finalgrid", np.int32), - ("k", np.int32), - ("i", np.int32), - ("j", np.int32), - ("cellface", np.int32), - ("zone", np.int32), - ("xloc", np.float32), - ("yloc", np.float32), - ("zloc", np.float32), - ("x", np.float32), - ("y", np.float32), - ("z", np.float32), - ("label", "|S40"), - ] - return np.dtype(dtype) - - def _get_mp7_dtype(self): - dtype = [ - ("particleid", np.int32), - ("particlegroup", np.int32), - ("particleidloc", np.int32), - ("status", np.int32), - ("time0", np.float32), - ("time", np.float32), - ("node0", np.int32), - ("k0", np.int32), - ("xloc0", np.float32), - ("yloc0", np.float32), - ("zloc0", np.float32), - ("x0", np.float32), - ("y0", np.float32), - ("z0", np.float32), - ("zone0", np.int32), - ("initialcellface", np.int32), - ("node", np.int32), - ("k", np.int32), - ("xloc", np.float32), - ("yloc", np.float32), - ("zloc", np.float32), - ("x", np.float32), - ("y", np.float32), - ("z", np.float32), - ("zone", np.int32), - ("cellface", np.int32), - ] - return np.dtype(dtype) - - def _add_particleid(self): # add particle ids for earlier version of MODPATH if self.version < 6: - # create particle ids - shaped = self._data.shape[0] - pids = np.arange(1, shaped + 1, 1, dtype=np.int32) - - # for numpy version 1.14 and higher - self._data = append_fields(self._data, "particleid", pids) - return - - def get_maxid(self): - """ - Get the maximum endpoint particle id in the file endpoint file - - Returns - ---------- - out : int - Maximum endpoint particle id. - - """ - return np.unique(self._data["particleid"]).max() - - def get_maxtime(self): - """ - Get the maximum time in the endpoint file - - Returns - ---------- - out : float - Maximum endpoint time. + shape = data.shape[0] + pids = np.arange(1, shape + 1, 1, dtype=np.int32) + data = append_fields(data, "particleid", pids) - """ - return self._data["time"].max() + return dtype, data def get_maxtraveltime(self): """ - Get the maximum travel time in the endpoint file + Get the maximum travel time. Returns ---------- out : float - Maximum endpoint travel time. + Maximum travel time. """ return (self._data["time"] - self._data["time0"]).max() - def get_data(self, partid=0): - """ - Get endpoint data from the endpoint file for a single particle. - - Parameters - ---------- - partid : int - The zero-based particle id. The first record is record 0. - (default is 0) - - Returns - ---------- - ra : numpy record array - A numpy recarray with the endpoint particle data for - endpoint partid. - - - See Also - -------- - - Notes - ----- - - Examples - -------- - - >>> import flopy - >>> endobj = flopy.utils.EndpointFile('model.mpend') - >>> e1 = endobj.get_data(partid=1) - - """ - idx = self._data["particleid"] == partid - ra = self._data[idx] - return ra - def get_alldata(self): """ Get endpoint data from the endpoint file for all endpoints. @@ -1099,7 +538,7 @@ def get_alldata(self): Returns ---------- - ra : numpy record array + data : numpy record array A numpy recarray with the endpoint particle data @@ -1134,7 +573,7 @@ def get_destination_endpoint_data(self, dest_cells, source=False): Returns ------- - epdest : np.recarray + data : np.recarray Slice of endpoint data array (e.g. EndpointFile.get_alldata) containing only endpoint data with final locations in (k,i,j) or (node) dest_cells. @@ -1150,7 +589,7 @@ def get_destination_endpoint_data(self, dest_cells, source=False): """ # create local copy of _data - ra = self.get_alldata() + data = self.get_alldata() # find the intersection of endpoints and dest_cells # convert dest_cells to same dtype for comparison @@ -1160,7 +599,7 @@ def get_destination_endpoint_data(self, dest_cells, source=False): else: keys = ["k", "i", "j"] try: - raslice = repack_fields(ra[keys]) + raslice = repack_fields(data[keys]) except (KeyError, ValueError): raise KeyError( "could not extract " @@ -1173,7 +612,7 @@ def get_destination_endpoint_data(self, dest_cells, source=False): else: keys = ["node"] try: - raslice = repack_fields(ra[keys]) + raslice = repack_fields(data[keys]) except (KeyError, ValueError): msg = f"could not extract '{keys[0]}' key from endpoint data" raise KeyError(msg) @@ -1192,8 +631,7 @@ def get_destination_endpoint_data(self, dest_cells, source=False): dest_cells = np.array(dest_cells, dtype=dtype) inds = np.in1d(raslice, dest_cells) - epdest = ra[inds].copy().view(np.recarray) - return epdest + return data[inds].copy().view(np.recarray) def write_shapefile( self, @@ -1273,25 +711,81 @@ def write_shapefile( recarray2shp(epd, geoms, shpname=shpname, crs=crs, **kwargs) -class TimeseriesFile(_ModpathSeries): +class TimeseriesFile(ModpathFile): """ - TimeseriesFile Class. + Particle timeseries file. Parameters ---------- - filename : string - Name of the timeseries file + filename : str or PathLike + Path of the timeseries file verbose : bool - Write information to the screen. Default is False. + Show verbose output. Default is False. Examples -------- >>> import flopy - >>> tsobj = flopy.utils.TimeseriesFile('model.timeseries') + >>> ts_file = flopy.utils.TimeseriesFile('model.timeseries') >>> ts1 = tsobj.get_data(partid=1) """ + dtypes = { + **dict.fromkeys( + [3, 5], + np.dtype( + [ + ("timestepindex", np.int32), + ("particleid", np.int32), + ("node", np.int32), + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("zloc", np.float32), + ("time", np.float32), + ("timestep", np.int32), + ] + ), + ), + 6: np.dtype( + [ + ("timepointindex", np.int32), + ("timestep", np.int32), + ("time", np.float32), + ("particleid", np.int32), + ("particlegroup", np.int32), + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("grid", np.int32), + ("k", np.int32), + ("i", np.int32), + ("j", np.int32), + ("xloc", np.float32), + ("yloc", np.float32), + ("zloc", np.float32), + ] + ), + 7: np.dtype( + [ + ("timepointindex", np.int32), + ("timestep", np.int32), + ("time", np.float32), + ("particleid", np.int32), + ("particlegroup", np.int32), + ("particleidloc", np.int32), + ("node", np.int32), + ("xloc", np.float32), + ("yloc", np.float32), + ("zloc", np.float32), + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("k", np.int32), + ] + ), + } + kijnames = [ "k", "i", @@ -1306,262 +800,40 @@ class TimeseriesFile(_ModpathSeries): ] def __init__(self, filename, verbose=False): - """ - Class constructor. - - """ - super().__init__(filename, verbose=verbose, output_type="timeseries") - - # set dtype - self.dtype = self._get_dtypes() - - # read data - self._data = loadtxt( - self.file, dtype=self.dtype, skiprows=self.skiprows - ) - - # convert layer, row, and column indices; particle id and group; and - # line segment indices to zero-based - for n in self.kijnames: - if n in self._data.dtype.names: - self._data[n] -= 1 - - # set number of particle ids + super().__init__(filename, verbose) + self.dtype, self._data = self._load() self.nid = np.unique(self._data["particleid"]) - # sort data - self._data.sort(order=["particleid", "time"]) - - # close the input file - self.file.close() - return - - def _build_index(self): - """ - Set position of the start of the timeseries data. - """ - compact = False - self.skiprows = 0 - self.file = open(self.fname, "r") - while True: - line = self.file.readline() - if isinstance(line, bytes): - line = line.decode() - if self.skiprows < 1: - if "MODPATH_TIMESERIES_FILE 6" in line.upper(): - self.version = 6 - elif "MODPATH_TIMESERIES_FILE 7" in line.upper(): - self.version = 7 - elif "MODPATH 5.0" in line.upper(): - self.version = 5 - if "COMPACT" in line.upper(): - compact = True - elif "MODPATH Version 3.00" in line.upper(): - self.version = 3 - else: - self.version = None - if self.version is None: - raise Exception( - f"{self.fname} is not a valid timeseries file" - ) - self.skiprows += 1 - if self.version == 6 or self.version == 7: - if "end header" in line.lower(): - break - elif self.version == 3 or self.version == 5: - break - - # set compact - self.compact = compact - - # return to the top of the file - self.file.seek(0) - - def _get_dtypes(self): - """ - Build numpy dtype for the MODPATH 6 timeseries file. - """ - if self.version == 3 or self.version == 5: - if self.compact: - dtype = np.dtype( - [ - ("timestepindex", np.int32), - ("particleid", np.int32), - ("node", np.int32), - ("x", np.float32), - ("y", np.float32), - ("z", np.float32), - ("zloc", np.float32), - ("time", np.float32), - ("timestep", np.int32), - ] - ) - else: - dtype = np.dtype( - [ - ("timestepindex", np.int32), - ("particleid", np.int32), - ("j", np.int32), - ("i", np.int32), - ("k", np.int32), - ("x", np.float32), - ("y", np.float32), - ("z", np.float32), - ("zloc", np.float32), - ("time", np.float32), - ("timestep", np.int32), - ] - ) - elif self.version == 6: + def _load(self) -> Tuple[np.dtype, np.ndarray]: + dtype = self.dtypes[self.version] + if not self.compact: dtype = np.dtype( [ - ("timepointindex", np.int32), - ("timestep", np.int32), - ("time", np.float32), + ("timestepindex", np.int32), ("particleid", np.int32), - ("particlegroup", np.int32), + ("j", np.int32), + ("i", np.int32), + ("k", np.int32), ("x", np.float32), ("y", np.float32), ("z", np.float32), - ("grid", np.int32), - ("k", np.int32), - ("i", np.int32), - ("j", np.int32), - ("xloc", np.float32), - ("yloc", np.float32), ("zloc", np.float32), - ] - ) - elif self.version == 7: - dtype = np.dtype( - [ - ("timepointindex", np.int32), - ("timestep", np.int32), ("time", np.float32), - ("particleid", np.int32), - ("particlegroup", np.int32), - ("particleidloc", np.int32), - ("node", np.int32), - ("xloc", np.float32), - ("yloc", np.float32), - ("zloc", np.float32), - ("x", np.float32), - ("y", np.float32), - ("z", np.float32), - ("k", np.int32), + ("timestep", np.int32), ] ) - return dtype - - def _get_outdtype(self): - outdtype = np.dtype( - [ - ("x", np.float32), - ("y", np.float32), - ("z", np.float32), - ("time", np.float32), - ("k", np.int32), - ("particleid", np.int32), - ] - ) - return outdtype - - def get_maxid(self): - """ - Get the maximum timeseries number in the file timeseries file - - Returns - ---------- - out : int - Maximum pathline number. - - """ - return super().get_maxid() - - def get_maxtime(self): - """ - Get the maximum time in timeseries file - - Returns - ---------- - out : float - Maximum pathline time. - - """ - return super().get_maxtime() - - def get_data(self, partid=0, totim=None, ge=True): - """ - get timeseries data from the timeseries file for a single timeseries - particleid. - Parameters - ---------- - partid : int - The zero-based particle id. The first record is record 0. - totim : float - The simulation time. All timeseries points for particle partid - that are greater than or equal to (ge=True) or less than or - equal to (ge=False) totim will be returned. Default is None - ge : bool - Boolean that determines if timeseries times greater than or equal - to or less than or equal to totim is used to create a subset - of timeseries. Default is True. - - Returns - ---------- - ra : numpy record array - A numpy recarray with the x, y, z, time, k, and particleid for - timeseries partid. - - - See Also - -------- - - Notes - ----- - - Examples - -------- - - >>> import flopy - >>> tsobj = flopy.utils.TimeseriesFile('model.timeseries') - >>> ts1 = tsobj.get_data(partid=1) - - """ - return super().get_data(partid=partid, totim=totim, ge=ge) - - def get_alldata(self, totim=None, ge=True): - """ - get timeseries data from the timeseries file for all timeseries - and all times. - - Parameters - ---------- - totim : float - The simulation time. All timeseries points for timeseries partid - that are greater than or equal to (ge=True) or less than or - equal to (ge=False) totim will be returned. Default is None - ge : bool - Boolean that determines if timeseries times greater than or equal - to or less than or equal to totim is used to create a subset - of timeseries. Default is True. + data = loadtxt(self.fname, dtype=dtype, skiprows=self.skiprows) - Returns - ---------- - tlist : a list of numpy record array - A list of numpy recarrays with the x, y, z, time, k, and - particleid for all timeseries. - - Examples - -------- + # convert indices to zero-based + for n in self.kijnames: + if n in data.dtype.names: + data[n] -= 1 - >>> import flopy - >>> tsobj = flopy.utils.TimeseriesFile('model.timeseries') - >>> ts = tsobj.get_alldata() + # sort by particle ID and time + data.sort(order=["particleid", "time"]) - """ - return super().get_alldata(totim=totim, ge=ge) + return dtype, data def get_destination_timeseries_data(self, dest_cells): """ @@ -1575,7 +847,7 @@ def get_destination_timeseries_data(self, dest_cells): Returns ------- - tsdest : np.recarray + data : np.recarray Slice of timeseries data array (e.g. TmeseriesFile._data) containing only timeseries that pass through (k,i,j) or (node) dest_cells. @@ -1590,59 +862,3 @@ def get_destination_timeseries_data(self, dest_cells): """ return super().get_destination_data(dest_cells=dest_cells) - - def write_shapefile( - self, - timeseries_data=None, - one_per_particle=True, - direction="ending", - shpname="pathlines.shp", - mg=None, - crs=None, - **kwargs, - ): - """ - Write pathlines to a shapefile - - timeseries_data : np.recarray - Record array of same form as that returned by - Timeseries.get_alldata(). (if none, Timeseries.get_alldata() - is exported). - one_per_particle : boolean (default True) - True writes a single LineString with a single set of attribute - data for each particle. False writes a record/geometry for each - pathline segment (each row in the Timeseries file). This option can - be used to visualize attribute information (time, model layer, - etc.) across a pathline in a GIS. - direction : str - String defining if starting or ending particle locations should be - included in shapefile attribute information. Only used if - one_per_particle=False. (default is 'ending') - shpname : str - File path for shapefile - mg : flopy.discretization.grid instance - Used to scale and rotate Global x,y,z values in MODPATH Timeseries - file. - crs : pyproj.CRS, int, str, optional - Coordinate reference system (CRS) for the model grid - (must be projected; geographic CRS are not supported). - The value can be anything accepted by - :meth:`pyproj.CRS.from_user_input() `, - such as an authority string (eg "EPSG:26916") or a WKT string. - kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp - - .. deprecated:: 3.5 - The following keyword options will be removed for FloPy 3.6: - - - ``epsg`` (int): use ``crs`` instead. - - """ - super().write_shapefile( - data=timeseries_data, - one_per_particle=one_per_particle, - direction=direction, - shpname=shpname, - mg=mg, - crs=crs, - **kwargs, - ) diff --git a/flopy/utils/particletracking.py b/flopy/utils/particletracking.py new file mode 100644 index 0000000000..8804471d25 --- /dev/null +++ b/flopy/utils/particletracking.py @@ -0,0 +1,333 @@ +""" +Utilities for parsing particle tracking output files. +""" + +from typing import Union +import numpy as np +from numpy.lib.recfunctions import stack_arrays + + +import os +from abc import ABC, abstractmethod +from pathlib import Path + + +class ParticleTrackFile(ABC): + """ + Abstract base class for particle track output files. Exposes a unified API + supporting MODPATH versions 3, 5, 6 and 7, as well as MODFLOW 6 PRT models. + + Notes + ----- + + + Parameters + ---------- + filename : str or PathLike + Path of output file + verbose : bool + Show verbose output. Default is False. + + """ + + outdtype = np.dtype( + [ + ("x", np.float32), + ("y", np.float32), + ("z", np.float32), + ("time", np.float32), + ("k", np.int32), + ("particleid", np.int32), + ] + ) + """ + Minimal information shared by all particle track file formats. + Track data are converted to this dtype for internal storage + and for return from (sub-)class methods. + """ + + def __init__( + self, + filename: Union[str, os.PathLike], + verbose: bool = False, + ): + self.fname = Path(filename).expanduser().absolute() + self.verbose = verbose + + def get_maxid(self) -> int: + """ + Get the maximum particle ID. + + Returns + ---------- + out : int + Maximum particle ID. + + """ + return self._data["particleid"].max() + + def get_maxtime(self) -> float: + """ + Get the maximum tracking time. + + Returns + ---------- + out : float + Maximum tracking time. + + """ + return self._data["time"].max() + + def get_data( + self, partid=0, totim=None, ge=True, minimal=False + ) -> np.recarray: + """ + Get a single particle track, optionally filtering by time. + + Parameters + ---------- + partid : int + Zero-based particle id. Default is 0. + totim : float + The simulation time. Default is None. + ge : bool + Filter tracking times greater than or equal + to or less than or equal to totim. + Only used if totim is not None. + minimal : bool + Whether to return only the minimal, canonical fields. Default is False. + + Returns + ---------- + data : np.recarray + Recarray with dtype ParticleTrackFile.outdtype + + """ + data = self._data[list(self.outdtype.names)] if minimal else self._data + idx = ( + np.where(data["particleid"] == partid)[0] + if totim is None + else ( + np.where( + (data["time"] >= totim) & (data["particleid"] == partid) + )[0] + if ge + else np.where( + (data["time"] <= totim) & (data["particleid"] == partid) + )[0] + ) + ) + + return data[idx] + + def get_alldata(self, totim=None, ge=True, minimal=False): + """ + Get all particle tracks separately, optionally filtering by time. + + Parameters + ---------- + totim : float + The simulation time. + ge : bool + Boolean that determines if pathline times greater than or equal + to or less than or equal to totim. + minimal : bool + Whether to return only the minimal, canonical fields. Default is False. + + Returns + ---------- + data : list of numpy record arrays + List of recarrays with dtype ParticleTrackFile.outdtype + + """ + nids = np.unique(self._data["particleid"]).size + data = self._data[list(self.outdtype.names)] if minimal else self._data + if totim is not None: + idx = ( + np.where(data["time"] >= totim)[0] + if ge + else np.where(data["time"] <= totim)[0] + ) + if len(idx) > 0: + data = data[idx] + return [data[data["particleid"] == i] for i in range(nids)] + + def get_destination_data( + self, dest_cells, to_recarray=True + ) -> np.recarray: + """ + Get data for set of destination cells. + + Parameters + ---------- + dest_cells : list or array of tuples + (k, i, j) of each destination cell for MODPATH versions less than + MODPATH 7 or node number of each destination cell. (zero based) + to_recarray : bool + Boolean that controls returned series. If to_recarray is True, + a single recarray with all of the pathlines that intersect + dest_cells are returned. If to_recarray is False, a list of + recarrays (the same form as returned by get_alldata method) + that intersect dest_cells are returned (default is False). + + Returns + ------- + data : np.recarray + Slice of data array (e.g. PathlineFile._data, TimeseriesFile._data) + containing endpoint, pathline, or timeseries data that intersect + (k,i,j) or (node) dest_cells. + + """ + + return self.intersect(dest_cells, to_recarray) + + @abstractmethod + def intersect(self, cells, to_recarray) -> np.recarray: + """Find intersection of pathlines with cells.""" + pass + + def write_shapefile( + self, + data=None, + one_per_particle=True, + direction="ending", + shpname="endpoints.shp", + mg=None, + crs=None, + **kwargs, + ): + """ + Write particle track data to a shapefile. + + data : np.recarray + Record array of same form as that returned by + get_alldata(). (if none, get_alldata() is exported). + one_per_particle : boolean (default True) + True writes a single LineString with a single set of attribute + data for each particle. False writes a record/geometry for each + pathline segment (each row in the PathLine file). This option can + be used to visualize attribute information (time, model layer, + etc.) across a pathline in a GIS. + direction : str + String defining if starting or ending particle locations should be + included in shapefile attribute information. Only used if + one_per_particle=False. (default is 'ending') + shpname : str + File path for shapefile + mg : flopy.discretization.grid instance + Used to scale and rotate Global x,y,z values. + crs : pyproj.CRS, int, str, optional + Coordinate reference system (CRS) for the model grid + (must be projected; geographic CRS are not supported). + The value can be anything accepted by + :meth:`pyproj.CRS.from_user_input() `, + such as an authority string (eg "EPSG:26916") or a WKT string. + kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp + + .. deprecated:: 3.5 + The following keyword options will be removed for FloPy 3.6: + + - ``epsg`` (int): use ``crs`` instead. + + """ + from ..export.shapefile_utils import recarray2shp + from ..utils import geometry + from ..utils.geometry import LineString + + series = data + if series is None: + series = self._data.view(np.recarray) + else: + # convert pathline list to a single recarray + if isinstance(series, list): + s = series[0] + print(s.dtype) + for n in range(1, len(series)): + s = stack_arrays((s, series[n])) + series = s.view(np.recarray) + + series = series.copy() + series.sort(order=["particleid", "time"]) + + if mg is None: + raise ValueError("A modelgrid object was not provided.") + + particles = np.unique(series.particleid) + geoms = [] + + # create dtype with select attributes in pth + names = series.dtype.names + dtype = [] + atts = ["particleid", "particlegroup", "time", "k", "i", "j", "node"] + for att in atts: + if att in names: + t = np.int32 + if att == "time": + t = np.float32 + dtype.append((att, t)) + dtype = np.dtype(dtype) + + # reset names to the selected names in the created dtype + names = dtype.names + + # 1 geometry for each path + if one_per_particle: + loc_inds = 0 + if direction == "ending": + loc_inds = -1 + + sdata = [] + for pid in particles: + ra = series[series.particleid == pid] + + x, y = geometry.transform( + ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians + ) + z = ra.z + geoms.append(LineString(list(zip(x, y, z)))) + + t = [pid] + if "particlegroup" in names: + t.append(ra.particlegroup[0]) + t.append(ra.time.max()) + if "node" in names: + t.append(ra.node[loc_inds]) + else: + if "k" in names: + t.append(ra.k[loc_inds]) + if "i" in names: + t.append(ra.i[loc_inds]) + if "j" in names: + t.append(ra.j[loc_inds]) + sdata.append(tuple(t)) + + sdata = np.array(sdata, dtype=dtype).view(np.recarray) + + # geometry for each row in PathLine file + else: + dtype = series.dtype + sdata = [] + for pid in particles: + ra = series[series.particleid == pid] + if mg is not None: + x, y = geometry.transform( + ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians + ) + else: + x, y = geometry.transform(ra.x, ra.y, 0, 0, 0) + z = ra.z + geoms += [ + LineString( + [(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])] + ) + for i in np.arange(1, (len(ra))) + ] + sdata += ra[1:].tolist() + sdata = np.array(sdata, dtype=dtype).view(np.recarray) + + # convert back to one-based + for n in set(self.kijnames).intersection(set(sdata.dtype.names)): + sdata[n] += 1 + + # write the final recarray to a shapefile + recarray2shp(sdata, geoms, shpname=shpname, crs=crs, **kwargs)