1
"""Classes and methods to create masks of planform features and attributes."""
2 9
import numpy as np
3 9
import matplotlib.pyplot as plt
4 9
from scipy.spatial import ConvexHull
5 9
from shapely.geometry import Point
6 9
from shapely.geometry.polygon import Polygon
7 9
from skimage import feature
8 9
from skimage import morphology
9 9
from skimage import measure
10

11 9
from . import utils
12

13

14 9
class BaseMask(object):
15
    """Low-level base class to be inherited by all mask objects."""
16

17 9
    def __init__(self, mask_type, data):
18
        """Initialize the base mask attributes and methods."""
19 9
        self.mask_type = mask_type
20
        # try to convert to np.ndarray if data not provided in that form
21 9
        if type(data) is not np.ndarray:
22 9
            try:
23 9
                self.data = data.__array__()
24 9
            except Exception:
25 9
                raise TypeError('Input data type must be numpy.ndarray,'
26
                                'but was ' + str(type(data)))
27
        else:
28 9
            self.data = data
29

30
        # set data to be 3D even if it is not (assuming t-x-y)
31 9
        if len(np.shape(self.data)) == 3:
32 9
            pass
33 9
        elif len(np.shape(self.data)) == 2:
34 9
            self.data = np.reshape(self.data, [1,
35
                                               np.shape(self.data)[0],
36
                                               np.shape(self.data)[1]])
37
        else:
38 9
            raise ValueError('Input data shape was not 2-D nor 3-D')
39

40 9
        self._mask = np.zeros(self.data.shape)
41

42 9
    @property
43 3
    def data(self):
44
        """ndarray : Values of the mask object.
45

46
        In setter, we should sanitize the inputs (enforce range 0-1) and
47
        convert everything to uints for speed and size.
48
        """
49 9
        return self._data
50

51 9
    @data.setter
52 3
    def data(self, var):
53 9
        self._data = var
54

55 9
    @property
56 3
    def mask(self):
57
        """ndarray : Binary mask values.
58

59
        Read-only mask attribute.
60
        """
61 9
        return self._mask
62

63 9
    def show(self, t=-1, ax=None, **kwargs):
64
        """Show the mask.
65

66
        Parameters
67
        ----------
68
        t : int, optional
69
            Value of time slice or integer from first dimension in 3D (t-x-y)
70
            convention to display the mask from. Default is -1 so the 'final'
71
            mask in time is displayed if this argument is not supplied.
72

73
        Passes `**kwargs` to ``matplotlib.imshow``.
74
        """
75 9
        if not ax:
76 9
            ax = plt.gca()
77 9
        cmap = kwargs.pop('cmap', 'gray')
78 9
        if hasattr(self, 'mask') and np.sum(self.mask) > 0:
79 9
            ax.imshow(self.mask[t, :, :], cmap=cmap, **kwargs)
80 9
            ax.set_title('A ' + self.mask_type + ' mask')
81
        else:
82 9
            raise AttributeError('No mask computed and/or mask is all 0s.')
83 9
        plt.draw()
84

85

86 9
class OAM(object):
87
    """Class for methods related to the Opening Angle Method [1]_."""
88

89 9
    def __init__(self, mask_type, data):
90
        """Require the same inputs as the :obj:`BaseMask` class."""
91 9
        self.mask_type = mask_type
92 9
        self.data = data
93

94 9
    def compute_shoremask(self, angle_threshold=75, **kwargs):
95
        """Compute the shoreline mask.
96

97
        Applies the opening angle method [1]_ to compute the shoreline mask.
98
        Particular method has been translated from [2]_.
99

100
        .. [1] Shaw, John B., et al. "An image‐based method for
101
           shoreline mapping on complex coasts." Geophysical Research Letters
102
           35.12 (2008).
103

104
        .. [2] Liang, Man, Corey Van Dyk, and Paola Passalacqua.
105
           "Quantifying the patterns and dynamics of river deltas under
106
           conditions of steady forcing and relative sea level rise." Journal
107
           of Geophysical Research: Earth Surface 121.2 (2016): 465-496.
108

109
        Parameters
110
        ----------
111
        angle_threshold : int, optional
112
            Threshold opening angle used in the OAM. Default is 75 degrees.
113

114
        Other Parameters
115
        ----------------
116
        topo_threshold : float, optional
117
            Threshold depth to use for the OAM. Default is -0.5.
118

119
        numviews : int, optional
120
            Defines the number of times to 'look' for the OAM. Default is 3.
121

122
        """
