@@ -2,6 +2,7 @@
Loading
2 2
import sys
3 3
4 4
import numpy as np
5 +
import matplotlib.pyplot as plt
5 6
import matplotlib.cm as cm
6 7
import matplotlib.colors as colors
7 8
from scipy import optimize
@@ -71,22 +72,6 @@
Loading
71 72
                      "or stratigraphy information."
72 73
        super().__init__(message)
73 74
74 -
"""
75 -
    yields an exception with:
76 -
77 -
    .. code::
78 -
79 -
        deltametrics.utils.NoStratigraphyError: 'DataCube' object
80 -
        has no preservation or stratigraphy information.
81 -
82 -
83 -
    .. code::
84 -
85 -
        deltametrics.utils.NoStratigraphyError: 'DataCube' object
86 -
        has no attribute 'strat_attr'.
87 -
88 -
89 -
"""
90 75
91 76
def needs_stratigraphy(func):
92 77
    """Decorator for properties requiring stratigraphy.
@@ -153,42 +138,40 @@
Loading
153 138
154 139
155 140
def curve_fit(data, fit='harmonic'):
156 -
    """
157 -
    Calculate curve fit given some data.
141 +
    """Calculate curve fit given some data.
158 142
159 143
    Several functional forms are available for fitting: exponential, harmonic,
160 144
    and linear. The input `data` can be 1-D, or 2-D, if it is 2-D, the data
161 145
    will be averaged. The expected 2-D shape is (Y-Values, # Values) where the
162 146
    data you wish to have fit is in the first dimension, and the second
163 -
    dimension is of len(# Values).
147 +
    dimension is of ``len(# Values)``.
164 148
165 149
    E.g. Given some mobility data output from one of the mobility metrics,
166 150
    fit a curve to the average of that data.
167 151
168 152
    Parameters
169 153
    ----------
170 -
    data : ndarray
154 +
    data : :obj:`ndarray`
171 155
        Data, either already averaged or a 2D array of of shape
172 156
        len(data values) x len(# values).
173 157
174 -
    fit : str, optional (default is 'harmonic')
158 +
    fit : :obj:`str`, optional
175 159
        A string specifying the type of function to be fit. Options are as
176 160
        follows:
177 -
            - 'exponential' : (a - b) * np.exp(-c * x) + b
178 -
            - 'harmonic' : a / (1 + b * x)
179 -
            - 'linear' : a * x + b
161 +
          * `exponential`, which evaluates :code:`(a - b) * np.exp(-c * x) + b`
162 +
          * `harmonic`, (default) which evaluates :code:`a / (1 + b * x)`
163 +
          * `linear`, which evaluates :code:`a * x + b`
180 164
181 165
    Returns
182 166
    -------
183 -
    yfit : ndarray
167 +
    yfit : :obj:`ndarray`
184 168
        y-values corresponding to the fitted function.
185 169
186 -
    pcov : ndarray
170 +
    pcov : :obj:`ndarray`
187 171
        Covariance associated with the fitted function parameters.
188 172
189 -
    perror : ndarray
190 -
        One standard deviation error for the parameters (from pcov)
191 -
173 +
    perror : :obj:`ndarray`
174 +
        One standard deviation error for the parameters (from pcov).
192 175
    """
193 176
    avail_fits = ['exponential', 'harmonic', 'linear']
194 177
    if fit not in avail_fits:
@@ -218,3 +201,280 @@
Loading
218 201
    perror = np.sqrt(np.diag(pcov))
219 202
220 203
    return yfit, pcov, perror
204 +
205 +
206 +
def guess_land_width_from_land(land_col_0):
207 +
    """Guess the land width from bed elevations.
208 +
209 +
    Utility to help autodetermine the domain setup. This utility should be
210 +
    replaced when possible by pulling domain setup variables directly from the
211 +
    netCDF file.
212 +
213 +
    Algortihm works by finding the point where the bed elevation is *flat*
214 +
    (i.e., where there is undisturbed basin).
215 +
    """
216 +
    i = 0
217 +
    delt = 10
218 +
    while i < len(land_col_0) - 1 and delt != 0:
219 +
        delt = land_col_0[i] - land_col_0[i+1]
220 +
        i += 1
221 +
    trim_idx = i - 1  # assign the trimming index
222 +
    return trim_idx
223 +
224 +
225 +
def coordinates_to_segments(coordinates):
226 +
    """Transform coordinate [x, y] array into line segments.
227 +
228 +
    Parameters
229 +
    ----------
230 +
    coordinates : :obj:`ndarray`
231 +
        `[N, 2]` array with `(x, y)` coordinate pairs in rows.
232 +
233 +
    Returns
234 +
    -------
235 +
    segments
236 +
        `[(N-1), 2, 2` array of line segments with dimensions of
237 +
        `each segment x x-coordinates x y-coordinates`.
238 +
    """
239 +
    _x = np.vstack((coordinates[:-1, 0],
240 +
                    coordinates[1:, 0])).T.reshape(-1, 2, 1)
241 +
    _y = np.vstack((coordinates[:-1, 1],
242 +
                    coordinates[1:, 1])).T.reshape(-1, 2, 1)
243 +
    return np.concatenate([_x, _y], axis=2)
244 +
245 +
246 +
def segments_to_cells(segments):
247 +
    """Transform a line segment (or array of segments) into integer coords.
248 +
249 +
    Helper function to convert a path into cells. This is generally used for
250 +
    determining the path of a `Section`, but is also available to use on any
251 +
    regularly gridded array by metric computation functions.
252 +
253 +
    Type and shape checking is performed, then the path (which may be
254 +
    composed of multiple vertices) is converted to a single path of cells.
255 +
    """
256 +
    _c = []  # append lists, do not know length of cells a priori
257 +
    for s in np.arange(segments.shape[0]):
258 +
        _c.append(line_to_cells(segments[s, ...]))
259 +
    _c = np.hstack(_c).T
260 +
    return _c
261 +
262 +
263 +
def line_to_cells(*args):
264 +
    """Convert a line to cell coordinates along the line.
265 +
266 +
    Takes as input the line to determine coordinates along. A line defined by
267 +
    two Cartesian coordinate endpoints `p1` and `p2`, may be specified  in one
268 +
    of three ways:
269 +
270 +
    * a single `ndarray` with fields ``[[x0, y0], [x1, y1]]``
271 +
    * a set of two-tuples with fields ``(x0, y0), (x1, y1)``
272 +
    * a set of four floats ``x0, y0, x1, y1``
273 +
274 +
    Cell coordinates will be ordered according to the order of the arguments
275 +
    (i.e., as though walking along the line from `p0` to `p1`).
276 +
277 +
    Parameters
278 +
    ----------
279 +
    line : four :obj:`float`, two :obj:`tuple`, one :obj:`ndarray`
280 +
        The line, defined in one of the options described above.
281 +
282 +
    Returns
283 +
    -------
284 +
    x : :obj:`ndarray`
285 +
        x-coordinates of cells intersected with line.
286 +
287 +
    y : :obj:`ndarray`
288 +
        y-coordinates of cells intersected with line.
289 +
290 +
    Examples
291 +
    --------
292 +
293 +
    .. plot::
294 +
        :include-source:
295 +
296 +
        >>> from deltametrics.utils import line_to_cells
297 +
        >>> p0 = (1, 6)
298 +
        >>> p1 = (6, 3)
299 +
        >>> x, y = line_to_cells(p0, p1)
300 +
301 +
        >>> fig, ax = plt.subplots(figsize=(2, 2))
302 +
        >>> _arr = np.zeros((10, 10))
303 +
        >>> _arr[y, x] = 1
304 +
        >>> ax.imshow(_arr, cmap='gray')
305 +
        >>> ax.plot([p0[0], p1[0]], [p0[1], p1[1]], 'r-')
306 +
        >>> plt.show()
307 +
    """
308 +
    # process the args
309 +
    if len(args) == 1:
310 +
        x0, y0, x1, y1 = args[0].ravel()
311 +
    elif len(args) == 2:
312 +
        x0, y0 = args[0]
313 +
        x1, y1 = args[1]
314 +
    elif len(args) == 4:
315 +
        x0, y0, x1, y1 = args
316 +
    else:
317 +
        raise TypeError
318 +
319 +
    # process the line to cells
320 +
    if np.abs(y1 - y0) < np.abs(x1 - x0):
321 +
        # if the line is "shallow" (dy < dx)
322 +
        if x0 > x1:
323 +
            # if the line is trending down (III)
324 +
            x, y = _walk_line((x1, y1), (x0, y0))
325 +
        else:
326 +
            # if the line is trending up (I)
327 +
            x, y = _walk_line((x0, y0), (x1, y1))
328 +
    else:
329 +
        # if the line is "steep" (dy >= dx)
330 +
        if y0 > y1:
331 +
            # if the line is trending down (IV)
332 +
            y, x = _walk_line((y1, x1), (y0, x0))
333 +
        else:
334 +
            # if the line is trending up (II)
335 +
            y, x = _walk_line((y0, x0), (y1, x1))
336 +
337 +
    return x, y
338 +
339 +
340 +
def _walk_line(p0, p1):
341 +
    """Walk a line to determine cells along path.
342 +
343 +
    Inputs depend on the steepness and direction of the input line. For a
344 +
    shallow line, where dx > dy, the input tuples should be in the form of
345 +
    `(x, y)` pairs. In contrast, for steep lines (where `dx < dy`), the input
346 +
    tuples should be `(y, x)` pairs.
347 +
348 +
    Additionally, the order of the tuples should be given based on whether the
349 +
    line is trending upward or downward, with increasing `x`.
350 +
351 +
    .. note:: TODO: finish descriptions for quadrants
352 +
353 +
    .. note::
354 +
355 +
        `x` in the code may not be cartesian x. Depends on quadrant of the
356 +
        line.
357 +
    """
358 +
    # unpack the point tuples
359 +
    x0, y0 = p0
360 +
    x1, y1 = p1
361 +
362 +
    dx, dy = x1 - x0, y1 - y0
363 +
    yi = 1
364 +
    if dy < 0:
365 +
        yi = -1
366 +
        dy = -dy
367 +
368 +
    D = 2 * dy - dx
369 +
    x = np.arange(x0, x1 + 1, dtype=np.int).T
370 +
    y = np.zeros((len(x),), dtype=np.int)
371 +
372 +
    yy = y0
373 +
    for i in np.arange(len(x)):
374 +
        y[i] = yy
375 +
        if D > 0:
376 +
            yy = yy + yi
377 +
            D = D - 2 * dx
378 +
379 +
        D = D + 2 * dy
380 +
381 +
    # sort by major axis, and index the cells
382 +
    xI = np.argsort(x)
383 +
    x = x[xI]
384 +
    y = y[xI]
385 +
386 +
    return x, y
