hechms_curve_numbers.txt

import geopandas as gpd
import pandas as pd
import os

SOIL_SHAPEFILE_PATH = r"<UPDATE>"
LAND_USE_SHAPEFILE_PATH = r"<UPDATE>"
SUBBASINS_SHAPEFILE_PATH = r"<UPDATE>"
OUTPUT_SHAPEFILE_PATH = r"<UPDATE>"
OUTPUT_SPREADSHEET_PATH = r"<UPDATE>"
HSG_FIELD_NAME = "<UPDATE>"
LAND_USE_FIELD_NAME = "<UPDATE>"
SOIL_TYPE_FIELD_NAME = "<UPDATE>"
PERCENT_IMPERVIOUS_FIELD_NAME = "<UPDATE>"
CN_A_FIELD_NAME = "<UPDATE>"
CN_B_FIELD_NAME = "<UPDATE>"
CN_C_FIELD_NAME = "<UPDATE>"
CN_D_FIELD_NAME = "<UPDATE>"
SUBBASINS_ID_FIELD_NAME = "<UPDATE>"

# Function to intersect geometries and combine attributes from two dataframes
def intersect_geometries(df1, df2):
    intersections = []
    for index1, row1 in df1.iterrows():
        for index2, row2 in df2.iterrows():
            if row1['geometry'].intersects(row2['geometry']):
                intersected_geom = row1['geometry'].intersection(row2['geometry'])
                combined_attributes = row1.to_dict()
                combined_attributes.update(row2.to_dict())
                combined_attributes['geometry'] = intersected_geom
                intersections.append(combined_attributes)
    return gpd.GeoDataFrame(intersections)

# Calculate area-weighted CN for each subbasin
def weighted_cn_calculation(group):
    weighted_sum = (group['intersected_area'] * group['selected_CN']).sum()
    total_area = group['intersected_area'].sum()
    return weighted_sum / total_area

# Load the shapefiles
subbasins_gdf = gpd.read_file(SUBBASINS_SHAPEFILE_PATH)
soils_gdf = gpd.read_file(SOILS_SHAPEFILE_PATH)
land_use_gdf = gpd.read_file(LAND_USE_PATH)

# Ensure CRS compatibility and convert if necessary
if subbasins_gdf.crs != soils_gdf.crs:
    subbasins_gdf = subbasins_gdf.to_crs(soils_gdf.crs)

if land_use_gdf.crs != soils_gdf.crs:
    land_use_gdf = land_use_gdf.to_crs(soils_gdf.crs)

# Intersect subbasins with soils
subbasin_soil_intersection = intersect_geometries(subbasins_gdf, soils_gdf)

# Intersect the above result with land use
final_intersection = intersect_geometries(subbasin_soil_intersection, land_use_gdf)

# Calculate the area of each intersection
final_intersection['intersected_area'] = final_intersection.geometry.area

final_intersection['selected_CN'] = final_intersection.apply(select_cn, axis=1)

subbasin_cn = final_intersection.groupby(subbasin_name_field).apply(weighted_cn_calculation)

# Convert to a DataFrame for display
subbasin_cn_df = subbasin_cn.reset_index()
subbasin_cn_df.columns = ['Subbasin', 'Area-Weighted Average CN']

print(subbasin_cn_df)