123
        # pop the kwargs
124 9
        self.topo_threshold = kwargs.pop('topo_threshold', -0.5)
125 9
        self.numviews = kwargs.pop('numviews', 3)
126

127
        # write the parameters to self
128 9
        self.angle_threshold = angle_threshold
129

130
        # initialize arrays
131 9
        self.oceanmap = np.zeros_like(self.data)
132 9
        self.shore_image = np.zeros_like(self.data)
133

134
        # loop through the time dimension
135 9
        for tval in range(0, self.data.shape[0]):
136 9
            trim_idx = utils.guess_land_width_from_land(self.data[tval, :, 0])
137 9
            data_trim = self.data[tval, trim_idx:, :]
138
            # use topo_threshold to identify oceanmap
139 9
            omap = (data_trim < self.topo_threshold) * 1.
140
            # if all ocean, there is no shore to be found - do next t-slice
141 9
            if (omap == 1).all():
142 9
                pass
143
            else:
144
                # apply seaangles function
145 9
                _, seaangles = self.Seaangles_mod(self.numviews, omap)
146
                # translate flat seaangles values to the shoreline image
147 9
                shore_image = np.zeros_like(data_trim)
148 9
                flat_inds = list(map(lambda x: np.ravel_multi_index(x,
149
                                                            shore_image.shape),
150
                                     seaangles[:2, :].T.astype(int)))
151 9
                shore_image.flat[flat_inds] = seaangles[-1, :]
152
                # grab contour from seaangles corresponding to angle threshold
153 9
                cs = measure.find_contours(shore_image,
154
                                           np.float(self.angle_threshold))
155 9
                C = cs[0]
156
                # convert this extracted contour to the shoreline mask
157 9
                shoremap = np.zeros_like(data_trim)
158 9
                flat_inds = list(map(lambda x: np.ravel_multi_index(x,
159
                                                            shoremap.shape),
160
                                     np.round(C).astype(int)))
161 9
                shoremap.flat[flat_inds] = 1
162
                # write shoreline map out to data.mask
163 9
                self._mask[tval, trim_idx:, :] = shoremap
164
                # assign shore_image to the mask object with proper size
165 9
                self.shore_image[tval, trim_idx:, :] = shore_image
166
                # properly assign the oceanmap to the self.oceanmap
167 9
                self.oceanmap[tval, trim_idx:, :] = omap
168

169 9
    def Seaangles_mod(self, numviews, thresholdimg):
170
        """Extract the opening angle map from an image.
171

172
        Adapted from the Matlab implementation in [2]_. Takes an image
173
        and extracts its opening angle map.
174

175
        Parameters
176
        ----------
177
        numviews : int
178
            Defines the number of times to 'look' for the opening angle map.
179

180
        thresholdimg : ndarray
181
            Binary image that has been thresholded to split deep water/land.
182

183
        Returns
184
        -------
185
        shoreangles : ndarray
186
            Flattened values corresponding to the shoreangle detected for each
187
            'look' of the opening angle method
188

189
        seaangles : ndarray
190
            Flattened values corresponding to the 'sea' angle detected for each
191
            'look' of the opening angle method. The 'sea' region is the convex
192
            hull which envelops the shoreline as well as the delta interior.
193

194
        """
195 9
        Sx, Sy = np.gradient(thresholdimg)
196 9
        G = (Sx**2 + Sy**2)**.5
197

198
        # threshold the gradient to produce edges
199 9
        edges = (G > 0) & (thresholdimg > 0)
200

201
        # extract coordinates of the edge pixels and define convex hull
202 9
        bordermap = np.pad(np.zeros_like(edges), 1, 'edge')
203 9
        bordermap[:-2, 1:-1] = edges
204 9
        bordermap[0, :] = 1
205 9
        points = np.fliplr(np.array(np.where(edges > 0)).T)
206 9
        hull = ConvexHull(points, qhull_options='Qc')
207 9
        polygon = Polygon(points[hull.vertices]).buffer(0.01)
208

209
        # identify set of points to evaluate
210 9
        sea = np.fliplr(np.array(np.where(thresholdimg > 0.5)).T)
211 9
        points_to_test = [Point(i[0], i[1]) for i in sea]
212

213
        # identify set of points in both the convex hull polygon and
214
        # defined as points_to_test and put these binary points into seamap
215 9
        In = np.array(list(map(lambda pt: polygon.contains(pt),
216
                               points_to_test)))
