extract_Geol.py 2.23 KiB
# -*- coding: utf-8 -*-
"""
Created on Tue Apr  6 15:26:11 2021

@author: laura.lindeperg
"""

import rioxarray as rxr
import geopandas as gpd
import pandas as pd


# **************************** Data *****************************

# 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)

## Geologic data

BDLisa_filepath = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOL/BDLisa/'

BRGM_filepath = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOL/BRGM/'

# KS - Saturated hydraulic conductivity of topsoil - log10(cm/day)
KS_foldername = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOL/JRC-ESDAC/KS/'

# SATURATED - Water retention of topsoil: saturated water content - cm3/cm3 no units
SAT_foldername = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOL/JRC-ESDAC/SAT/'

# FC - Water retention of topsoil:  water content at  field capacity - cm3/cm3 no units
FC_foldername = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOL/JRC-ESDAC/FC/'



# 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]



# *************************** Extract indicators ************************************


from GeologicProperties import GeologicProperties
geol_indicators = pd.DataFrame()
# for i in code_for_test:
for i in watershed_code:
    watershed_contour_i = gpd.read_file(shp_foldername+i+'.shp')
    BDLisa_i = gpd.read_file(BDLisa_filepath+i+'_BDLisa.shp')
    ks_i = rxr.open_rasterio(KS_foldername+i+'_KS.tif')
    sat_i = rxr.open_rasterio(SAT_foldername+i+'_SAT.tif')
    fc_i = rxr.open_rasterio(FC_foldername+i+'_FC.tif')
    
    geol_i = GeologicProperties(i)
    geol_i.extract_main_geology(BDLisa_i)
    geol_i.extract_average_age_geology(BRGM_filepath+i+'_BRGM.shp')
    geol_i.compute_hydraulic_properties(ks_i, sat_i, fc_i)

    geol_indicators = geol_indicators.append(geol_i.__dict__, ignore_index=True)

# And save it
# geol_indicators.to_csv('616_stations_geol_df.csv', index=False)