diff --git a/GeologicProperties.py b/GeologicProperties.py index e1aaf14bd207ab8df636398c387af9ead640901a..4e74c977f45e18421d4b3f2ddd158f4847073f12 100644 --- a/GeologicProperties.py +++ b/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) diff --git a/TestData/W3315010_ZHp.tif b/TestData/W3315010_ZHp.tif new file mode 100644 index 0000000000000000000000000000000000000000..f6e33546f9dc584f81969bac3f38577f30f43672 Binary files /dev/null and b/TestData/W3315010_ZHp.tif differ diff --git a/read_ZHPdata.py b/read_ZHPdata.py new file mode 100644 index 0000000000000000000000000000000000000000..abe0939511790309a270bea4080e34f5298f0d89 --- /dev/null +++ b/read_ZHPdata.py @@ -0,0 +1,52 @@ +# -*- 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() + +