Update mosaic shp transects workflow

CameronBodine · Feb 20, 2024 · c7964e9 · c7964e9
1 parent e8cba92
commit c7964e9
Showing 1 changed file with 81 additions and 80 deletions.
diff --git a/utils/Substrate_Summaries/00_substrate_shps_mosaic_transects.py b/utils/Substrate_Summaries/00_substrate_shps_mosaic_transects.py
@@ -23,21 +23,15 @@
 ############
 # Parameters
 
-substrateOutput = 'EGN'
+substrateOutputs = ['Raw', 'EGN']
 topDir = r'E:/SynologyDrive/GulfSturgeonProject/SSS_Data_Processed'
-transectDir = os.path.join(topDir, substrateOutput)
-outDir = os.path.join(topDir, 'Substrate_Summaries', '02_Substrate_Shps_Mosaic_Transects', substrateOutput)
 
 subShpPattern = '*map_substrate*.shp'
-
 least2mostImport = ['Other', 'Fines Flat', 'Fines Ripple', 'Hard Bottom', 'Cobble Boulder', 'Wood']
-
 outEpsg = 32616
 
 # Normalize Paths
 topDir = os.path.normpath(topDir)
-transectDir = os.path.normpath(transectDir)
-outDir = os.path.normpath(outDir)
 
 ###########
 # Functions
@@ -162,116 +156,123 @@ def getBearingLine(x, y, rr, rl, d):
 
 #########
 # Do Work
+for substrateOutput in substrateOutputs:
 
-# Prepare directories
-transectDir = os.path.abspath(os.path.normpath(transectDir))
-
-if not os.path.exists(outDir):
-    os.makedirs(outDir)
-
-# Get shps
-shpFiles = glob(os.path.join(transectDir, '**', subShpPattern), recursive=True)
-
-# Remove projName from dirs
-shpFiles = [f for f in shpFiles if outDir not in f]
+    # Prepare Paths
+    transectDir = os.path.join(topDir, substrateOutput)
+    outDir = os.path.join(topDir, 'Substrate_Summaries', '02_Substrate_Shps_Mosaic_Transects', substrateOutput)
 
+    # Normalize Paths
+    transectDir = os.path.normpath(transectDir)
+    outDir = os.path.normpath(outDir)
 
+    # Prepare directories
+    transectDir = os.path.abspath(os.path.normpath(transectDir))
 
+    if not os.path.exists(outDir):
+        os.makedirs(outDir)
 
-#############
-# shpFiles = [f for f in shpFiles if 'PRL_' in f]
+    # Get shps
+    shpFiles = glob(os.path.join(transectDir, '**', subShpPattern), recursive=True)
 
-#############
+    # Remove projName from dirs
+    shpFiles = [f for f in shpFiles if outDir not in f]
 
 
 
+    ####################################################
+    # shpFiles = [f for f in shpFiles if 'PRL_' in f]
+    ####################################################
 
 
-# Store shps in dataframe
-shpDF = pd.DataFrame()
-shpDF['shp'] = shpFiles
 
-# Add projName to df
-shpDF['proj'] = shpDF['shp'].apply(lambda d: os.path.basename(d))
 
-# Add river and code
-shpDF[['river', 'river_code']] = shpDF.apply(lambda row: get_river_code(row), axis=1).values.tolist()
 
-# Get RKM
-shpDF['up_rkm'] = shpDF['proj'].apply(lambda d: d.split('_')[1]).astype(int)
-shpDF['dn_rkm'] = shpDF['proj'].apply(lambda d: d.split('_')[2]).astype(int)
+    # Store shps in dataframe
+    shpDF = pd.DataFrame()
+    shpDF['shp'] = shpFiles
 
-# Sort by river code and up rkm
-shpDF = shpDF.sort_values(['river_code', 'up_rkm'], ascending=False).reset_index()
+    # Add projName to df
+    shpDF['proj'] = shpDF['shp'].apply(lambda d: os.path.basename(d))
 
+    # Add river and code
+    shpDF[['river', 'river_code']] = shpDF.apply(lambda row: get_river_code(row), axis=1).values.tolist()
 
