<jupyter_start><jupyter_text>找出Susquehanna区域中的GAGES-II和CAMELS流域<jupyter_code>import geopandas as gpd
import pandas as pd
import os

fp = os.path.join("data", "3557basins_ID_NSE_DOR.csv")
# Read file using gpd.read_file()
data = pd.read_csv(fp, dtype={"GAUGE ID": str})
data.head()
data.shape<jupyter_output><empty_output><jupyter_text>下载流域边界shapefile: https://water.usgs.gov/GIS/dsdl/boundaries_shapefiles_by_aggeco.zip解压后文件夹放入 data 文件夹。下面先试试一个shpfile中筛选：<jupyter_code>shpfile_dir = os.path.join("data", "boundaries-shapefiles-by-aggeco")
shpfile1 = os.path.join(shpfile_dir, "bas_nonref_CntlPlains.shp")
shp1 = gpd.read_file(shpfile1)
shp1.head()
shp1_chosen = shp1.loc[shp1["GAGE_ID"].isin(data["GAUGE ID"].values)]
shp1_chosen<jupyter_output><empty_output><jupyter_text>下面看看两个：<jupyter_code>shpfile2 = os.path.join(shpfile_dir, "bas_nonref_EastHghlnds.shp")
shp2 = gpd.read_file(shpfile2)
shp2_chosen = shp2.loc[shp2["GAGE_ID"].isin(data["GAUGE ID"].values)]
chosen12 = [shp1_chosen, shp2_chosen]
rdf = gpd.GeoDataFrame( pd.concat( chosen12, ignore_index=True) )
rdf<jupyter_output><empty_output><jupyter_text>接下来就可以尝试多个了：<jupyter_code>shpfiles = ["bas_nonref_CntlPlains.shp", "bas_nonref_EastHghlnds.shp", "bas_nonref_MxWdShld.shp", "bas_nonref_NorthEast.shp", "bas_nonref_SECstPlain.shp", 
            "bas_nonref_SEPlains.shp", "bas_nonref_WestMnts.shp", "bas_nonref_WestPlains.shp", "bas_nonref_WestXeric.shp", "bas_ref_all.shp"]
chosen_lst = []
for shpfile_name in shpfiles:
    shpfile = os.path.join(shpfile_dir,shpfile_name)
    shp = gpd.read_file(shpfile)
    shp_chosen = shp.loc[shp["GAGE_ID"].isin(data["GAUGE ID"].values)]
    chosen_lst.append(shp_chosen)
merge_chosen = gpd.GeoDataFrame( pd.concat( chosen_lst, ignore_index=True) )
merge_chosen
shp2.crs<jupyter_output><empty_output><jupyter_text>最后按id排序，并将结果输出到一个新的shpfile中。<jupyter_code>sorted_merge_chosen = merge_chosen.sort_values(by=['GAGE_ID'], ignore_index=True)
sorted_merge_chosen
# Set the GeoDataFrame's coordinate system
sorted_merge_chosen.crs = shp2.crs
if not os.path.isdir("results"):
    os.mkdir("results")
sorted_merge_chosen.to_file(os.path.join("results", "3557basins.shp"))<jupyter_output><empty_output><jupyter_text>QGIS看了一下3557basins和Susquehanna的相交情况，就选择几个有相交的流域就行，流域90%以上在Susquehanna中的都算Susquehanna里的。<jupyter_code>shpfile_susquehanna = os.path.join("data", "HUC8_Susquehanna.shp")
shp_susquehanna = gpd.read_file(shpfile_susquehanna)
shp_susquehanna.head()
shp_susquehanna['new_column'] = 0
susquehanna_union = shp_susquehanna.dissolve(by='new_column')
susquehanna_union
susquehanna_union.to_file(os.path.join("results", "susquehanna_union.shp"))
print(susquehanna_union.crs)
print(sorted_merge_chosen.crs)
susquehanna_union = susquehanna_union.to_crs(sorted_merge_chosen.crs)
susquehanna_union.crs == sorted_merge_chosen.crs
basin3557_polygon_lst= [sorted_merge_chosen.geometry[i] for i in range(sorted_merge_chosen.shape[0])]<jupyter_output><empty_output><jupyter_text>找到相交的流域<jupyter_code>intersects_indices = [i for i in range(len(basin3557_polygon_lst)) if basin3557_polygon_lst[i].intersects(susquehanna_union.geometry[0])]
intersects_indices
intersect_gdf = sorted_merge_chosen.iloc[intersects_indices]
intersect_gdf
intersect_gdf["GAGE_ID"].values
intersect_gdf.to_file(os.path.join("results", "susquehanna_gages_intersect.shp"))<jupyter_output><empty_output><jupyter_text>CAMELS流域与Susquehanna的相交，从这里下载CAMELS数据：https://ral.ucar.edu/solutions/products/camels<jupyter_code>shpfile_camels = os.path.join("data", "HCDN_nhru_final_671.shp")
camels_gdf = gpd.read_file(shpfile_camels)
camels_gdf
camels_gdf.crs
susquehanna4camels = susquehanna_union.to_crs(camels_gdf.crs)
susquehanna4camels.crs == camels_gdf.crs
camels_susquehanna_intersects_indices = [i for i in range(camels_gdf.shape[0]) if camels_gdf.geometry[i].intersects(susquehanna4camels.geometry[0])]
camels_susquehanna_intersects_indices
camels_susquehanna_intersect_gdf = camels_gdf.iloc[camels_susquehanna_intersects_indices]
camels_susquehanna_intersect_gdf
camels_susquehanna_intersect_gdf.to_file(os.path.join("results", "susquehanna_camels_intersect.shp"))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12159884420900553, 'fraction_numerical': 0.014447387430772935, 'mean_word_length': 10.734463276836157, 'pattern_counts': {'":': 1, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '211599', 'license': 'MIT', 'n_tokens_mistral': 1959, 'n_tokens_neox': 1900, 'n_words': 163, 'path': '1-find-susquehanna-intersect.ipynb', 'repo_name': 'OuyangWenyu/susquehanna', 'size': '4153'}