Merge pull request #162 from mdbartos/numba

Add numba acceleration and major refactor

docs/accumulation.md

@@ -5,14 +5,15 @@
 The `grid.accumulation` method operates on a flow direction grid. This flow direction grid can be computed from a DEM, as shown in [flow directions](https://mdbartos.github.io/pysheds/flow-directions.html).
 
 ```python
->>> from pysheds.grid import Grid
+from pysheds.grid import Grid
 
 # Instantiate grid from raster
->>> grid = Grid.from_raster('../data/dem.tif', data_name='dem')
+grid = Grid.from_raster('./data/dem.tif')
+dem = grid.read_raster('./data/dem.tif')
 
 # Resolve flats and compute flow directions
->>> grid.resolve_flats(data='dem', out_name='inflated_dem')
->>> grid.flowdir('inflated_dem', out_name='dir')
+inflated_dem = grid.resolve_flats(dem)
+fdir = grid.flowdir(inflated_dem)
 ```
 
 ## Computing accumulation
@@ -21,29 +22,34 @@ Accumulation is computed using the `grid.accumulation` method.
 
 ```python
 # Compute accumulation
->>> grid.accumulation(data='dir', out_name='acc')
-
-# Plot accumulation
->>> acc = grid.view('acc')
->>> plt.imshow(acc)
+acc = grid.accumulation(fdir)
 ```
 
-![Full accumulation](https://s3.us-east-2.amazonaws.com/pysheds/img/full_accumulation.png)
-
-## Computing weighted accumulation
-
-Weights can be used to adjust the relative contribution of each cell.
+<details>
+<summary>Plotting code...</summary>
+<p>
 
 ```python
-import pyproj
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+%matplotlib inline
+
+fig, ax = plt.subplots(figsize=(8,6))
+fig.patch.set_alpha(0)
+plt.grid('on', zorder=0)
+im = ax.imshow(acc, extent=grid.extent, zorder=2,
+               cmap='cubehelix',
+               norm=colors.LogNorm(1, acc.max()),
+               interpolation='bilinear')
+plt.colorbar(im, ax=ax, label='Upstream Cells')
+plt.title('Flow Accumulation', size=14)
+plt.xlabel('Longitude')
+plt.ylabel('Latitude')
+plt.tight_layout()
+```
 
-# Compute areas of each cell in new projection
-new_crs = pyproj.Proj('+init=epsg:3083')
-areas = grid.cell_area(as_crs=new_crs, inplace=False)
+</p>
+</details>
 
-# Weight each cell by its relative area
-weights = (areas / areas.max()).ravel()
 
-# Compute accumulation with new weights
-grid.accumulation(data='dir', weights=weights, out_name='acc')
-```
+![Full accumulation](https://s3.us-east-2.amazonaws.com/pysheds/img/acc_acc.png)

docs/catchment.md

@@ -5,53 +5,139 @@
 The `grid.catchment` method operates on a flow direction grid. This flow direction grid can be computed from a DEM, as shown in [flow directions](https://mdbartos.github.io/pysheds/flow-directions.html).
 
 ```python
->>> from pysheds.grid import Grid
+from pysheds.grid import Grid
 
 # Instantiate grid from raster
->>> grid = Grid.from_raster('../data/dem.tif', data_name='dem')
+grid = Grid.from_raster('./data/dem.tif')
+dem = grid.read_raster('./data/dem.tif')
 
 # Resolve flats and compute flow directions
->>> grid.resolve_flats(data='dem', out_name='inflated_dem')
->>> grid.flowdir('inflated_dem', out_name='dir')
+inflated_dem = grid.resolve_flats(dem)
+fdir = grid.flowdir(inflated_dem)
 ```
 
 ## Delineating the catchment
 
-To delineate a catchment, first specify a pour point (the outlet of the catchment). If the x and y components of the pour point are spatial coordinates in the grid's spatial reference system, specify `xytype='label'`.
+To delineate a catchment, first specify a pour point (the outlet of the catchment). If the x and y components of the pour point are spatial coordinates in the grid's spatial reference system, specify `xytype='coordinate'`.
 
 ```python
 # Specify pour point
->>> x, y = -97.294167, 32.73750
+x, y = -97.294167, 32.73750
 
 # Delineate the catchment
->>> grid.catchment(data='dir', x=x, y=y, out_name='catch',
-                   recursionlimit=15000, xytype='label')
+catch = grid.catchment(x=x, y=y, fdir=fdir, xytype='coordinate')
 
 # Plot the result
->>> grid.clip_to('catch')
->>> plt.imshow(grid.view('catch'))
+grid.clip_to(catch)
+catch_view = grid.view(catch)
 ```
 
-![Delineated catchment](https://s3.us-east-2.amazonaws.com/pysheds/img/catchment.png)
+<details>
+<summary>Plotting code...</summary>
+<p>
+
+```python
+# Plot the catchment
+fig, ax = plt.subplots(figsize=(8,6))
+fig.patch.set_alpha(0)
+
+plt.grid('on', zorder=0)
+im = ax.imshow(np.where(catch_view, catch_view, np.nan), extent=grid.extent,
+               zorder=1, cmap='Greys_r')
+plt.xlabel('Longitude')
+plt.ylabel('Latitude')
+plt.title('Delineated Catchment', size=14)
+```
+
+</p>
+</details>
+
+
+![Delineated catchment](https://s3.us-east-2.amazonaws.com/pysheds/img/catch.png)
 
 If the x and y components of the pour point correspond to the row and column indices of the flow direction array, specify `xytype='index'`:
 
 ```python
 # Reset the view
->>> grid.clip_to('dir')
+grid.viewfinder = fdir.viewfinder
 
 # Find the row and column index corresponding to the pour point
->>> col, row = grid.nearest_cell(x, y)
->>> col, row
-(229, 101)
+col, row = grid.nearest_cell(x, y)
+
+# Delineate the catchment
+catch = grid.catchment(x=col, y=row, fdir=fdir, xytype='index')
+
+# Plot the result
+grid.clip_to(catch)
+catch_view = grid.view(catch)
+```
+
+<details>
+<summary>Plotting code...</summary>
+<p>
+
+```python
+# Plot the catchment
+fig, ax = plt.subplots(figsize=(8,6))
+fig.patch.set_alpha(0)
+
+plt.grid('on', zorder=0)
+im = ax.imshow(np.where(catch_view, catch_view, np.nan), extent=grid.extent,
+               zorder=1, cmap='Greys_r')
+plt.xlabel('Longitude')
+plt.ylabel('Latitude')
+plt.title('Delineated Catchment', size=14)
+```
+
+</p>
+</details>
+
+
+![Delineated catchment index](https://s3.us-east-2.amazonaws.com/pysheds/img/catch.png)
+
+## Snapping pour point to high accumulation cells
+
+Sometimes the pour point isn't known exactly. In this case, it can be helpful to first compute the accumulation and then snap a trial pour point to the nearest high accumulation cell.
+
+```python
+# Reset view
+grid.viewfinder = fdir.viewfinder
+
+# Compute accumulation
+acc = grid.accumulation(fdir)
+
+# Snap pour point to high accumulation cell
+x_snap, y_snap = grid.snap_to_mask(acc > 1000, (x, y))
+
 
 # Delineate the catchment
->>> grid.catchment(data=grid.dir, x=col, y=row, out_name='catch',
-                   recursionlimit=15000, xytype='index')
+catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir, xytype='coordinate')
 
 # Plot the result
->>> grid.clip_to('catch')
->>> plt.imshow(grid.view('catch'))
+grid.clip_to(catch)
+catch_view = grid.view(catch)
+```
+
+<details>
+<summary>Plotting code...</summary>
+<p>
+
+```python
+# Plot the catchment
+fig, ax = plt.subplots(figsize=(8,6))
+fig.patch.set_alpha(0)
+
+plt.grid('on', zorder=0)
+im = ax.imshow(np.where(catch_view, catch_view, np.nan), extent=grid.extent,
+               zorder=1, cmap='Greys_r')
+plt.xlabel('Longitude')
+plt.ylabel('Latitude')
+plt.title('Delineated Catchment', size=14)
 ```
 
-![Delineated catchment index](https://s3.us-east-2.amazonaws.com/pysheds/img/catchment.png)
+</p>
+</details>
+
+![Delineated catchment snap](https://s3.us-east-2.amazonaws.com/pysheds/img/catch.png)
+
+