adding connectivity option for islands

rivgraph/classes.py

@@ -123,7 +123,6 @@ def __init__(self, name, path_to_mask, results_folder=None,
         self.paths['input_mask'] = os.path.normpath(path_to_mask)
 
         # Handle georeferencing
-        print(self.paths['input_mask'])
         # GA_Update required for setting dummy projection/geotransform
         self.gdobj = gdal.Open(self.paths['input_mask'], gdal.GA_Update)
         self.imshape = (self.gdobj.RasterYSize, self.gdobj.RasterXSize)
@@ -240,7 +239,8 @@ def compute_junction_angles(self, weight=None):
 
 
     def get_islands(self, props=['area', 'maxwidth', 'major_axis_length',
-                                 'minor_axis_length', 'surrounding_links']):
+                                 'minor_axis_length', 'surrounding_links'],
+                          connectivity=2):
         """
         Finds all the islands in the binary mask and computes their morphological
         properties. Can be used to help "clean" masks of small islands. Must
@@ -252,7 +252,10 @@ def get_islands(self, props=['area', 'maxwidth', 'major_axis_length',
             Properties to compute for each island. Properties can be any of those
             provided by rivgraph.im_utils.regionprops.
             The default is ['area', 'maxwidth', 'major_axis_length', 'minor_axis_length'].
-
+        connectivity : int, optional
+            If 1, 4-connectivity will be used to determine connected blobs. If
+            2, 8-connectivity will be used. The default is 2. 
+            
         Returns
         -------
         islands : geopandas GeoDataFrame
@@ -272,7 +275,7 @@ def get_islands(self, props=['area', 'maxwidth', 'major_axis_length',
         if self.verbose is True:
             print('Getting island properties...', end='')
 
-        islands, Iislands = mu.get_island_properties(self.Imask, self.pixlen, self.pixarea, self.crs, self.gt, props)
+        islands, Iislands = mu.get_island_properties(self.Imask, self.pixlen, self.pixarea, self.crs, self.gt, props, connectivity=connectivity)
 
         if self.verbose is True:
             print('done.')

rivgraph/ln_utils.py

@@ -1518,9 +1518,6 @@ def links_to_gpd(links, gdobj):
     # Create geodataframe
     links_gpd = gpd.GeoDataFrame()
 
-    # Assign CRS
-    links_gpd.crs = CRS(gdobj.GetProjection())
-
     # Append geometries
     geoms = []
     for i, lidx in enumerate(links['idx']):
@@ -1532,5 +1529,9 @@ def links_to_gpd(links, gdobj):
     links_gpd['id'] = links['id']
     links_gpd['us node'] = [c[0] for c in links['conn']]
     links_gpd['ds node'] = [c[1] for c in links['conn']]
+
+    # Assign CRS - done last to avoid DeprecationWarning - need geometry 
+    # to exist before assigning CRS.
+    links_gpd.crs = CRS(gdobj.GetProjection())
 
     return links_gpd

rivgraph/mask_utils.py

@@ -38,10 +38,17 @@ def pixagon(c_cent, r_cent, pixlen):
     return pixgon
 
 
-def get_island_properties(Imask, pixlen, pixarea, crs, gt, props):
+def get_island_properties(Imask, pixlen, pixarea, crs, gt, props, connectivity=2):
     """Get island properties."""
+    # Imask = d.Imask.copy()
+    # pixlen = d.pixlen
+    # pixarea = d.pixarea
+    # crs = d.crs
+    # gt = d.gt
+    # props = ['area', 'maxwidth', 'major_axis_length',
+    #                              'minor_axis_length', 'surrounding_links']
+
     # maxwidth is an additional property
-
     if 'maxwidth' in props:
         props.remove('maxwidth')
         do_maxwidth = True
@@ -57,7 +64,7 @@ def get_island_properties(Imask, pixlen, pixarea, crs, gt, props):
     Imaskpad = np.array(np.pad(Imask, 1, mode='constant'), dtype=np.bool)
     Imp_invert = np.invert(Imaskpad)
 
-    rp_islands, Ilabeled = iu.regionprops(Imp_invert, props=props)
+    rp_islands, Ilabeled = iu.regionprops(Imp_invert, props=props, connectivity=connectivity)
 
     # Make polygons of the island perimeters
     # Also get ids to match the labeled image