217 9
        Shallowsea_ = sea[In]
218 9
        seamap = np.zeros(bordermap.shape)
219 9
        flat_inds = list(map(lambda x: np.ravel_multi_index(x, seamap.shape),
220
                             np.fliplr(Shallowsea_)))
221 9
        seamap.flat[flat_inds] = 1
222 9
        seamap[:3, :] = 0
223

224
        # define other points as these 'Deepsea' points
225 9
        Deepsea_ = sea[~In]
226 9
        Deepsea = np.zeros((numviews+2, len(Deepsea_)))
227 9
        Deepsea[:2, :] = np.flipud(Deepsea_.T)
228 9
        Deepsea[-1, :] = 200.  # 200 is a background value for waves1s later
229

230
        # define points for the shallow sea and the shoreborder
231 9
        Shallowsea = np.array(np.where(seamap > 0.5))
232 9
        shoreandborder = np.array(np.where(bordermap > 0.5))
233 9
        c1 = len(Shallowsea[0])
234 9
        maxtheta = np.zeros((numviews, c1))
235

236
        # compute angle between each shallowsea and shoreborder point
237 9
        for i in range(c1):
238

239 9
            diff = shoreandborder - Shallowsea[:, i, np.newaxis]
240 9
            x = diff[0]
241 9
            y = diff[1]
242

243 9
            angles = np.arctan2(x, y)
244 9
            angles = np.sort(angles) * 180. / np.pi
245

246 9
            dangles = angles[1:] - angles[:-1]
247 9
            dangles = np.concatenate((dangles,
248
                                      [360 - (angles.max() - angles.min())]))
249 9
            dangles = np.sort(dangles)
250

251 9
            maxtheta[:, i] = dangles[-numviews:]
252

253
        # set up arrays for tracking the shore points and  their angles
254 9
        allshore = np.array(np.where(edges > 0))
255 9
        c3 = len(allshore[0])
256 9
        maxthetashore = np.zeros((numviews, c3))
257

258
        # get angles between the shore points and shoreborder points
259 9
        for i in range(c3):
260

261 9
            diff = shoreandborder - allshore[:, i, np.newaxis]
262 9
            x = diff[0]
263 9
            y = diff[1]
264

265 9
            angles = np.arctan2(x, y)
266 9
            angles = np.sort(angles) * 180. / np.pi
267

268 9
            dangles = angles[1:] - angles[:-1]
269 9
            dangles = np.concatenate((dangles,
270
                                      [360 - (angles.max() - angles.min())]))
271 9
            dangles = np.sort(dangles)
272

273 9
            maxthetashore[:, i] = dangles[-numviews:]
274

275
        # define the shoreangles and seaangles identified
276 9
        shoreangles = np.vstack([allshore, maxthetashore])
277 9
        seaangles = np.hstack([np.vstack([Shallowsea, maxtheta]), Deepsea])
278

279 9
        return shoreangles, seaangles
280

281

282 9
class ChannelMask(BaseMask, OAM):
283
    """Identify a binary channel mask.
284

285
    A channel mask object, helps enforce valid masking of channels.
286

287
    Examples
288
    --------
289
    Initialize the channel mask
290
        >>> velocity = rcm8cube['velocity'][-1, :, :]
291
        >>> topo = rcm8cube['eta'][-1, :, :]
292
        >>> cmsk = dm.mask.ChannelMask(velocity, topo)
293

294
    And visualize the mask:
295
        >>> cmsk.show()
296

297
    .. plot:: mask/channelmask.py
298

299
    """
300

301 9
    def __init__(self,
302
                 velocity,
303
                 topo,
304
                 velocity_threshold=0.3,
305
                 angle_threshold=75,
306
                 is_mask=False,
307
                 **kwargs):
308
        """Initialize the ChannelMask.
309

310
        Intializing the channel mask requires a flow velocity field and an
311
        array of the delta topography.
312

313
        Parameters
314
        ----------
315
        velocity : ndarray
316
            The velocity array to be used for mask creation.
317

318
        topo : ndarray
319
            The model topography to be used for mask creation.
320

321
        velocity_threshold : float, optional
322
            Threshold velocity above which flow is considered 'channelized'.
323
            The default value is 0.3 m/s based on DeltaRCM default parameters
324
            for sediment transport.
325

326
        angle_threshold : int, optional
327
            Threshold opening angle used in the OAM. Default is 75 degrees.
328

329
        is_mask : bool, optional
330
            Whether the data in :obj:`arr` is already a binary mask. Default is
331
            False. This should be set to True, if you have already binarized
332
            the data yourself, using custom routines, and want to just store
333
            the data in the ChannelMask object.
334

335
        Other Parameters
336
        ----------------
337
        landmask : :obj:`LandMask`, optional
338
            A :obj:`LandMask` object with a defined binary land mask.
339
            If given, it will be used to help define the channel mask.
340

341
        wetmask : :obj:`WetMask`, optional
342
            A :obj:`WetMask` object with a defined binary wet mask.
343
            If given, the landmask attribute it contains will be used to
344
            determine the channel mask.
345

346
        kwargs : optional
347
            Keyword arguments for :obj:`compute_shoremask`.
348

349
        """
