Merge branch 'master' of https://gitlab.irstea.fr/jean.marcais/fhysa

GeologicProperties.py

@@ -25,8 +25,14 @@ class GeologicProperties(object):
 
     def extract_average_age_geology(self, geol_foldername, catchment_contour_shp):
         geol_shp = gpd.read_file(geol_foldername)
-        geol_shp = gpd.clip(geol_shp, catchment_contour_shp)
-
+        
+        try:
+            gpd.clip(geol_shp, catchment_contour_shp)
+        except:
+            gpd.overlay(geol_shp, catchment_contour_shp)
+        
+    
+    
         geol_shp['area'] = geol_shp.geometry.area / 1e6
         geol_shp['prop'] = geol_shp.area / sum(geol_shp.area)
         age_min = sum(geol_shp.prop*geol_shp.AGE_SUP)

read_ZHPdata.py

+# -*- coding: utf-8 -*-
+"""
+Created on Wed Mar 17 15:35:20 2021
+
+@author: laura.lindeperg
+"""
+
+import matplotlib.pyplot as plt
+import seaborn as sns
+import numpy as np
+from shapely.geometry import mapping
+import rioxarray as rxr
+import xarray as xr
+import geopandas as gpd
+import pandas as pd
+
+
+# Carte des Zones Humides (raster)
+ZHp_map_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/ZonesHumidesPotentielles/mph-fr.tif'
+ZHp_map = rxr.open_rasterio(ZHp_map_path).squeeze()
+
+# Contour du bassin versant test (shapefile)
+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)
+
+# extraction du raster propre au contour du bassin versant et enregistrement du raster obtenu
+ZHp_clipped = ZHp_map.rio.clip(BV_W3315010.geometry.apply(mapping), BV_W3315010.crs)
+ZHp_clipped.rio.to_raster('C:/Users/laura.lindeperg/Documents/DonneesLaura/ZonesHumidesPotentielles/W3315010_ZHp.tif')
+
+
+# Raster test
+ZHp_test_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/ZonesHumidesPotentielles/W3315010_ZHp.tif'
+ZHp_test = rxr.open_rasterio(ZHp_test_path).squeeze()
+
+# ZHp_test.values.flatten() -> fonction "unique" s'en occupe si ce n'est pas 1D
+classes, counts = np.unique(ZHp_test, return_counts=True)
+
+classesO, countsO = np.unique(ZHp_map, return_counts=True)
+
+# Quelques tests
+
+ZHp_test.plot.hist()
+plt.show()
+
+f, ax = plt.subplots(figsize=(10, 4))
+ZHp_test.plot(ax=ax,
+                 cmap='Greys')
+ax.set(title="Final Clipped ZHp")
+ax.set_axis_off()
+plt.show()
+
+