keep working on ZHp ratio extraction and BDTopage

read_ZHPdata.py

@@ -21,6 +21,23 @@ import pandas as pd
 BV_W3315010_path = 'C:/Users/laura.lindeperg/Documents/INRAE_2021/CODE/fhysa/TestData/BH/W3315010_contour.shp'
 BV_W3315010 = gpd.read_file(BV_W3315010_path)
 
+bv_for_test_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMETRY/J4742020.shp'
+bv_for_test = gpd.read_file(bv_for_test_path)
+
+## BDTopage
+
+# Tronçons hydro
+troncons_hydro_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Topage/TronconHydrographique_FXX-shp/TronconHydrographique_FXX.shp'
+troncons_hydro = gpd.read_file(troncons_hydro_path)
+
+# Extraction d'un shp propre au bv test et enregistrement
+troncons_hydro_clipped = gpd.clip(troncons_hydro, bv_for_test)
+troncons_hydro_clipped.to_file('C:/Users/laura.lindeperg/Documents/DonneesLaura/Topage/TronconHydrographique_FXX-shp/bv_for_test_TronconHydro.shp')
+
+
+troncon_hydro_test_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Topage/TronconHydrographique_FXX-shp/bv_for_test_TronconHydro.shp'
+troncon_hydro_test = gpd.read_file(troncon_hydro_test_path)
+
 ## ZHp
 
 # Carte des Zones Humides (raster)
@@ -36,6 +53,19 @@ ZHp_test_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/ZonesHumidesPot
 ZHp_test = rxr.open_rasterio(ZHp_test_path).squeeze()  # masked=True pour masquer les noData
 
 
+# Watersheds
+shp_watersheds_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/complete_df_wrong_geometries.shp'
+shp_foldername = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMETRY/'
+
+shp_watersheds = gpd.read_file(shp_watersheds_path)
+
+# List of the stations'codes
+watershed_code = shp_watersheds.loc[:,'Code']
+# Get a sample of them for test
+code_for_test = watershed_code.loc[0:3]
+
+
+
 ## Hydraulic conductivity of topsoil
 
 # Raster file
@@ -54,8 +84,38 @@ KS_test = rxr.open_rasterio(KS_test_path).squeeze()
 
 # *************************** Extract indicators ************************************
 
-# ZHp_test.values.flatten() -> fonction "unique" s'en occupe si ce n'est pas 1D
-classes, counts = np.unique(ZHp_test, return_counts=True)
+
+# Tronçons hydro
+total_length_bv = troncon_hydro_test.length
+total_length_bv = total_length_bv.sum()
+
+# Keep class 3 only (intermittent)
+intermittent = troncon_hydro_test[troncon_hydro_test.loc[:, 'Persistanc'] == '3']
+intermittent_length = intermittent.length.sum()
+
+
+
+
+
+def compute_ZHp_ratio(ZHp_test):
+    # ZHp_test.values.flatten() -> fonction "unique" s'en occupe si ce n'est pas 1D
+    classes, counts = np.unique(ZHp_test, return_counts=True)
+    my_np_array = np.array([classes, counts]).transpose()
+    ZHp_ClassAndCounts = pd.DataFrame(my_np_array, columns=['Classes', 'Counts'])
+    ZHp_ClassAndCounts = ZHp_ClassAndCounts[ZHp_ClassAndCounts.loc[:, 'Classes'] != 65535]
+    total_classes = ZHp_ClassAndCounts.loc[:,'Counts'].sum()
+    ZHp_classes = ZHp_ClassAndCounts[ZHp_ClassAndCounts.loc[:, 'Classes'].isin([2,3])]
+    ZHp_classes = ZHp_classes.loc[:,'Counts'].sum()
+    ZHp_ratio = ZHp_classes/total_classes
+    return ZHp_ratio
+
+geomorpho_indicators = pd.DataFrame()
+for i in code_for_test:
+    watershed_contour_i = gpd.read_file(shp_foldername+i+'.shp')
+    ZHp_within_bv = ZHp_map.rio.clip(watershed_contour_i.geometry.apply(mapping), watershed_contour_i.crs)
+    ratio = compute_ZHp_ratio(ZHp_within_bv)
+    geomorpho_indicators = geomorpho_indicators.append({'Code':i, 'ZHp_ratio': ratio}, ignore_index=True)
+
 
 classesO, countsO = np.unique(ZHp_map, return_counts=True)