In [1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from pysheds.grid import Grid
from shapely import geometry, ops
import fiona
import geopandas as gpd
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
Data from USGS hydrosheds project: https://hydrosheds.cr.usgs.gov/datadownload.php
In [2]:
grid = Grid.from_raster('../data/n30w100_con', data_name='dem')
grid.read_raster('../data/n30w100_dir', data_name='dir')
In [3]:
#N NE E SE S SW W NW
dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
In [4]:
# Specify pour point
x, y = -97.294167, 32.73750
# Delineate the catchment
grid.catchment(data='dir', x=x, y=y, dirmap=dirmap, out_name='catch',
recursionlimit=15000, xytype='label', nodata_out=0)
# Clip to catchment
grid.clip_to('catch')
In [5]:
# Calling grid.polygonize without arguments will default to the catchment mask
shapes = grid.polygonize()
# Is equivalent to...
shapes = grid.polygonize(grid.mask.astype(np.uint8))
In [6]:
schema = {
'geometry': 'Polygon',
'properties': {'LABEL': 'float:16'}
}
with fiona.open('catchment.shp', 'w',
driver='ESRI Shapefile',
crs=grid.crs.srs,
schema=schema) as c:
i = 0
for shape, value in shapes:
rec = {}
rec['geometry'] = shape
rec['properties'] = {'LABEL' : str(value)}
rec['id'] = str(i)
c.write(rec)
i += 1
In [7]:
shp = gpd.read_file('catchment.shp')
In [8]:
shp
Out[8]:
In [9]:
fig, ax = plt.subplots(figsize=(6,6))
shp.plot(ax=ax)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Catchment polygon')
Out[9]:
In [10]:
# grid.extract_river_network outputs geojson, which is a little different
grid.accumulation(data='catch', dirmap=dirmap, pad_inplace=False, out_name='acc')
branches = grid.extract_river_network('catch', 'acc', threshold=50, dirmap=dirmap)
In [11]:
schema = {
'geometry': 'LineString',
'properties': {}
}
with fiona.open('rivers.shp', 'w',
driver='ESRI Shapefile',
crs=grid.crs.srs,
schema=schema) as c:
i = 0
for branch in branches['features']:
rec = {}
rec['geometry'] = branch['geometry']
rec['properties'] = {}
rec['id'] = str(i)
c.write(rec)
i += 1
In [12]:
shp = gpd.read_file('rivers.shp')
In [13]:
shp.head()
Out[13]:
In [14]:
fig, ax = plt.subplots(figsize=(6,6))
shp.plot(ax=ax)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('River network')
Out[14]:
In [ ]: