Commit 7bd02e73 authored by Laura LINDEPERG's avatar Laura LINDEPERG
Browse files

Mean DTB investigation

parent 72985b91
......@@ -8,17 +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, intermittent_ratio=-1, drainage_density=-1, mean_elevation = -1, median_slope=-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.intermittent_ratio = intermittent_ratio
self.drainage_density = drainage_density
self.mean_elevation = mean_elevation
self.median_slope = median_slope
self.drainage_density = drainage_density # [m-1]
self.mean_elevation = mean_elevation # [m]
self.median_slope = median_slope # [%]
self.mean_DTB = mean_DTB # [cm]
......@@ -74,26 +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_mean_elevation(self, topographic_property_map):
import numpy.ma as ma
mean_elevation = np.mean(ma.masked_values(topographic_property_map, -99999.))
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, -99999.)
slope = slope.filled(np.nan)
median_slope = np.nanmedian(slope)
self.median_slope = median_slope
# 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):
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_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)
......
......@@ -36,6 +36,9 @@ bdalti_foldername = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/
## Slope
slope_foldername = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/Slope/'
## DTB
DTB_foldername = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/DTB/'
# List of the stations'codes
watershed_code = shp_watersheds.loc[:,'Code']
......@@ -54,15 +57,18 @@ for i in watershed_code:
watershed_contour_i = gpd.read_file(shp_foldername+i+'.shp')
ZHp_i = rxr.open_rasterio(ZHp_foldername+i+'_ZHp.tif')
hydro_sections_Topage_i = gpd.read_file(topage_foldername+i+'_BDTopage.shp')
# hydro_sections_Carthage_i = gpd.read_file(carthage_foldername+i+'_BDCarthage.shp')
bdalti_i = rxr.open_rasterio(bdalti_foldername+i+'_BDAlti_25m.tif')
slope_i = rxr.open_rasterio(slope_foldername+i+'_slope_25m.tif')
# hydro_sections_Carthage_i = gpd.read_file(carthage_foldername+i+'_BDCarthage.shp')
# bdalti_i = rxr.open_rasterio(bdalti_foldername+i+'_BDAlti_25m.tif')
# slope_i = rxr.open_rasterio(slope_foldername+i+'_slope_25m.tif')
bdalti_path_i = bdalti_foldername+i+'_BDAlti_25m.tif'
slope_path_i = slope_foldername+i+'_slope_25m.tif'
DTB_path_i = DTB_foldername+i+'_DTB_250m.tif'
geomorpho_i = GeomorphologicProperties(i)
# geomorpho_i.compute_ZHp_ratio(ZHp_i)
geomorpho_i.extract_geomorphologic_properties(watershed_contour_i, ZHp_i, hydro_sections_Topage_i)
geomorpho_i.extract_mean_elevation(bdalti_i)
geomorpho_i.extract_median_slope(slope_i)
geomorpho_i.extract_geomorphologic_properties(watershed_contour_i, ZHp_i, hydro_sections_Topage_i, bdalti_path_i, slope_path_i, DTB_path_i)
# geomorpho_i.extract_mean_elevation(bdalti_i)
# geomorpho_i.extract_median_slope(slope_i)
geomorpho_indicators = geomorpho_indicators.append(geomorpho_i.__dict__, ignore_index=True)
......@@ -74,3 +80,18 @@ for i in watershed_code:
# ************************** Merge with extra geomorpho data **************************
delta_v_path = 'DeltaV_qualitycheck.csv'
delta_v = pd.read_csv(delta_v_path)
delta_v = delta_v[['Code', 'DeltaV', 'alpha']]
delta_v = delta_v.rename(columns = {'Code':'code'})
geomorpho_indicators = pd.read_csv('616_stations_geomorpho_df.csv')
geomorpho_indicators = pd.merge(geomorpho_indicators, delta_v, on = 'code')
# Save
geomorpho_indicators.to_csv('616_stations_8geomorpho_df.csv', index=False)
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