Commit 9a5d35a6 authored by Marcais Jean's avatar Marcais Jean
Browse files
parents 5c9a8f89 6a76a76f
......@@ -4,7 +4,7 @@ import pandas as pd
class GeologicProperties(object):
def __init__(self, code=-1, maingeol_id=-1, maingeol_description=-1, age=-1, proportion=-1, maingeol_age=-1, maingeol_age_proportion=-1, ks=-1, fc=-1):
def __init__(self, code=-1, maingeol_id=-1, maingeol_description=-1, age=-1, proportion=-1, maingeol_age=-1, maingeol_age_proportion=-1, ks=-1, theta_d=-1, k_glhymps = -1):
self.code = code
self.maingeol_id = maingeol_id
self.maingeol_description = maingeol_description
......@@ -12,8 +12,9 @@ class GeologicProperties(object):
self.age = age # [min, mean, max] age averaged from the 1/1000000e brgm map
self.maingeol_age = maingeol_age
self.maingeol_age_proportion = maingeol_age_proportion
self.ks = ks
self.fc = fc
self.ks = ks # [cm/day]
self.theta_d = theta_d
self.k_glhymps = k_glhymps # [cm/day]
def extract_main_geology_from_filename(self, geol_foldername, catchment_contour_shp):
geol_shp = gpd.read_file(geol_foldername)
......@@ -55,17 +56,36 @@ class GeologicProperties(object):
ValAndCounts = ValAndCounts[ValAndCounts.loc[:, 'Values'] != -1.7e+308] # remove pixels outside the studied area (value -1.7e+308)
weight = ValAndCounts.Counts/(ValAndCounts.Counts.sum())
average = (ValAndCounts.Values*weight).sum()
return average
return average
def compute_drainage_porosity(self, saturated_water_content_map, field_capacity_map):
theta_s = self.extract_mean_hydraulic_property(saturated_water_content_map)
theta_fc = self.extract_mean_hydraulic_property(field_capacity_map)
self.theta_d = theta_s-theta_fc
def compute_KS_mean(self, hydraulic_property_map):
KS_mean_log10 = self.extract_mean_hydraulic_property(hydraulic_property_map)
self.ks = pow(10,KS_mean_log10)
def compute_hydraulic_properties(self, ks_map, fc_map):
def compute_hydraulic_properties(self, ks_map, sat_map, fc_map):
self.compute_KS_mean(ks_map)
self.fc = self.extract_mean_hydraulic_property(fc_map)
self.compute_drainage_porosity(sat_map, fc_map)
def extract_K_GLHYMPS(self, glhymps_file):
where_gum_k_0 = glhymps_file['GUM_K'].eq(0)
log_k = glhymps_file['logK_Ferr_']
log_k = log_k.where(~where_gum_k_0, np.nan)
k = pow(10,log_k/100) # k [m2] permeability
K = k*999.97*9.8/(1e-3)
K = np.nanmean(K)
# k = np.nanmean(k)
# K = k*999.97*9.8/(1e-3) # K [m/s] hydraulic conductivity
K = K*1e2*24*3600
self.k_glhymps = K
## From original BDLisa and BRGM files using clip function
......
......@@ -8,15 +8,19 @@ Created on Fri Apr 2 14:38:24 2021
import numpy as np
import pandas as pd
import rioxarray as rxr
class GeomorphologicProperties(object):
def __init__(self, code=-1, ZHp_ratio=-1, order_1_ratio=-1, intermittent_ratio=-1, drainage_density=-1):
def __init__(self, code=-1, ZHp_ratio=-1, intermittent_ratio=-1, drainage_density=-1, mean_elevation = -1, median_slope=-1, mean_DTB = -1):
self.code = code
self.ZHp_ratio = ZHp_ratio
self.order_1_ratio = order_1_ratio
# self.order_1_ratio = order_1_ratio
self.intermittent_ratio = intermittent_ratio
self.drainage_density = drainage_density
self.drainage_density = drainage_density # [m-1]
self.mean_elevation = mean_elevation # [m]
self.median_slope = median_slope # [%]
self.mean_DTB = mean_DTB # [cm]
......@@ -72,13 +76,74 @@ class GeomorphologicProperties(object):
drainage_density = float(total_length/surface_area)
self.drainage_density = drainage_density
def extract_mean_value_from_raster(self, raster_path):
raster = rxr.open_rasterio(raster_path, masked=True).squeeze()
mean_value = np.nanmean(raster) # if masked = True when reading raster
return mean_value
def extract_median_value_from_raster(raster_path):
raster = rxr.open_rasterio(raster_path, masked=True).squeeze()
median_value = np.nanmedian(raster) # if masked = True when reading raster
return median_value
# def extract_mean_elevation(self, raster_path):
# self.mean_elevation = self.extract_mean_value_from_raster(raster_path)
# def extract_median_slope(self, raster_path):
# self.median_slope = self.extract_median_value_from_raster(raster_path)
# def extract_mean_DTB(self, raster_path):
# self.mean_DTB = self.extract_mean_value_from_raster(raster_path)
def extract_mean_elevation(self, raster_path):
raster = rxr.open_rasterio(raster_path, masked=True).squeeze()
mean_value = np.nanmean(raster) # if masked = True when reading raster
self.mean_elevation = mean_value
def extract_median_slope(self, raster_path):
raster = rxr.open_rasterio(raster_path, masked=True).squeeze()
median_value = np.nanmedian(raster) # if masked = True when reading raster
self.median_slope = median_value
def extract_mean_DTB(self, raster_path):
raster = rxr.open_rasterio(raster_path, masked=True).squeeze()
mean_value = np.nanmean(raster) # if masked = True when reading raster
self.mean_DTB = mean_value
# def extract_mean_elevation(self, topographic_property_map):
# import numpy.ma as ma
# mean_elevation = np.mean(ma.masked_values(topographic_property_map, topographic_property_map._FillValue))
# self.mean_elevation = mean_elevation
# def extract_median_slope(self, topographic_property_map):
# import numpy.ma as ma
# slope = ma.masked_values(topographic_property_map, topographic_property_map._FillValue)
# slope = slope.filled(np.nan)
# median_slope = np.nanmedian(slope)
# self.median_slope = median_slope
# def extract_mean_DTB(self, DTB_map):
# import numpy.ma as ma
# mean_DTB = np.mean(ma.masked_values(DTB_map, DTB_map._FillValue)) # if masked = False when reading raster
# return mean_DTB
# def extract_mean_DTB(self, DTB_map):
# mean_DTB = np.nanmean(DTB_map) # if masked = True when reading raster
# return mean_DTB
def extract_geomorphologic_properties(self, contour, ZHp, hydro_sections_Topage, hydro_sections_Carthage):
def extract_geomorphologic_properties(self, contour, ZHp, hydro_sections_Topage, alti_path, slope_path, DTB_path):
self.compute_ZHp_ratio(ZHp)
self.compute_order_1_ratio(hydro_sections_Carthage)
# self.compute_order_1_ratio(hydro_sections_Carthage)
self.compute_intermittent_ratio(hydro_sections_Topage)
self.compute_drainage_density(contour, hydro_sections_Topage)
self.extract_mean_elevation(alti_path)
self.extract_median_slope(slope_path)
self.extract_mean_DTB(DTB_path)
......
......@@ -4,10 +4,13 @@ import pandas as pd
class HydroClimaticFluxes(object):
def __init__(self, code, discharge_df=-1, safran_df=-1):
def __init__(self, code, discharge_df=-1, safran_df=-1, discharge_timeseries_quality = -1, deleted_Hydro_Year = -1, entire_q_timeseries = -1):
self.code = code # identifier of the site (e.g. gauge number)
self.discharge_timeseries = discharge_df # dataframe containing discharge time series
self.safran_timeseries = safran_df # dataframe containing P, ET0 and Tair reanalysis over 1958 to present
self.discharge_timeseries_quality = discharge_timeseries_quality
self.deleted_Hydro_Year = deleted_Hydro_Year
self.entire_q_timeseries = entire_q_timeseries
def extract_discharge_timeseries(self, filepath, type_data='txt_file'):
if type_data == 'BanqueHydro':
......@@ -19,7 +22,12 @@ class HydroClimaticFluxes(object):
def keep_clean_HydroYear(self):
import read_BankHydro as readBH
self.discharge_timeseries = readBH.delete_HydroYear_withNan(self.discharge_timeseries)
q_timeseries, quality, removed_years, entire_timeseries = readBH.delete_HydroYear_withNan(self.discharge_timeseries)
self.discharge_timeseries = q_timeseries
self.discharge_timeseries_quality = quality
self.deleted_Hydro_Year = removed_years
self.entire_q_timeseries = entire_timeseries
def extract_safran_variable(self, foldername, quantity_to_retrieve='Ptot'):
import netCDF4 as nc
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
# -*- coding: utf-8 -*-
"""
Created on Fri May 21 15:45:10 2021
@author: laura.lindeperg
"""
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
# **************************** Data *****************************
# Load hydrological signatures df
df_hydro_sig_path = 'C:/Users/laura.lindeperg/Documents/INRAE_2021/616_stations_hydrosig_dim_df.csv'
df_hydro_sig = pd.read_csv(df_hydro_sig_path)
# Watersheds
shp_watersheds_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/616_Catchments.shp'
shp_watersheds = gpd.read_file(shp_watersheds_path)
rrse_241_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/STATIONS/Stations241_RRSE.csv'
rrse_241 = pd.read_csv(rrse_241_path)
rrse_241 = rrse_241.rename(columns = {'levels.RRSE.':'Code'})
rrse_241 = rrse_241.loc[:, 'Code']
# Geol
df_geol_path = 'C:/Users/laura.lindeperg/Documents/INRAE_2021/CODE/fhysa/616_stations_geol_df.csv'
df_geol = pd.read_csv(df_geol_path)
# Geomorpho
df_geomorpho_path = 'C:/Users/laura.lindeperg/Documents/INRAE_2021/CODE/fhysa/616_stations_geomorpho_df.csv'
df_geomorpho = pd.read_csv(df_geomorpho_path)
# Border catchments
border_catchments = ['A3792010', 'A3902010', 'B4224310', 'B4601010', 'B4631010', 'D0156510', 'O1900010', 'O2620010', 'O6140010', 'U2512010', 'U2542010', 'U2722010','V1000010', 'V1020010']
# Catchments size filter
size = shp_watersheds.loc[shp_watersheds.loc[:, 'S_km2'] < 5000]
# Main geol proportion filter
proportion_geol = df_geol.loc[df_geol.loc[:, 'maingeol_proportion'] > 0.70]
# Weird hydrographs filter
weird_catchments = pd.DataFrame(np.array(['H6201010', 'L4010710', 'Y1232010']), columns=['Code'])
weird_catchments=weird_catchments.loc[:,'Code']
geol = df_geol.loc[:, ['code', 'maingeol_description']]
my_data = pd.merge(df_hydro_sig, geol, on='code')
elevation = df_geomorpho.loc[:, ['code', 'mean_elevation']]
my_data = pd.merge(df_hydro_sig, elevation, on='code')
# Filtering workflow
filtered_watersheds = my_data.loc[df_hydro_sig.loc[:, 'code'].isin(rrse_241)==True]
filtered_watersheds = filtered_watersheds.loc[filtered_watersheds.loc[:, 'code'].isin(border_catchments)==False]
filtered_watersheds = filtered_watersheds.loc[filtered_watersheds.loc[:, 'code'].isin(size.loc[:,'Code'])==True]
filtered_watersheds = filtered_watersheds.loc[filtered_watersheds.loc[:, 'code'].isin(proportion_geol.loc[:,'code'])==True]
filtered_watersheds = filtered_watersheds.loc[filtered_watersheds.loc[:, 'code'].isin(weird_catchments)==False]
def budyko_equation(x):
# x=np.array(x)
y = np.sqrt(x*np.tanh(1/x)*(1-np.exp(-x)))
return y
x = np.linspace(0, 1.2, 100)
ET0_P = np.array(1/filtered_watersheds.aridity_ratio)
E_P = np.array(1-filtered_watersheds.runoff_ratio)
ET0_P = np.array(1/df_hydro_sig.aridity_ratio)
E_P = np.array(1-df_hydro_sig.runoff_ratio)
fig, ax = plt.subplots(figsize=(8,5)) # Create a figure containing a single axes.
ax.plot([0, 1, 1.2], [0, 1, 1]) # Plot some data on the axes.
ax.plot(x, budyko_equation(x))
sns.scatterplot(ax=ax, x=1/filtered_watersheds.aridity_ratio, y=1-filtered_watersheds.runoff_ratio, hue = 'mean_elevation', palette = 'viridis', data = filtered_watersheds).set(xlabel = 'ET0/P = 1/aridity_ratio', ylabel = 'E/P = 1-runoff_ratio')
plt.legend(title = 'mean elevation [m]')
for i in filtered_watersheds.code:
watershed_i = filtered_watersheds.loc[filtered_watersheds.code == i]
if float(1-watershed_i.runoff_ratio) > float(1/watershed_i.aridity_ratio):
plt.text(1/watershed_i.aridity_ratio, 1-watershed_i.runoff_ratio, i, horizontalalignment='left', size='medium', color='black')
print(i)
# ax.scatter(ET0_P, E_P)
# plt.legend(loc='center left', bbox_to_anchor=(1.04, 0.5))
sns.relplot(x=1/filtered_watersheds.aridity_ratio, y=1-filtered_watersheds.runoff_ratio)
sns.relplot(x=1/df_hydro_sig.aridity_ratio, y=1-df_hydro_sig.runoff_ratio)
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -113,8 +113,13 @@ def delete_HydroYear_withNan(df_obs, NaNvalues_number=20):
df_obs['HydroYear'] = pd.DatetimeIndex(df_obs['HydroDate']).year
# get yearly average with number of NaN values per year
df_obs_year = df_obs.resample('AS', on='HydroDate').sum()
Hydro_Year_to_Delete = df_obs_year[df_obs_year.NaN > NaNvalues_number ].index.year.values
full_discharge_timeseries = pd.DataFrame(df_obs)
Hydro_Year_to_Delete = df_obs_year[df_obs_year.NaN > NaNvalues_number].index.year.values
Hydro_Year_Information = df_obs_year['NaN']
Hydro_Year_Information.index = Hydro_Year_Information.index.year
if not Hydro_Year_to_Delete.any():
print('Not too much Nan Values per year. So no Hydrological Year have been deleted from datasets \n')
......@@ -125,4 +130,4 @@ def delete_HydroYear_withNan(df_obs, NaNvalues_number=20):
Q_mean_annual = df_obs['Q'].mean()
df_obs['Q'] = df_obs['Q'].fillna(Q_mean_annual)
return df_obs
\ No newline at end of file
return df_obs, Hydro_Year_Information, Hydro_Year_to_Delete, full_discharge_timeseries
\ No newline at end of file
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment