Commit 56cbcad1 authored by Laura LINDEPERG's avatar Laura LINDEPERG
Browse files

First tests with whitebox tools

parent 056e2081
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 21 14:19:12 2021
@author: laura.lindeperg
"""
import os
import whitebox
import rioxarray as rxr
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
wbt = whitebox.WhiteboxTools()
# print(wbt.version())
# print(wbt.help())
# print(wbt.list_tools())
# # Set the working directory, i.e. the folder containing the data,
# # to the 'data' sub-directory.
# wbt.set_working_dir('C:/Users/laura.lindeperg/Documents/INRAE_2021/CODE/fhysa/Ressources')
# wbt.verbose = False
# path = 'C:/Users/laura.lindeperg/Documents/INRAE_2021/CODE/fhysa'
# os.chdir(path)
# data_dir = './Ressources/'
# print(os.listdir(data_dir))
# # set whitebox working directory
# wbt.set_working_dir(data_dir)
# 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)
## TOPO
# # 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_filepath = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/TOPO/BDAlti_25m/'
# *************************** Create watersheds' rasters ************************************
# 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:1]
# Set the working directory, i.e. the folder containing the data,
# to the 'data' sub-directory.
# wbt.set_working_dir('C:/Users/laura.lindeperg/Documents/DonneesLaura/BDAlti_25m')
wbt.verbose = False
## TOPOGRAPHY
# BDAlti 25m
# for i in code_for_test:
for i in watershed_code:
raster_input_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/TOPO/BDAlti_25m/'+i+'_BDAlti_25m.tif'
raster_output_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/TOPO/Slope/'+i+"_slope_25m.tif"
wbt.slope(raster_input_path, raster_output_path, units="percent")
for i in code_for_test:
topage_shp = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/BDTopage/'+i+'_BDTopage.shp'
topage_raster = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/BDTopage/Raster/BDTopage_temp.tif'
model_raster = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/BDAlti_25m/'+i+'_BDAlti_25m.tif'
# wbt.vector_lines_to_raster(topage_shp, topage_raster, field = -99999., cell_size=1.0)
wbt.vector_lines_to_raster(topage_shp, topage_raster, field = -99999., base = model_raster)
dem_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/BDAlti_25m/'+i+'_BDAlti_25m.tif'
breached_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/BDAlti_25m/Breached/'+i+'_breached.tif'
wbt.breach_depressions(dem_path, breached_path)
d8pointer_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/BDAlti_25m/D8pointer/'+i+'_d8_pointer.tif'
wbt.d8_pointer(breached_path, d8pointer_path)
strahler_stream_order_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/BDAlti_25m/StrahlerStreamOrder/'+i+'_strahler_stream_order.tif'
wbt.strahler_stream_order(d8pointer_path, topage_raster, strahler_stream_order_path)
topage_raster_test = rxr.open_rasterio('C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/BDTopage/Raster/A1050310_BDTopage.tif', masked=True).squeeze()
array_box = np.array(topage_raster_test)
topage_raster_test.plot.hist()
plt.show()
f, ax = plt.subplots(figsize=(10, 4))
topage_raster_test.plot(ax=ax,
cmap='Greys')
ax.set(title="Topage raster")
ax.set_axis_off()
plt.show()
slope_test = rxr.open_rasterio('C:/Users/laura.lindeperg/Documents/INRAE_2021/CODE/fhysa/Ressources/slope.tif').squeeze()
array = np.array(slope_test)
import numpy.ma as ma
median_slope = np.median(ma.masked_values(slope_test, -99999.))
# raster_BDAlti = BDAlti_map.rio.clip(shp_contour_i.geometry.apply(mapping), shp_contour_i.crs)
# raster_BDAlti.rio.to_raster('C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/TOPO/BDAlti_25m/'+i+'_BDAlti_25m.tif')
wbt.slope("clipped.tif", "slope_from_clipped.tif", units="percent")
slope_test = rxr.open_rasterio('C:/Users/laura.lindeperg/Documents/DonneesLaura/BDAlti_25m/slope_from_clipped.tif').squeeze()
array_box = np.array(slope_test)
bv_for_test_path = 'C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMETRY/J4742020.shp'
bv_for_test = gpd.read_file(bv_for_test_path)
from shapely.geometry import mapping
slope_clipped = slope_test.rio.write_crs(bv_for_test.crs)
slope_clipped = slope_clipped.rio.clip(bv_for_test.geometry.apply(mapping), bv_for_test.crs)
array_clipped = np.array(slope_clipped)
import numpy.ma as ma
median_slope = np.median(ma.masked_values(slope_clipped, -99999.))
slope_test = rxr.open_rasterio('C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/Slope/A2332110_slope_25m.tif').squeeze()
array = np.array(slope_test)
elevation_test = rxr.open_rasterio('C:/Users/laura.lindeperg/Documents/DonneesLaura/Watersheds/GEOMORPHO/BDAlti_25m/A2332110_BDAlti_25m.tif').squeeze()
array_elevation = np.array(elevation_test)
wbt.breach_depressions("smoothed.tif", "breached.tif")
wbt.d_inf_flow_accumulation("breached.tif", "flow_accum.tif")
\ No newline at end of file
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