350 9
        super().__init__(mask_type='channel', data=topo)
351 9
        if type(velocity) is np.ndarray:
352 9
            self.velocity = velocity
353
        else:
354 9
            try:
355 9
                self.velocity = velocity.__array__()
356 9
            except Exception:
357 9
                raise TypeError('Input velocity parameter is invalid. Must be'
358
                                'a np.ndarray or a'
359
                                'deltametrics.cube.CubeVariable object but it'
360
                                'was a: ' + str(type(velocity)))
361

362
        # assign **kwargs to self incase there are masks
363 9
        for key, value in kwargs.items():
364 9
            setattr(self, key, value)
365

366
        # write parameter values to self so they don't get lost
367 9
        self.velocity_threshold = velocity_threshold
368 9
        self.angle_threshold = angle_threshold
369

370 9
        if is_mask is False:
371 9
            self.compute_channelmask(**kwargs)
372 9
        elif is_mask is True:
373 9
            self._mask = self.data
374
        else:
375 9
            raise TypeError('is_mask must be a `bool`,'
376
                            'but was: ' + str(type(is_mask)))
377

378 9
    def compute_channelmask(self, **kwargs):
379
        """Compute the ChannelMask.
380

381
        Either uses an existing landmask, or recomputes it to use with the flow
382
        velocity threshold to calculate a binary channelmask.
383

384
        Other Parameters
385
        ----------------
386
        kwargs : optional
387
            Keyword arguments for :obj:`compute_shoremask`.
388

389
        """
390
        # check if a 'landmask' is available
391 9
        if hasattr(self, 'landmask'):
392 9
            try:
393 9
                self.oceanmap = self.landmask.oceanmap
394 9
                self.landmask = self.landmask.mask
395 9
            except Exception:
396 9
                self.compute_shoremask(self.angle_threshold, **kwargs)
397 9
                self.landmask = (self.shore_image < self.angle_threshold) * 1
398 9
        elif hasattr(self, 'wetmask'):
399 9
            try:
400 9
                self.landmask = self.wetmask.landmask
401 9
            except Exception:
402 9
                self.compute_shoremask(self.angle_threshold, **kwargs)
403 9
                self.landmask = (self.shore_image < self.angle_threshold) * 1
404
        else:
405 9
            self.compute_shoremask(self.angle_threshold, **kwargs)
406 9
            self.landmask = (self.shore_image < self.angle_threshold) * 1
407

408
        # compute a flowmap of cells where flow exceeds the velocity threshold
409 9
        self.flowmap = (self.velocity > self.velocity_threshold) * 1.
410

411
        # calculate the channelmask as the cells exceeding the threshold
412
        # within the topset of the delta (ignoring flow in ocean)
413 9
        self._mask = np.minimum(self.landmask, self.flowmap)
414

415

416 9
class WetMask(BaseMask, OAM):
417
    """Compute the wet mask.
418

419
    A wet mask object, identifies all wet pixels on the delta topset. Starts
420
    with the land mask and then uses the topo_threshold defined for the
421
    shoreline computation to add the wet pixels on the topset back to the mask.
422

423
    If a land mask has already been computed, then it can be used to define the
424
    wet mask. Otherwise the wet mask can be computed from scratch.
425

426
    Examples
427
    --------
428
    Initialize the wet mask
429
        >>> arr = rcm8cube['eta'][-1, :, :]
430
        >>> wmsk = dm.mask.WetMask(arr)
431

432
    And visualize the mask:
433
        >>> wmsk.show()
434

435
    .. plot:: mask/wetmask.py
436

437
    """
438

439 9
    def __init__(self,
440
                 arr,
441
                 angle_threshold=75,
442
                 is_mask=False,
443
                 **kwargs):