-# shpDF = shpDF[1:4]  ####### For testing
+    # Get RKM
+    shpDF['up_rkm'] = shpDF['proj'].apply(lambda d: d.split('_')[1]).astype(int)
+    shpDF['dn_rkm'] = shpDF['proj'].apply(lambda d: d.split('_')[2]).astype(int)
 
+    # Sort by river code and up rkm
+    shpDF = shpDF.sort_values(['river_code', 'up_rkm'], ascending=False).reset_index()
 
 
-# Create a map based on river
-for name, group in shpDF.groupby('river_code'):
+    # shpDF = shpDF[1:4]  ####### For testing
 
-    print('\n\nWorking On:', name)
 
-    # Get shapefiles
-    shpFiles = group.shp.values.tolist()
 
-    # Iterate classes
-    for subClassName in least2mostImport:
-        print('\t', subClassName)
+    # Create a map based on river
+    for name, group in shpDF.groupby('river_code'):
 
-        allClassPolys = gpd.GeoDataFrame()
+        print('\n\nWorking On:', name)
 
-        for shpFile in shpFiles:
+        # Get shapefiles
+        shpFiles = group.shp.values.tolist()
 
-            # Open the shapefile
-            shp = gpd.read_file(shpFile)
+        # Iterate classes
+        for subClassName in least2mostImport:
+            print('\t', substrateOutput, ':', subClassName)
 
-            print(shp.crs)
+            allClassPolys = gpd.GeoDataFrame()
 
-            if shp.crs is None:
-                shp = shp.set_crs(outEpsg)
+            for shpFile in shpFiles:
 
-            # Reproject
-            shp = shp.to_crs(outEpsg)
+                print('\t\t', os.path.basename(shpFile))
 
-            # Get the class polys
-            classPolys = shp.loc[shp['Name'] == subClassName]
+                # Open the shapefile
+                shp = gpd.read_file(shpFile)
 
-            # Explode
-            classPolys = classPolys.explode()
+                if shp.crs is None:
+                    shp = shp.set_crs(outEpsg)
 
-            # Concatenate with allClass Polys
-            allClassPolys = pd.concat([allClassPolys, classPolys])
+                # Reproject
+                shp = shp.to_crs(outEpsg)
 
-            del shp
+                # Get the class polys
+                classPolys = shp.loc[shp['Name'] == subClassName]
 
-        # Calculate area
-        allClassPolys['Area_m'] = np.around(allClassPolys.geometry.area, 2)
+                # Explode
+                classPolys = classPolys.explode(index_parts=False)
 
-        # Dissolve
-        allClassPolys = allClassPolys.dissolve()
+                # Concatenate with allClass Polys
+                allClassPolys = pd.concat([allClassPolys, classPolys])
 
-        # Explode
-        allClassPolys = allClassPolys.explode().reset_index(drop=True)
+                del shp
 
-        # Overlay (opposite of clip)
-        if 'finalSubMap' not in locals():
-            finalSubMap = allClassPolys
-        else:
-            finalSubMap = finalSubMap.overlay(allClassPolys, how='difference', keep_geom_type=True)
-            finalSubMap = pd.concat([finalSubMap, allClassPolys])
+            # Calculate area
+            allClassPolys['Area_m'] = np.around(allClassPolys.geometry.area, 2)
 
-        del allClassPolys
-
-
-
-    outShp = name+'_substrate_shps_mosaic.shp' 
-    outShp = os.path.join(outDir, outShp)
-    finalSubMap.to_file(outShp)
+            # Dissolve
+            allClassPolys = allClassPolys.dissolve()
 
-    del finalSubMap
-
-    print(outShp)
+            # Explode
+            allClassPolys = allClassPolys.explode(index_parts=False).reset_index(drop=True)
+
+            # Overlay (opposite of clip)
+            if 'finalSubMap' not in locals():
+                finalSubMap = allClassPolys
+            else:
+                finalSubMap = finalSubMap.overlay(allClassPolys, how='difference', keep_geom_type=True)
+                finalSubMap = pd.concat([finalSubMap, allClassPolys])
+
+            del allClassPolys
+
+
+
+        outShp = name+'_substrate_shps_mosaic.shp' 
+        outShp = os.path.join(outDir, outShp)
+        finalSubMap.to_file(outShp)
+
+        del finalSubMap
+
+        print(outShp)