387 +
388 +
389 +
def circle_to_cells(origin, radius, remove_duplicates=True):
390 +
    """The x, y pixel coordinates of a circle
391 +
392 +
    Use the mid-point circle algorithm is used for computation
393 +
    (http://en.wikipedia.org/wiki/Midpoint_circle_algorithm).
394 +
395 +
    Compute in advance the number of points that will be generated by the
396 +
    algorithm, to pre-allocate the coordinates arrays. Also, this function
397 +
    ensures that sorted coordinates are returned.
398 +
399 +
    This implementation removes duplicate points. Optionally, specify
400 +
    `remove_duplicates=False` to keep duplicates and achieve faster execution.
401 +
    Note that keeping duplicates is fine for drawing a circle with colored
402 +
    pixels, but for slicing an array along the arc, we need to have only
403 +
    unique and sorted indices.
404 +
405 +
    Original implementation from Jean-Yves Tinevez.
406 +
    <jeanyves.tinevez@gmail.com> - Nov 2011 - Feb 2012
407 +
    """
408 +
    x0, y0 = origin
409 +
410 +
    # Compute first the number of points
411 +
    octant_size = np.int((np.sqrt(2) * (radius - 1) + 4) / 2)
412 +
    n_points = 4 * octant_size
413 +
    xc = np.zeros((n_points,), dtype=np.int)
414 +
    yc = np.zeros((n_points,), dtype=np.int)
415 +
416 +
    x = 0
417 +
    y = radius
418 +
    f = 1 - radius
419 +
    dx = 1
420 +
    dy = - 2 * radius
421 +
422 +
    # 7th octant -- driver
423 +
    xc[0 * octant_size] = x0 - y
424 +
    yc[0 * octant_size] = y0 + x
425 +
    # 8th octant
426 +
    xc[2 * octant_size - 1] = x0 - x
427 +
    yc[2 * octant_size - 1] = y0 + y
428 +
    # 1st octant
429 +
    xc[2 * octant_size] = x0 + x
430 +
    yc[2 * octant_size] = y0 + y
431 +
    # 2nd octant
432 +
    xc[4 * octant_size - 1] = x0 + y
433 +
    yc[4 * octant_size - 1] = y0 + x
434 +
435 +
    for i in np.arange(1, n_points / 4, dtype=np.int):
436 +
        # update x and y, follwing midpoint algo
437 +
        if f > 0:
438 +
            y = y - 1
439 +
            dy = dy + 2
440 +
            f = f + dy
441 +
        x = x + 1
442 +
        dx = dx + 2
443 +
        f = f + dx
444 +
445 +
        # 7th octant
446 +
        xc[i] = x0 - y
447 +
        yc[i] = y0 + x
448 +
        # 8th octant
449 +
        xc[2 * octant_size - i - 1] = x0 - x
450 +
        yc[2 * octant_size - i - 1] = y0 + y
451 +
        # 1st octant
452 +
        xc[2 * octant_size + i] = x0 + x
453 +
        yc[2 * octant_size + i] = y0 + y
454 +
        # 2nd octant
455 +
        xc[4 * octant_size - i - 1] = x0 + y
456 +
        yc[4 * octant_size - i - 1] = y0 + x
457 +
458 +
    # There may be some duplicate entries
459 +
    #     We loop through to remove duplicates. This is slow, but necessary in
460 +
    #     most of our applications. We have to use something custom, rather
461 +
    #     than np.unique() because we need to preserve the ordering of the
462 +
    #     octants.
463 +
    if remove_duplicates:
464 +
        xyc = np.column_stack((xc, yc))
465 +
        keep = np.ones((n_points,), dtype=np.bool)
466 +
        for i in np.arange(1, 4):
467 +
            prv = xyc[(i-1)*octant_size:i*octant_size, :]
468 +
            nxt = xyc[i*octant_size:(i+1)*octant_size, :]
469 +
            dupe = np.nonzero(np.all(prv == nxt[:, np.newaxis], axis=2))[0]
470 +
            keep[(i*octant_size)+dupe] = False
471 +
        xyc = xyc[keep]
472 +
        xc = xyc[:, 0]
473 +
        yc = xyc[:, 1]
474 +
475 +
    # limit to positive indices (no wrapping)
476 +
    _and = np.logical_and(xc >= 0, yc >= 0)
477 +
    xc = xc[_and]
478 +
    yc = yc[_and]
479 +
480 +
    return xc, yc

@@ -1,4 +1,5 @@
Loading
1 1
import abc
2 +
import warnings
2 3
3 4
import numpy as np
4 5
from scipy import stats, sparse
@@ -178,34 +179,39 @@
Loading
178 179
179 180
    """
180 181
181 -
    def __init__(self, section_type, *args):
182 +
    def __init__(self, section_type, *args, name=None):
182 183
        """
183 184
        Identify coordinates defining the section.
184 185
185 186
        Parameters
186 187
        ----------
187 188
        CubeInstance : :obj:`~deltametrics.cube.Cube` subclass instance, optional
188 -
            Connect to this cube. No connection is made if cube is not provided.
189 +
            Connect to this cube. No connection is made if cube is not
190 +
            provided.
189 191
190 192
        Notes
191 193
        -----
192 194
193 195
        If no arguments are passed, an empty section not connected to any cube
194 196
        is returned. This cube will will need to be manually connected to have
195 -
        any functionality (via the :meth:`connect` method.
197 +
        any functionality (via the :meth:`connect` method).
196 198
        """
197 199
        # begin unconnected
198 200
        self._s = None
199 201
        self._z = None
200 202
        self._x = None
201 203
        self._y = None
204 +
        self._trace = None
205 +
        self._shape = None
202 206
        self._variables = None
203 207
        self.cube = None
204 208
205 209
        self.section_type = section_type
210 +
        self._name = name
206 211
207 212
        if len(args) > 1:
208 -
            raise ValueError('Expected single argument to %s instantiation.'
213 +
            raise ValueError('Expected single positional argument to \
214 +
                             %s instantiation.'
209 215
                             % type(self))
210 216
211 217
        if len(args) > 0:
@@ -213,7 +219,7 @@
Loading
213 219
        else:
214 220
            pass
215 221
216 -
    def connect(self, CubeInstance):
222 +
    def connect(self, CubeInstance, name=None):
217 223
        """Connect this Section instance to a Cube instance.
218 224
        """
219 225
        if not issubclass(type(CubeInstance), cube.BaseCube):
@@ -223,9 +229,29 @@
Loading
223 229
                                _gottype=type(CubeInstance)))
224 230
        self.cube = CubeInstance
225 231
        self._variables = self.cube.variables
232 +
        # self._name = self._name or name or self.section_type
233 +
        self.name = name  # use the setter to determine the _name
226 234
        self._compute_section_coords()
227 235
        self._compute_section_attrs()
228 236
237 +
    @property
238 +
    def name(self):
239 +
        return self._name
240 +
241 +
    @name.setter
242 +
    def name(self, var):
243 +
        if (self._name is None):
244 +
            # _name is not yet set
245 +
            self._name = var or self.section_type
246 +
        else:
247 +
            # _name is already set
248 +
            if not (var is None):
249 +
                warnings.warn(UserWarning("`name` argument supplied to \
250 +
                    instantiated `Section` object. To change the name of \
251 +
                    a Section, you must set the attribute directly with \
252 +
                    `section._name = 'name'`."))
253 +
            # do nothing
254 +
229 255
    @abc.abstractmethod
230 256
    def _compute_section_coords(self):
231 257
        """Should calculate x-y coordinates of the section.
@@ -250,12 +276,14 @@
Loading
250 276
        self._s = np.cumsum(np.hstack((0, np.sqrt((self._x[1:] - self._x[:-1])**2
251 277
                                                  + (self._y[1:] - self._y[:-1])**2))))
252 278
        self._z = self.cube.z
279 +
        self._shape = (len(self._z), len(self._s))
280 +
        self._trace = np.column_stack((self._x, self._y))
253 281
254 282
    @property
255 283
    def trace(self):
256 284
        """Coordinates of the section in the x-y plane.
257 285
        """
258 -
        return np.column_stack((self._x, self._y))
286 +
        return self._trace
259 287
260 288
    @property
261 289
    def s(self):
@@ -267,6 +295,14 @@
Loading
267 295
        """Up-section (vertical) coordinate."""
268 296
        return self._z
269 297
298 +
    @property
299 +
    def shape(self):
300 +
        """Section shape.
301 +
302 +
        Simply a `tuple` equivalent to ``(len(z), len(s))``
303 +
        """
304 +
        return self._shape
305 +
270 306
    @property
271 307
    def variables(self):
272 308
        """List of variables.
@@ -275,6 +311,13 @@
Loading
275 311
276 312
    @property
277 313
    def strat_attr(self):
314 +
        """Stratigraphic attributes data object.
315 +
316 +
        Raises
317 +
        ------
318 +
        NoStratigraphyError
319 +
            If no stratigraphy information is found for the section.
320 +
        """
278 321
        if self.cube._knows_stratigraphy:
279 322
            return self.cube.strat_attr
280 323
        else:
@@ -312,12 +355,12 @@
Loading
312 355
                                            'section', self._y, self._x))
313 356
            else:
314 357
                return DataSectionVariable(_data=self.cube[var].data.values[:,
315 -
                                                                     self.y,
358 +
                                                                     self._y,
316 359
                                                                     self._x],
317 360
                                           _s=self.s, _z=self.z)
318 361
        elif type(self.cube) is cube.StratigraphyCube:
319 362
            return StratigraphySectionVariable(_data=self.cube[var].data.values[:,
320 -
                                                                    self.y,
363 +
                                                                    self._y,
321 364
                                                                    self._x],
322 365
                                               _s=self.s, _z=self.z)
323 366
        elif self.cube is None:
@@ -365,6 +408,11 @@
Loading
365 408
            :obj:`~deltametrics.plot.VariableSet` is used. Other arguments are
366 409
            attempted to coerce to `str`, and the literal is diplayed.
367 410
411 +
        ax : :obj:`~matplotlib.pyplot.Axes` object, optional
412 +
            A `matplotlib` `Axes` object to plot the section. Optional; if not
413 +
            provided, a call is made to ``plt.gca()`` to get the current (or
414 +
            create a new) `Axes` object.
415 +
368 416
        Examples
369 417
        --------
370 418
        *Example 1:* Display the `velocity` spacetime section of a DataCube.
@@ -441,43 +489,116 @@
Loading
441 489
        ax.set_xlim(xmin, xmax)
442 490
        ax.set_ylim(ymin, ymax)
443 491
492 +
    def show_trace(self, *args, ax=None, **kwargs):
493 +
        """Plot section trace (x-y plane path).
444 494
445 -
class PathSection(BaseSection):
446 -
    """Path section object.
495 +
        Plot the section trace (:obj:`trace`) onto an x-y planview.
447 496
448 -
    Create a Section along user-specified path.
497 +
        Parameters
498 +
        ----------
499 +
        *args
500 +
            Passed to `matplotlib` :obj:`~matplotlib.pyplot.plot()`.
449 501
450 -
    .. note::
502 +
        ax : :obj:`~matplotlib.pyplot.Axes` object, optional
503 +
            A `matplotlib` `Axes` object to plot the trace. Optional; if not
504 +
            provided, a call is made to ``plt.gca()`` to get the current (or
505 +
            create a new) `Axes` object.
506 +
507 +
        **kwargs
508 +
            Passed to `matplotlib` :obj:`~matplotlib.pyplot.plot()`.
509 +
        """
510 +
        if not ax:
511 +
            ax = plt.gca()
512 +
513 +
        _label = kwargs.pop('label', self.name)
514 +
515 +
        ax.plot(self._x, self._y, label=_label, *args, **kwargs)
516 +
517 +
518 +
class PathSection(BaseSection):
519 +
    """Path section object.
451 520
452 -
        Currently, this extracts only *at* the points specified. A good
453 -
        improvement would be to interpolate along the path defined, and
454 -
        extract the section everywhere the path intersects within 50% of the
455 -
        center of the surface area of a grid cell.
521 +
    Create a Section along user-specified path. Specify the section location
522 +
    as an `(N, 2)` `ndarray` of x-y pairs of coordinates that define the
523 +
    verticies of the path. All coordinates along the path will be included in
524 +
    the section.
525 +
526 +
    .. important::
527 +
528 +
        The vertex coordinates must be specified as cell indices (not
529 +
        actual x and y coordinate values). This is a needed patch.
530 +
531 +
    Parameters
532 +
    ----------
533 +
    *args : :obj:`DataCube` or `StratigraphyCube`
534 +
        The `Cube` object to link for underlying data. This option should be
535 +
        ommitted if using the :obj:`register_section` method of a `Cube`.
536 +
537 +
    path : :obj:`ndarray`
538 +
        An `(N, 2)` `ndarray` specifying the x-y pairs of coordinates that
539 +
            define the verticies of the path to extract the section from.
540 +
541 +
    **kwargs
542 +
        Keyword arguments are passed to `BaseSection.__init__()`. Supported
543 +
        options are `name`.
544 +
545 +
    Returns
546 +
    -------
547 +
    section : :obj:`PathSection`
548 +
        `PathSection` object with specified parameters. The section is
549 +
        automatically connected to the underlying `Cube` data source if the
550 +
        :obj:`register_section` method of a `Cube` is used to set up the
551 +
        section, or the `Cube` is passed as the first positional argument
552 +
        during instantiation.
553 +
554 +
    Examples
555 +
    --------
556 +
557 +
    To create a `PathSection` that is registered to a `DataCube` at
558 +
    specified coordinates:
559 +
560 +
    .. plot::
561 +
        :include-source:
562 +
563 +
        >>> rcm8cube = dm.sample_data.cube.rcm8()
564 +
        >>> rcm8cube.register_section('path', section.PathSection(
565 +
        ...     path=np.array([[50, 3], [65, 17], [130, 10]])))
566 +
        >>>
567 +
        >>> # show the location and the "velocity" variable
568 +
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
569 +
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
570 +
        >>> rcm8cube.sections['path'].show_trace('r--', ax=ax[0])
571 +
        >>> rcm8cube.sections['path'].show('velocity', ax=ax[1])
572 +
        >>> plt.show()
456 573
    """
457 574
458 -
    def __init__(self, *args, path):
575 +
    def __init__(self, *args, path, **kwargs):
459 576
        """Instantiate.
460 577
461 578
        Parameters
462 579
        ----------
463 580
        path : :obj:`ndarray`
464 -
            An Mx2 `ndarray` specifying the x-y pairs of coordinates to
465 -
            extract the section from.
581 +
            An `(N, 2)` `ndarray` specifying the x-y pairs of coordinates that
582 +
            define the verticies of the path to extract the section from.
466 583
467 -
        Notes
468 -
        -----
584 +
        .. note::
469 585
470 -
        `path` must be supplied as a keyword argument.
586 +
            :obj:`path` must be supplied as a keyword argument.
471 587
472 588
        """
473 589
        self._input_path = path
474 -
        super().__init__('path', *args)
590 +
        super().__init__('path', *args, **kwargs)
475 591
476 592
    def _compute_section_coords(self):
477 593
        """Calculate coordinates of the strike section.
478 594
        """
595 +
        # convert the points into segments into lists of cells
596 +
        _segs = utils.coordinates_to_segments(self._input_path)
597 +
        _cell = utils.segments_to_cells(_segs)
598 +
479 599
        # determine only unique coordinates along the path
480 -
        self._path = np.unique(self._input_path, axis=0)
600 +
        self._path = np.unique(_cell, axis=0)
601 +
        self._vertices = np.unique(self._input_path, axis=0)
481 602
482 603
        self._x = self._path[:, 0]
483 604
        self._y = self._path[:, 1]
@@ -493,18 +614,98 @@
Loading
493 614
494 615
class StrikeSection(BaseSection):
495 616
    """Strike section object.
617 +
618 +
    Section oriented along the delta strike (i.e., perpendicular to an inlet
619 +
    channel). Specify the location of the strike section with :obj`y` and
620 +
    :obj:`x` keyword parameter options.
621 +
622 +
    .. important::
623 +
624 +
        The `y` and `x` parameters must be specified as cell indices (not
625 +
        actual x and y coordinate values). This is a needed patch.
626 +
627 +
    Parameters
628 +
    ----------
629 +
    *args : :obj:`DataCube` or `StratigraphyCube`
630 +
        The `Cube` object to link for underlying data. This option should be
631 +
        ommitted if using the :obj:`register_section` method of a `Cube`.
632 +
633 +
    y : :obj:`int`, optional
634 +
        The `y` location of the section. This is the distance to locate the
635 +
        section from the domain edge with a channel inlet. Defaults to ``0``
636 +
        if no value is given.
637 +
638 +
    x : :obj:`int`, optional
639 +
        The `x` limits for the section. Defaults to the full domain width.
640 +
        Specify as a two-element `tuple` or `list` of `int`, giving the lower
641 +
        and upper bounds of `x` values to span the section.
642 +
643 +
    **kwargs
644 +
        Keyword arguments are passed to `BaseSection.__init__()`. Supported
645 +
        options are `name`.
646 +
647 +
    Returns
648 +
    -------
649 +
    section : :obj:`StrikeSection`
650 +
        `StrikeSection` object with specified parameters. The section is
651 +
        automatically connected to the underlying `Cube` data source if the
652 +
        :obj:`register_section` method of a `Cube` is used to set up the
653 +
        section, or the `Cube` is passed as the first positional argument
654 +
        during instantiation.
655 +
656 +
    Examples
657 +
    --------
658 +
659 +
    To create a `StrikeSection` that is registered to a `DataCube` at
660 +
    specified `y` coordinate ``=10``, and spans the entire model domain:
661 +
662 +
    .. plot::
663 +
        :include-source:
664 +
665 +
        >>> rcm8cube = dm.sample_data.cube.rcm8()
666 +
        >>> rcm8cube.register_section('strike', dm.section.StrikeSection(y=10))
667 +
        >>>
668 +
        >>> # show the location and the "velocity" variable
669 +
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
670 +
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
671 +
        >>> rcm8cube.sections['strike'].show_trace('r--', ax=ax[0])
672 +
        >>> rcm8cube.sections['strike'].show('velocity', ax=ax[1])
673 +
        >>> plt.show()
674 +
675 +
    Similarly, create a `StrikeSection` that is registered to a
676 +
    `StratigraphyCube` at specified `y` coordinate ``=20``, and spans only the
677 +
    left side of the model domain:
678 +
679 +
    .. plot::
680 +
        :include-source:
681 +
682 +
        >>> rcm8cube = dm.sample_data.cube.rcm8()
683 +
        >>> sc8cube = dm.cube.StratigraphyCube.from_DataCube(rcm8cube)
684 +
        >>> sc8cube.register_section('strike_half', dm.section.StrikeSection(y=20, x=[0, 120]))
685 +
        >>>
686 +
        >>> # show the location and the "velocity" variable
687 +
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
688 +
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
689 +
        >>> sc8cube.sections['strike_half'].show_trace('r--', ax=ax[0])
690 +
        >>> sc8cube.sections['strike_half'].show('velocity', ax=ax[1])
691 +
        >>> plt.show()
496 692
    """
497 693
498 -
    def __init__(self, *args, y=0):
694 +
    def __init__(self, *args, y=None, x=None, **kwargs):
499 695
500 696
        self.y = y  # strike coord scalar
501 -
        super().__init__('strike', *args)
697 +
        self._input_xlim = x  # the input x lims
698 +
        super().__init__('strike', *args, **kwargs)
502 699
503 700
    def _compute_section_coords(self):
504 701
        """Calculate coordinates of the strike section.
505 702
        """
506 -
        _nx = self.cube.shape[2]
507 -
        self._x = np.arange(_nx)
703 +
        if self._input_xlim is None:
704 +
            _nx = self.cube['eta'].shape[2]
705 +
            self._x = np.arange(_nx)
706 +
        else:
707 +
            self._x = np.arange(self._input_xlim[0], self._input_xlim[1])
708 +
            _nx = len(self._x)
508 709
        self._y = np.tile(self.y, (_nx))
509 710
510 711
@@ -518,10 +719,254 @@
Loading
518 719
        # choose center point if x=-1
519 720
520 721
722 +
class CircularSection(BaseSection):
723 +
    """Circular section object.
724 +
725 +
    Section drawn as a circular cut, located a along the arc a specified
726 +
    radius from specified origin.  Specify the location of the circular section
727 +
    with :obj`radius` and :obj:`origin` keyword parameter options. The
728 +
    circular section trace is interpolated to the nearest integer model domain
729 +
    cells, following the mid-point circle algorithm
730 +
    (:obj:`~deltametrics.utils.circle_to_cells`).
731 +
732 +
    .. important::
733 +
734 +
        The `radius` and `origin` parameters must be specified as cell indices
735 +
        (not actual x and y coordinate values). This is a needed patch.
736 +
737 +
    .. important::
738 +
739 +
        The `origin` attempts to detect the land width from bed elevation
740 +
        changes, but should use the value of ``L0`` recorded in the netcdf
741 +
        file, or defined in the cube.
742 +
743 +
    Parameters
744 +
    ----------
745 +
    *args : :obj:`DataCube` or `StratigraphyCube`
746 +
        The `Cube` object to link for underlying data. This option should be
747 +
        ommitted if using the :obj:`register_section` method of a `Cube`.
748 +
749 +
    radius : :obj:`float`, `int`, optional
750 +
        The `radius` of the section. This is the distance to locate the
751 +
        section from the :obj:`origin`. If no value is given, the `radius`
752 +
        defaults to half of the minimum model domain edge length if it can be
753 +
        determined, otherwise defaults to ``1``.
754 +
755 +
    origin : :obj:`tuple` or `list` of `int`, optional
756 +
        The `origin` of the circular section. This is the center of the
757 +
        circle. If no value is given, the origin defaults to the center of the
758 +
        x-direction of the model domain, and offsets into the domain a
759 +
        distance of ``y == L0``, if these values can be determined. I.e., the
760 +
        origin defaults to be centered over the channel inlet. If no value is
761 +
        given, and these values cannot be determined, the origin defaults to
762 +
        ``(0, 0)``.
763 +
764 +
    **kwargs
765 +
        Keyword arguments are passed to `BaseSection.__init__()`. Supported
766 +
        options are `name`.
767 +
768 +
    Returns
769 +
    -------
770 +
    section : :obj:`CircularSection`
771 +
        `CircularSection` object with specified parameters. The section is
772 +
        automatically connected to the underlying `Cube` data source if the
773 +
        :obj:`register_section` method of a `Cube` is used to set up the
774 +
        section, or the `Cube` is passed as the first positional argument
775 +
        during instantiation.
776 +
777 +
    Examples
778 +
    --------
779 +
780 +
    To create a `CircularSection` that is registered to a `DataCube` with
781 +
    radius ``=30``, and using the default `origin` options:
782 +
783 +
    .. plot::
784 +
        :include-source:
785 +
786 +
        >>> rcm8cube = dm.sample_data.cube.rcm8()
787 +
        >>> rcm8cube.register_section('circular', dm.section.CircularSection(radius=30))
788 +
        >>>
789 +
        >>> # show the location and the "velocity" variable
790 +
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
791 +
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
792 +
        >>> rcm8cube.sections['circular'].show_trace('r--', ax=ax[0])
793 +
        >>> rcm8cube.sections['circular'].show('velocity', ax=ax[1])
794 +
        >>> plt.show()
795 +
    """
796 +
797 +
    def __init__(self, *args, radius=None, origin=None, **kwargs):
798 +
799 +
        self._input_radius = radius
800 +
        self._input_origin = origin
801 +
        super().__init__('circular', *args, **kwargs)
802 +
803 +
    def _compute_section_coords(self):
804 +
        if (self._input_radius is None):
805 +
            self.radius = int(np.min(self.cube.shape[1:]) / 2)
806 +
        else:
807 +
            self.radius = self._input_radius
808 +
809 +
        if (self._input_origin is None):
810 +
            land_width = np.minimum(utils.guess_land_width_from_land(
811 +
                self.cube['eta'][-1, :, 0]), 5)
812 +
            self.origin = (int(self.cube.shape[2] / 2),
813 +
                           land_width)
814 +
        else:
815 +
            self.origin = self._input_origin
816 +
817 +
        xy = utils.circle_to_cells(self.origin, self.radius)
818 +
        self._x = xy[0]
819 +
        self._y = xy[1]
820 +
821 +
521 822
class RadialSection(BaseSection):
522 823
    """Radial section object.
523 824
825 +
    Section drawn as a radial cut, located a along the line starting from
826 +
    `origin` and proceeding away in direction specified by azimuth. Specify
827 +
    the location of the radial section with :obj`azimuth` and :obj:`origin`
828 +
    keyword parameter options. The radial section trace is interpolated to the
829 +
    nearest integer model domain cells, following the a line-walking algorithm
830 +
    (:obj:`~deltametrics.utils.line_to_cells`).
831 +
832 +
    .. important::
833 +
834 +
        The `origin` parameter must be specified as cell indices (not actual x
835 +
        and y coordinate values). This is a needed patch.
836 +
837 +
    .. important::
838 +
839 +
        The `origin` attempts to detect the land width from bed elevation
840 +
        changes, but should use the value of ``L0`` recorded in the netcdf
841 +
        file, or defined in the cube.
842 +
843 +
    .. important::
844 +
845 +
        This Section type will only work for deltas with an inlet along the
846 +
        ``y=0`` line. For other delta configurations, specify a radial
847 +
        section by defining two end points and instantiating a `Section` with
848 +
        the :obj:`PathSection`.
849 +
850 +
    Parameters
851 +
    ----------
852 +
    *args : :obj:`DataCube` or `StratigraphyCube`
853 +
        The `Cube` object to link for underlying data. This option should be
854 +
        ommitted if using the :obj:`register_section` method of a `Cube`.
855 +
856 +
    azimuth : :obj:`float`, `int`, optional
857 +
        The `azimuth` of the section, directed away from the origin. If no
858 +
        value is given, the `azimuth` defaults to ``90``.
859 +
860 +
    origin : :obj:`tuple` or `list` of `int`, optional
861 +
        The `origin` of the radial section. This is the "start" of the radial
862 +
        line. If no value is given, the origin defaults to the center of the
863 +
        x-direction of the model domain, and offsets into the domain a
864 +
        distance of ``y == L0``, if these values can be determined. I.e., the
865 +
        origin defaults to be centered over the channel inlet. If no value is
866 +
        given and these values cannot be determined, the origin defaults to
867 +
        ``(0, 0)``.
868 +
869 +
    length : :obj:`float`, `int`, optional
870 +
        The length of the section (note this must be given in pixel length).
871 +
        If no value is given, the length defaults to the length required to
872 +
        reach a model boundary (if a connection to underlying `Cube` exists).
873 +
        Otherwise, length is set to ``1``.
874 +
875 +
    **kwargs
876 +
        Keyword arguments are passed to `BaseSection.__init__()`. Supported
877 +
        options are `name`.
878 +
879 +
    Returns
880 +
    -------
881 +
    section : :obj:`RadialSection`
882 +
        `RadialSection` object with specified parameters. The section is
883 +
        automatically connected to the underlying `Cube` data source if the
884 +
        :obj:`register_section` method of a `Cube` is used to set up the
885 +
        section, or the `Cube` is passed as the first positional argument
886 +
        during instantiation.
887 +
888 +
    Examples
889 +
    --------
890 +
891 +
    To create a `RadialSection` that is registered to a `DataCube` at
892 +
    specified `origin` coordinate, and spans the entire model domain:
893 +
894 +
    .. plot::
895 +
        :include-source:
896 +
897 +
        >>> rcm8cube = dm.sample_data.cube.rcm8()
898 +
        >>> rcm8cube.register_section('radial', dm.section.RadialSection(azimuth=45))
899 +
        >>>
900 +
        >>> # show the location and the "velocity" variable
901 +
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
902 +
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
903 +
        >>> rcm8cube.sections['radial'].show_trace('r--', ax=ax[0])
904 +
        >>> rcm8cube.sections['radial'].show('velocity', ax=ax[1])
905 +
        >>> plt.show()
524 906
    """
907 +
    def __init__(self, *args, azimuth=None, origin=None, length=None,
908 +
                 **kwargs):
909 +
        self._input_azimuth = azimuth
910 +
        self._input_origin = origin
911 +
        self._input_length = length
912 +
        super().__init__('radial', *args, **kwargs)
525 913
526 -
    def __init__(self, apex, radius):
527 -
        raise NotImplementedError
914 +
    def _compute_section_coords(self):
915 +
916 +
        # determine the azimuth
917 +
        if (self._input_azimuth is None):
918 +
            self.azimuth = 90
919 +
        else:
920 +
            self.azimuth = self._input_azimuth
921 +
922 +
        # determine the origin of the line
923 +
        if (self._input_origin is None):
924 +
            land_width = np.minimum(utils.guess_land_width_from_land(
925 +
                self.cube['eta'][-1, :, 0]), 5)
926 +
            self.origin = (int(self.cube.shape[2] / 2),
927 +
                           land_width)
928 +
        else:
929 +
            self.origin = self._input_origin
930 +
931 +
        # determine the length of the line to travel
932 +
        # find the line of the azimuth
933 +
        theta = self.azimuth
934 +
        m = np.tan(theta * np.pi / 180)
935 +
        b = self.origin[1] - m * self.origin[0]
936 +
        if (self._input_length is None):
937 +
            # if no input
938 +
            # find the intersection with an edge
939 +
            if self.azimuth <= 90.0 and self.azimuth >= 0:
940 +
                dx = (self.cube.W - self.origin[0])
941 +
                dy = (np.tan(theta * np.pi / 180) * dx)
942 +
                if dy <= self.cube.L:
943 +
                    end_y = int(np.minimum(m * (self.cube.W) + b, self.cube.L - 1))
944 +
                    end_point = (self.cube.W - 1, end_y)
945 +
                else:
946 +
                    end_x = int(np.minimum((self.cube.L - b) / m, self.cube.W - 1))
947 +
                    end_point = (end_x, self.cube.L - 1)
948 +
            elif self.azimuth > 90 and self.azimuth <= 180:
949 +
                dx = (self.origin[0])
950 +
                dy = (np.tan(theta * np.pi / 180) * dx)
951 +
                if np.abs(dy) <= self.cube.L:
952 +
                    end_y = b
953 +
                    end_point = (0, end_y)
954 +
                else:
955 +
                    end_x = int(np.maximum((self.cube.L - b) / m,
956 +
                                           0))
957 +
                    end_point = (end_x, self.cube.L - 1)
958 +
            else:
959 +
                raise ValueError('Azimuth must be in range (0, 180).')
960 +
        else:
961 +
            # if input length
962 +
            _len = self._input_length
963 +
            # use vector math to determine end point len along azimuth
964 +
            #   vector is from (0, b) to (origin)
965 +
            vec = np.array([self.origin[0] - 0, self.origin[1] - b])
966 +
            vec_norm = vec / np.sqrt(vec**2)
967 +
            end_point = (self.origin[0] + _len*vec_norm[0],
968 +
                         self.origin[1] + _len*vec_norm[1])
969 +
970 +
        xy = utils.line_to_cells(self.origin, end_point)
971 +
        self._x = xy[0]
972 +
        self._y = xy[1]
Files Coverage
deltametrics 91.02%
Project Totals (8 files) 91.02%
codecov-umbrella
Build #329817275 -
PYTHON=3.6
OS=ubuntu-latest
codecov-umbrella
Build #329817275 -
PYTHON=3.8
OS=ubuntu-latest
codecov-umbrella
Build #329817275 -
PYTHON=3.8
OS=macos-latest
codecov-umbrella
Build #329817275 -
PYTHON=3.7
OS=ubuntu-latest
codecov-umbrella
Build #329817275 -
PYTHON=3.7
OS=macos-latest
codecov-umbrella
Build #329817275 -
PYTHON=3.8
OS=windows-latest
codecov-umbrella
Build #329817275 -
PYTHON=3.6
OS=macos-latest

No yaml found.

Create your codecov.yml to customize your Codecov experience

Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading