From c7964e9167d36824a5bbb5ca2bb78b0ff631a2e2 Mon Sep 17 00:00:00 2001 From: Cameron Bodine Date: Tue, 20 Feb 2024 08:59:37 -0700 Subject: [PATCH] Update mosaic shp transects workflow --- .../00_substrate_shps_mosaic_transects.py | 161 +++++++++--------- 1 file changed, 81 insertions(+), 80 deletions(-) diff --git a/utils/Substrate_Summaries/00_substrate_shps_mosaic_transects.py b/utils/Substrate_Summaries/00_substrate_shps_mosaic_transects.py index a1ec3c0..71c9e8d 100644 --- 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)