Commit bf951972 authored by Laura LINDEPERG's avatar Laura LINDEPERG
Browse files

Tests for median slope extraction from pysheds library

parent 697f2be8
......@@ -6,14 +6,13 @@ Created on Wed Mar 17 15:35:20 2021
"""
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
from pysheds.grid import Grid
# **************************** Data *****************************
......@@ -111,36 +110,181 @@ FC_test = rxr.open_rasterio(FC_test_path).squeeze()
# *************************** Extract indicators ************************************
topographic_property_map = BDAlti_test
def extract_mean_hydraulic_property(hydraulic_property_map):
values, counts = np.unique(hydraulic_property_map[0], return_counts=True)
def extract_mean_elevation(topographic_property_map):
values, counts = np.unique(topographic_property_map, return_counts=True)
my_np_array = np.array([values, counts]).transpose()
ValAndCounts = pd.DataFrame(my_np_array, columns=['Values', 'Counts'])
ValAndCounts = ValAndCounts[ValAndCounts.loc[:, 'Values'] != -1.7e+308] # remove pixels outside the studied area (value -1.7e+308)
ValAndCounts = ValAndCounts[ValAndCounts.loc[:, 'Values'] != -99999] # remove pixels outside the studied area (value -99999)
weight = ValAndCounts.Counts/(ValAndCounts.Counts.sum())
average = (ValAndCounts.Values*weight).sum()
return average
# def extract_mean_hydraulic_property(hydraulic_property_map):
# values, counts = np.unique(hydraulic_property_map[0], return_counts=True)
# ValAndCounts = pd.DataFrame()
# ValAndCounts['Values'] = values
# ValAndCounts['Counts'] = counts
# 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
import numpy.ma as ma
mean_elevation = np.mean(ma.masked_values(topographic_property_map, -99999.))
def compute_KS_mean(hydraulic_property_map):
KS_mean_log10 = extract_mean_hydraulic_property(hydraulic_property_map)
KS_mean = pow(10,KS_mean_log10)
return KS_mean
# BDAlti 25m
BDAlti_map_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/BDAlti_25m/France_BDAlti_25m.tif'
BDAlti_map = rxr.open_rasterio(BDAlti_map_path).squeeze()
BDAlti_clipped = BDAlti_map.rio.clip_box(
minx=2.018e+05,
miny=6.772e+06,
maxx=2.342e+05,
maxy=6.809e+06,)
BDAlti_clipped.rio.to_raster('C:/Users/laura.lindeperg/Documents/DonneesLaura/BDAlti_25m/clipped.tif')
Alticlip = rxr.open_rasterio('C:/Users/laura.lindeperg/Documents/DonneesLaura/BDAlti_25m/clipped.tif').squeeze()
grid = Grid.from_raster('C:/Users/laura.lindeperg/Documents/DonneesLaura/BDAlti_25m/clipped.tif', data_name='dem')
grid.read_raster(data = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/BDAlti_25m/clipped.tif', data_name='dem', nodata=-99999.)
elevation = grid.dem
# Determine D8 flow directions from DEM
# ----------------------
# Fill depressions in DEM
grid.fill_depressions('dem', out_name='flooded_dem')
# Resolve flats in DEM
grid.resolve_flats('flooded_dem', out_name='inflated_dem')
# Specify directional mapping
dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
# Compute flow directions
# -------------------------------------
grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap)
fdir = grid.dir
grid.cell_slopes(fdir = 'dir', dem = 'dem', out_name = 'slopes', nodata_in=0)
slope = grid.slopes
grid = Grid.from_raster('C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/TOPO/BDAlti_25m/J4742020_BDAlti_25m.tif', data_name='dem')
grid.read_raster(data = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/TOPO/BDAlti_25m/J4742020_BDAlti_25m.tif', data_name='dem', nodata=-99999.)
import numpy.ma as ma
crs = BDAlti_test.rio.crs
BDAlti_newtest = np.array(BDAlti_test)
masked = ma.masked_values(BDAlti_test, -99999.)
mask = np.isfinite(masked)
new = mask.filled(False)
new.rio.to_raster('test.tif')
# x = ma.array(BDAlti_test, mask=masked)
y= masked.filled(np.nan)
'+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
grid.add_gridded_data(y, data_name = 'test', crs={'SRS' :'EPSG:2154'} )
# Determine D8 flow directions from DEM
# ----------------------
# Fill depressions in DEM
grid.fill_depressions('dem', out_name='flooded_dem', nodata_in=-99999.)
# Resolve flats in DEM
grid.resolve_flats('flooded_dem', out_name='inflated_dem', nodata_in=-99999.)
# Specify directional mapping
#N NE E SE S SW W NW
dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
KS_bv_test = compute_KS_mean(KS_test)
FC_bv_test = extract_mean_hydraulic_property(FC_test)
# Compute flow directions
# -------------------------------------
grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap, nodata_in = -99999.)
fdir = grid.dir
dirmap = (70, 100, 21, 20, 47, 82, 1, 2)
grid.flowdir(data='inflated_dem', out_name='dir_1', dirmap=dirmap)
fdir_1 = grid.dir_1
grid.cell_slopes(fdir = 'dir', dem = 'dem', out_name = 'slopes', nodata_in=0)
slope = grid.slopes
median_slope = np.median(slope[np.isfinite(slope)])
grid.cell_slopes(fdir = 'dir', dem = 'dem', dirmap = dirmap, out_name = 'slopes_dirmap')
slope_dirmap = grid.slopes_dirmap
median_slope_dirmap = np.median(slope_dirmap[np.isfinite(slope_dirmap)])
grid.cell_slopes(fdir = 'dir_1', dem = 'dem', out_name = 'slopes_1')
slope_1 = grid.slopes_1
median_slope_1 = np.median(slope_1[np.isfinite(slope_1)])
grid.flowdir(data='inflated_dem', out_name='dir')
grid.flowdir(data='inflated_dem_1', out_name='dir_1', nodata_in= -99999)
grid.view('dir')
grid.cell_slopes(fdir = 'dir', dem = 'BDAlti_raster', out_name = 'slopes')
grid.cell_slopes(fdir = 'dir_1', dem = 'BDAlti_raster', out_name = 'slopes_1', nodata_in = 0)
grid.cell_area()
array_area = grid.area
dem = grid.BDAlti_raster
slope = grid.slopes
slope_1 = grid.slopes_1
median_slope = np.median(slope[np.isfinite(slope)])
median_slope_1 = np.median(slope[~np.isnan(slope)])
from numpy import nanmedian, NaN
median_slope_2 = nanmedian(slope)
fdir_1 = grid.dir_1
#N NE E SE S SW W NW
dirmap = (70, 100, 10, 20, 44, 8, 16, 32)
# Access the flow direction grid
fdir = grid.dir
# Plot the flow direction grid
fig = plt.figure(figsize=(8,6))
plt.imshow(fdir, extent=grid.extent, cmap='viridis', zorder=1)
boundaries = ([0] + sorted(list(dirmap)))
plt.colorbar(boundaries= boundaries,
values=sorted(dirmap))
plt.title('Flow direction grid')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
# Plot the result
fig, ax = plt.subplots(figsize=(8,6))
plt.imshow(raster_slopes, cmap='bone', zorder=1,
extent=grid.extent)
plt.colorbar(boundaries = [0, 1], label='Slopes (m - cells)')
plt.title('Slopes')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
......@@ -194,7 +338,30 @@ plt.show()
def extract_mean_hydraulic_property(hydraulic_property_map):
values, counts = np.unique(hydraulic_property_map[0], return_counts=True)
my_np_array = np.array([values, counts]).transpose()
ValAndCounts = pd.DataFrame(my_np_array, columns=['Values', 'Counts'])
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
# def extract_mean_hydraulic_property(hydraulic_property_map):
# values, counts = np.unique(hydraulic_property_map[0], return_counts=True)
# ValAndCounts = pd.DataFrame()
# ValAndCounts['Values'] = values
# ValAndCounts['Counts'] = counts
# 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
def compute_KS_mean(hydraulic_property_map):
KS_mean_log10 = extract_mean_hydraulic_property(hydraulic_property_map)
KS_mean = pow(10,KS_mean_log10)
return KS_mean
......
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