444
        """Initialize the WetMask.
445

446
        Intializing the wet mask requires either a 2-D array of data, or it
447
        can be computed if a :obj:`LandMask` has been previously computed.
448

449
        Parameters
450
        ----------
451
        arr : ndarray
452
            The data array to make the mask from.
453

454
        angle_threshold : int, optional
455
            Threshold opening angle used in the OAM. Default is 75 degrees.
456

457
        is_mask : bool, optional
458
            Whether the data in :obj:`arr` is already a binary mask. Default is
459
            False. This should be set to True, if you have already binarized
460
            the data yourself, using custom routines, and want to just store
461
            the data in the WetMask object.
462

463
        Other Parameters
464
        ----------------
465
        landmask : :obj:`LandMask`, optional
466
            A :obj:`LandMask` object with a defined binary shoreline mask.
467
            If given, the :obj:`LandMask` object will be checked for the
468
            `shore_image` and `angle_threshold` attributes.
469

470
        kwargs : optional
471
            Keyword arguments for :obj:`compute_shoremask`.
472

473
        """
474 9
        super().__init__(mask_type='wet', data=arr)
475

476
        # assign **kwargs to self in case the landmask was passed
477 9
        for key, value in kwargs.items():
478 9
            setattr(self, key, value)
479

480
        # put angle_threshold in self
481 9
        self.angle_threshold = angle_threshold
482

483 9
        if is_mask is False:
484 9
            self.compute_wetmask(**kwargs)
485 9
        elif is_mask is True:
486 9
            self._mask = self.data
487
        else:
488 9
            raise TypeError('is_mask must be a `bool`,'
489
                            'but was: ' + str(type(is_mask)))
490

491 9
    def compute_wetmask(self, **kwargs):
492
        """Compute the WetMask.
493

494
        Either recomputes the landmask, or uses precomputed information from
495
        the landmask to create the wetmask.
496

497
        """
498
        # check if a 'landmask' was passed in
499 9
        if hasattr(self, 'landmask'):
500 9
            try:
501 9
                self.oceanmap = self.landmask.oceanmap
502 9
                self.landmask = self.landmask.mask
503 9
            except Exception:
504 9
                self.compute_shoremask(self.angle_threshold, **kwargs)
505 9
                self.landmask = (self.shore_image < self.angle_threshold) * 1
506
        else:
507 9
            self.compute_shoremask(self.angle_threshold, **kwargs)
508 9
            self.landmask = (self.shore_image < self.angle_threshold) * 1
509
        # set wet mask as data.mask
510 9
        self._mask = self.oceanmap * self.landmask
511

512

513 9
class LandMask(BaseMask, OAM):
514
    """Identify a binary mask of the delta topset.
515

516
    A land mask object, helps enforce valid masking of delta topset.
517

518
    If a shoreline mask has been computed, it can be used to help compute the
519
    land mask, otherwise it will be computed from scratch.
520

521
    Examples
522
    --------
523
    Initialize the mask.
524
        >>> arr = rcm8cube['eta'][-1, :, :]
525
        >>> lmsk = dm.mask.LandMask(arr)
526

527
    And visualize the mask:
528
        >>> lmsk.show()
529

530
    .. plot:: mask/landmask.py
531

532
    """
533

534 9
    def __init__(self,
535
                 arr,
536
                 angle_threshold=75,
537
                 is_mask=False,
538
                 **kwargs):
539
        """Initialize the LandMask.
540

541
        Intializing the land mask requires an array of data, should be
542
        two-dimensional. If a shoreline mask (:obj:`ShoreMask`) has been
543
        computed, then this can be used to define the land mask. Otherwise the
544
        necessary computations will occur from scratch.
545

546
        Parameters
547
        ----------
548
        arr : ndarray
549
            2-D topographic array to make the mask from.
550

551
        angle_threshold : int, optional
552
            Threshold opening angle used in the OAM. Default is 75 degrees.
553

554
        is_mask : bool, optional
555
            Whether the data in :obj:`arr` is already a binary mask. Default
556
            value is False. This should be set to True, if you have already
557
            binarized the data yourself, using custom routines, and want to
558
            just store the data in the LandMask object.
559

560
        Other Parameters
561
        ----------------
562
        shoremask : :obj:`ShoreMask`, optional
563
            A :obj:`ShoreMask` object with a defined binary shoreline mask.
564
            If given, the :obj:`ShoreMask` object will be checked for the
565
            `shore_image` and `angle_threshold` attributes.
566

567
        kwargs : optional
568
            Keyword arguments for :obj:`compute_shoremask`.
569

570
        """
