pysheds_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.morphology import skeletonize\n",
    "from scipy import ndimage\n",
    "import pyproj\n",
    "import numpy as np\n",
    "import fiona\n",
    "from rasterio import features\n",
    "from shapely import geometry\n",
    "from pysheds.grid import Grid\n",
    "\n",
    "\n",
    "id_column = \"station_id\"\n",
    "id_value = \"08LF061\"\n",
    "pourpoint_coord = [1310870.36276, 647616.36174]\n",
    "grid = Grid.from_raster(\"dem_08LF061.tif\", data_name='dem')\n",
    "\n",
    "# load upstream FWA stream segment\n",
    "with fiona.open(\"wsdrefine_streams.shp\") as src:\n",
    "    stream_features = list(filter(lambda f: f['properties'][id_column] == id_value, src))\n",
    "\n",
    "# convert stream geojson features to shapely shapes\n",
    "stream_shapes = [geometry.shape(f['geometry']) for f in stream_features]\n",
    "\n",
    "# convert shapes to raster\n",
    "stream_raster = features.rasterize(\n",
    "    ((g, 1) for g in stream_shapes),\n",
    "    out_shape=grid.shape,\n",
    "    transform=grid.affine,\n",
    "    all_touched=False\n",
    ")\n",
    "stream_raster = skeletonize(stream_raster).astype(np.uint8)\n",
    "\n",
    "# Create boolean mask based on rasterized river shapes\n",
    "mask = stream_raster.astype(np.bool)\n",
    "\n",
    "# Blur mask using a gaussian filter\n",
    "blurred_mask = ndimage.filters.gaussian_filter(mask.astype(np.float64), sigma=2.5)\n",
    "\n",
    "# Set central channel to max to prevent pits\n",
    "blurred_mask[mask.astype(np.bool)] = blurred_mask.max()\n",
    "\n",
    "# Set elevation change for burned cells\n",
    "dz = 16.5\n",
    "\n",
    "# Set mask to blurred mask\n",
    "mask = blurred_mask\n",
    "\n",
    "# Create a view onto the DEM array\n",
    "dem = grid.view('dem', dtype=np.float64, nodata=np.nan)\n",
    "\n",
    "# Subtract dz where mask is nonzero\n",
    "dem[(mask > 0)] -= dz*mask[(mask > 0)]\n",
    "\n",
    "# defining the crs used improves results\n",
    "new_crs = pyproj.Proj('+init=epsg:3005')\n",
    "\n",
    "#         N    NE    E    SE    S    SW    W    NW\n",
    "dirmap = (64,  128,  1,   2,    4,   8,    16,  32)\n",
    "\n",
    "# fill / resolve flats / flow direction / accumulation\n",
    "grid.fill_depressions(data=dem, out_name='flooded_dem')\n",
    "grid.resolve_flats(data='flooded_dem', out_name='inflated_dem')\n",
    "grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap, as_crs=new_crs)\n",
    "grid.accumulation(data='dir', dirmap=dirmap, out_name='acc', apply_mask=False)\n",
    "\n",
    "# snap pour point to higher accumulation cells\n",
    "# (in theory, this shouldn't really be necesssary after burning in the\n",
    "# streams, but the DEM and streams are not exact matches)\n",
    "x, y = pourpoint_coord\n",
    "xy_snapped = grid.snap_to_mask(grid.acc > 50, [[x, y]], return_dist=False)\n",
    "x, y = xy_snapped[0][0], xy_snapped[0][1]\n",
    "\n",
    "# create catchment\n",
    "grid.catchment(data='dir', x=x, y=y, dirmap=dirmap, out_name='catch',\n",
    "               recursionlimit=15000, xytype='label', nodata_out=0)\n",
    "\n",
    "# Clip the bounding box to the catchment\n",
    "grid.clip_to('catch')\n",
    "\n",
    "polys = [geometry.shape(shape) for shape, value in grid.polygonize()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"300\" height=\"300\" viewBox=\"1310844.8238010006 647548.1671408139 923.0131980003789 492.88842207426205\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,1295589.222703702)\"><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"6.153421320002526\" opacity=\"0.6\" d=\"M 1311633.1052250008,648006.8698888881 L 1311633.1052250008,647981.8982962955 L 1311607.9687000008,647981.8982962955 L 1311457.1495500007,647981.8982962955 L 1311457.1495500007,647956.9267037029 L 1311432.0130250007,647956.9267037029 L 1311381.7399750007,647956.9267037029 L 1311381.7399750007,647931.9551111104 L 1311356.6034500008,647931.9551111104 L 1311356.6034500008,647906.9835185177 L 1311331.4669250008,647906.9835185177 L 1311256.0573500008,647906.9835185177 L 1311256.0573500008,647882.0119259252 L 1311230.9208250009,647882.0119259252 L 1311205.7843000006,647882.0119259252 L 1311205.7843000006,647857.0403333325 L 1311180.6477750007,647857.0403333325 L 1311180.6477750007,647832.06874074 L 1311155.5112500007,647832.06874074 L 1311155.5112500007,647807.0971481473 L 1311130.3747250007,647807.0971481473 L 1311130.3747250007,647782.1255555548 L 1311105.2382000007,647782.1255555548 L 1311080.1016750007,647782.1255555548 L 1311080.1016750007,647757.1539629622 L 1311054.9651500008,647757.1539629622 L 1311054.9651500008,647732.1823703696 L 1311029.8286250008,647732.1823703696 L 1311029.8286250008,647707.210777777 L 1311004.6921000008,647707.210777777 L 1311004.6921000008,647682.2391851844 L 1310979.5555750006,647682.2391851844 L 1310979.5555750006,647657.2675925918 L 1310954.4190500006,647657.2675925918 L 1310954.4190500006,647632.2959999992 L 1310929.2825250006,647632.2959999992 L 1310929.2825250006,647607.3244074066 L 1310904.1460000006,647607.3244074066 L 1310879.0094750007,647607.3244074066 L 1310879.0094750007,647582.352814814 L 1310954.4190500006,647582.352814814 L 1310954.4190500006,647607.3244074066 L 1311029.8286250008,647607.3244074066 L 1311029.8286250008,647582.352814814 L 1311130.3747250007,647582.352814814 L 1311130.3747250007,647607.3244074066 L 1311155.5112500007,647607.3244074066 L 1311155.5112500007,647632.2959999992 L 1311180.6477750007,647632.2959999992 L 1311180.6477750007,647657.2675925918 L 1311205.7843000006,647657.2675925918 L 1311205.7843000006,647682.2391851844 L 1311230.9208250009,647682.2391851844 L 1311230.9208250009,647707.210777777 L 1311256.0573500008,647707.210777777 L 1311256.0573500008,647732.1823703696 L 1311281.1938750008,647732.1823703696 L 1311281.1938750008,647757.1539629622 L 1311306.3304000008,647757.1539629622 L 1311306.3304000008,647782.1255555548 L 1311331.4669250008,647782.1255555548 L 1311331.4669250008,647807.0971481473 L 1311356.6034500008,647807.0971481473 L 1311356.6034500008,647832.06874074 L 1311532.5591250008,647832.06874074 L 1311557.6956500008,647832.06874074 L 1311557.6956500008,647857.0403333325 L 1311633.1052250008,647857.0403333325 L 1311658.2417500007,647857.0403333325 L 1311658.2417500007,647931.9551111104 L 1311683.3782750007,647931.9551111104 L 1311683.3782750007,647956.9267037029 L 1311708.514800001,647956.9267037029 L 1311708.514800001,647981.8982962955 L 1311733.651325001,647981.8982962955 L 1311733.651325001,648006.8698888881 L 1311633.1052250008,648006.8698888881 z\" /></g></svg>"
      ],
      "text/plain": [
       "<shapely.geometry.polygon.Polygon at 0x109935ba8>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polys[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}