571 9
        super().__init__(mask_type='land', data=arr)
572

573
        # assign **kwargs to self in event a mask was passed
574 9
        for key, value in kwargs.items():
575 9
            setattr(self, key, value)
576

577
        # set parameter value to self so it is kept
578 9
        self.angle_threshold = angle_threshold
579

580 9
        if not is_mask:
581 9
            self.compute_landmask(**kwargs)
582 9
        elif is_mask is True:
583 9
            self._mask = self.data
584
        else:
585 9
            raise TypeError('is_mask must be a `bool`,'
586
                            'but was: ' + str(type(is_mask)))
587

588 9
    def compute_landmask(self, **kwargs):
589
        """Compute the LandMask.
590

591
        Uses data from the shoreline computation or recomputes the shoreline
592
        using the opening angle method, then an angle threshold is used to
593
        identify the land.
594

595
        """
596
        # check if a `shoremask' was passed in
597 9
        if hasattr(self, 'shoremask'):
598 9
            try:
599 9
                self.shore_image = self.shoremask.shore_image
600 9
                self.angle_threshold = self.shoremask.angle_threshold
601 9
            except Exception:
602 9
                self.compute_shoremask(self.angle_threshold, **kwargs)
603
        else:
604 9
            self.compute_shoremask(self.angle_threshold, **kwargs)
605
        # set the land mask to data.mask
606 9
        self._mask = (self.shore_image < self.angle_threshold) * 1
607

608

609 9
class ShorelineMask(BaseMask, OAM):
610
    """Identify the shoreline as a binary mask.
611

612
    A shoreline mask object, provides a binary identification of shoreline
613
    pixels.
614

615
    Examples
616
    --------
617
    Initialize the mask.
618
        >>> arr = rcm8cube['eta'][-1, :, :]
619
        >>> shrmsk = dm.mask.ShorelineMask(arr)
620

621
    Visualize the mask:
622
        >>> shrmsk.show()
623

624
    .. plot:: mask/shoremask.py
625

626
    """
627

628 9
    def __init__(self,
629
                 arr,
630
                 angle_threshold=75,
631
                 is_mask=False,
632
                 **kwargs):
633
        """Initialize the ShorelineMask.
634

635
        Initializing the shoreline mask requires a 2-D array of data. The
636
        shoreline mask is computed using the opening angle method (OAM)
637
        described in [1]_. The default parameters are designed for
638
        DeltaRCM and follow the implementation described in [2]_.
639

640
        Parameters
641
        ----------
642
        arr : ndarray
643
            2-D topographic array to make the mask from.
644

645
        angle_threshold : int, optional
646
            Threshold opening angle used in the OAM. Default is 75 degrees.
647

648
        is_mask : bool, optional
649
            Whether the data in :obj:`arr` is already a binary mask. Default
650
            value is False. This should be set to True, if you have already
651
            binarized the data yourself, using custom routines, and want to
652
            just store the data in the ShorelineMask object.
653

654
        Other Parameters
655
        ----------------
656
        kwargs : optional
657
            Keyword arguments for :obj:`compute_shoremask`.
658

659
        """
660 9
        super().__init__(mask_type='shore', data=arr)
661

662 9
        if is_mask is False:
663 9
            self.compute_shoremask(angle_threshold, **kwargs)
664 9
        elif is_mask is True:
665 9
            self._mask = self.data
666
        else:
667 9
            raise TypeError('is_mask must be a `bool`,'
668
                            'but was: ' + str(type(is_mask)))
669

670

671 9
class EdgeMask(BaseMask, OAM):
672
    """Identify the land-water edges.
673

674
    An edge mask object, delineates the boundaries between land and water.
675

676
    Examples
677
    --------
678
    Initialize the edge mask
679
        >>> arr = rcm8cube['eta'][-1, :, :]
680
        >>> emsk = dm.mask.EdgeMask(arr)
681

682
    Visualize the mask:
683
        >>> emsk.show()
684

685
    .. plot:: mask/edgemask.py
686

687
    """
688

689 9
    def __init__(self,
690
                 arr,
691
                 angle_threshold=75,
692
                 is_mask=False,
693
                 **kwargs):
694
        """Initialize the EdgeMask.
695

696
        Initializing the edge mask requires either a 2-D array of topographic
697
        data, or it can be computed using the :obj:`LandMask` and the
698
        :obj:`WetMask`.
699

700
        Parameters
701
        ----------
702
        arr : ndarray
703
            The data array to make the mask from.
704

705
        angle_threshold : int, optional
706
            Threshold opening angle used in the OAM. Default is 75 degrees.
707

708
        is_mask : bool, optional
709
            Whether the data in :obj:`arr` is already a binary mask. Default
710
            value is False. This should be set to True, if you have already
711
            binarized the data yourself, using custom routines, and want to
712
            just store the data in the EdgeMask object.
713

714
        Other Parameters
715
        ----------------
716
        landmask : :obj:`LandMask`, optional
717
            A :obj:`LandMask` object with the land identified
718

719
        wetmask : :obj:`WetMask`, optional
720
            A :obj:`WetMask` object with the surface water identified
721

722
        kwargs : optional
723
            Keyword arguments for :obj:`compute_shoremask`.
724

725
        """
726 9
        super().__init__(mask_type='edge', data=arr)
727

728
        # assign **kwargs to self incase there are masks
729 9
        for key, value in kwargs.items():
730 9
            setattr(self, key, value)
731

732
        # write parameter to self so it is not lost
733 9
        self.angle_threshold = angle_threshold
734

735 9
        if is_mask is False:
736 9
            self.compute_edgemask(**kwargs)
737 9
        elif is_mask is True:
738 9
            self._mask = self.data
739
        else:
740 9
            raise TypeError('is_mask must be a `bool`,'
741
                            'but was: ' + str(type(is_mask)))
742

743 9
    def compute_edgemask(self, **kwargs):
744
        """Compute the EdgeMask.
745

746
        Either computes it from just the topographic array, or uses information
747
        from :obj:`LandMask` and :obj:`WetMask` to create the edge mask.
748

749
        """
750
        # if landmask and edgemask exist it is quick
751 9
        if hasattr(self, 'landmask') and hasattr(self, 'wetmask'):
752 9
            try:
753 9
                self.landmask = self.landmask.mask
754 9
                self.wetmask = self.wetmask.mask
755 9
            except Exception:
756 9
                self.compute_shoremask(self.angle_threshold, **kwargs)
757 9
                self.landmask = (self.shore_image < self.angle_threshold) * 1
758 9
                self.wetmask = self.oceanmap * self.landmask
759 9
        elif hasattr(self, 'wetmask'):
760 9
            try:
761 9
                self.wetmask = self.wetmask.mask
762 9
                self.landmask = self.wetmask.landmask
763 9
            except Exception:
764 9
                self.compute_shoremask(self.angle_threshold, **kwargs)
765 9
                self.landmask = (self.shore_image < self.angle_threshold) * 1
766 9
                self.wetmask = self.oceanmap * self.landmask
767
        else:
768 9
            self.compute_shoremask(self.angle_threshold, **kwargs)
769 9
            self.landmask = (self.shore_image < self.angle_threshold) * 1
770 9
            self.wetmask = self.oceanmap * self.landmask
771
        # compute the edge mask
772 9
        for i in range(0, np.shape(self._mask)[0]):
773 9
            self._mask[i, :, :] = np.maximum(0,
774
                                             feature.canny(self.wetmask[i,
775
                                                                        :,
776
                                                                        :])*1 -
777
                                             feature.canny(self.landmask[i,
778
                                                                         :,
779
                                                                         :])*1)
780

781

782 9
class CenterlineMask(BaseMask):
783
    """Identify channel centerline mask.
784

785
    A centerline mask object, provides the location of channel centerlines.
786

787
    Examples
788
    --------
789
    Initialize the centerline mask
790
        >>> channelmask = dm.mask.ChannelMask(rcm8cube['velocity'][-1, :, :],
791
                                              rcm8cube['eta'][-1, :, :])
792
        >>> clmsk = dm.mask.CenterlineMask(channelmask)
793

794
    Visualize the mask
795
        >>> clmsk.show()
796

797
    .. plot:: mask/centerlinemask.py
798

799
    """
800

801 9
    def __init__(self,
802
                 channelmask,
803
                 is_mask=False,
804
                 method='skeletonize',
805
                 **kwargs):
806
        """Initialize the CenterlineMask.
807

808
        Initialization of the centerline mask object requires a 2-D channel
809
        mask (can be the :obj:`ChannelMask` object or a binary 2-D array).
810

811
        Parameters
812
        ----------
813
        channelmask : :obj:`ChannelMask` or ndarray
814
            The channel mask to derive the centerlines from
815

816
        is_mask : bool, optional
817
            Whether the data in :obj:`arr` is already a binary mask. Default
818
            value is False. This should be set to True, if you have already
819
            binarized the data yourself, using custom routines, and want to
820
            just store the data in the CenterlineMask object.
821

822
        method : str, optional
823
            The method to use for the centerline mask computation. The default
824
            method ('skeletonize') is a morphological skeletonization of the
825
            channel mask.
826

827
        Other Parameters
828
        ----------------
829
        kwargs : optional
830
            Keyword arguments for the 'rivamap' functionality.
831

832
        """
833 9
        if isinstance(channelmask, np.ndarray):
834 9
            super().__init__(mask_type='centerline', data=channelmask)
835 9
        elif isinstance(channelmask, ChannelMask):
836 9
            super().__init__(mask_type='centerline',
837
                             data=channelmask.mask)
838
        else:
839 9
            raise TypeError('Input channelmask parameter is invalid. Must be'
840
                            'a np.ndarray or ChannelMask object but it was'
841
                            'a: ' + str(type(channelmask)))
842

843
        # save method type value to self
844 9
        self.method = method
845

846 9
        if is_mask is False:
847 9
            self.compute_centerlinemask(**kwargs)
848 9
        elif is_mask is True:
849 9
            self._mask = self.data
850
        else:
851 9
            raise TypeError('is_mask must be a `bool`,'
852
                            'but was: ' + str(type(is_mask)))
853

854 9
    def compute_centerlinemask(self, **kwargs):
855
        """Compute the centerline mask.
856

857
        Function for computing the centerline mask. The default implementation
858
        is a morphological skeletonization operation using the
859
        `skimage.morphology.skeletonize` function.
860

861
        Alternatively, the method of  centerline extraction based on non-maxima
862
        suppression of the singularity index, as described in [3]_ can be
863
        specified. This requires the optional dependency `RivaMap`_.
864

865
        .. [3] Isikdogan, Furkan, Alan Bovik, and Paola Passalacqua. "RivaMap:
866
               An automated river analysis and mapping engine." Remote Sensing
867
               of Environment 202 (2017): 88-97.
868

869
        .. _Rivamap: https://github.com/isikdogan/rivamap
870

871
        Other Parameters
872
        ----------------
873
        minScale : float, optional
874
            Minimum scale to use for the singularity index extraction, see [3]_
875

876
        nrScales : int, optional
877
            Number of scales to use for singularity index, see [3]_
878

879
        nms_threshold : float between 0 and 1, optional
880
            Threshold to convert the non-maxima suppression results into a
881
            binary mask. Default value is 0.1 which means that the lowest 10%
882
            non-maxima suppression values are ignored when making the binary
883
            centerline mask.
884

885
        """
886
        # skimage.morphology.skeletonize() method
887 9
        if self.method == 'skeletonize':
888 9
            for i in range(0, np.shape(self._mask)[0]):
889 9
                self._mask[i, :, :] = morphology.skeletonize(self.data[i, :, :])
890

891 9
        if self.method == 'rivamap':
892
            # rivamap based method - first check for import error
893 9
            try:
894 9
                from rivamap.singularity_index import applyMMSI as MMSI
895 0
                from rivamap.singularity_index import SingularityIndexFilters as SF
896 0
                from rivamap.delineate import extractCenterlines as eCL
897 9
            except Exception:
898 9
                raise ImportError('You must install the optional dependency:'
899
                                  ' rivamap, to use this centerline extraction'
900
                                  ' method')
901

902
            # pop the kwargs
903 0
            self.minScale = kwargs.pop('minScale', 1.5)
904 0
            self.nrScales = kwargs.pop('nrScales', 12)
905 0
            self.nms_threshold = kwargs.pop('nms_threshold', 0.1)
906

907
            # now do the computation - first change type and do psi extraction
908 0
            if self.data.dtype == 'int64':
909 0
                self.data = self.data.astype('float')/(2**64 - 1)
910 0
            self.psi, widths, orient = MMSI(self.data,
911
                                            filters=SF(minScale=self.minScale,
912
                                                       nrScales=self.nrScales))
913
            # compute non-maxima suppresion then normalize/threshold to
914
            # make binary
915 0
            self.nms = eCL(orient, self.psi)
916 0
            nms_norm = self.nms/self.nms.max()
917
            # compute mask
918 0
            self._mask = (nms_norm > self.nms_threshold) * 1

Read our documentation on viewing source code .

Loading