# Read elevation and flow direction rasters
    # ----------------------------
    from pysheds.grid import Grid

    grid = Grid.from_raster('n30w100_con', data_name='dem')
    grid.read_raster('n30w100_dir', data_name='dir')
    grid.view('dem')

    # Delineate a catchment
    # ---------------------
    # Specify pour point
    x, y = -97.294167, 32.73750
    # Specify directional mapping
    dirmap=(64, 128, 1, 2, 4, 8, 16, 32)

    # Delineate the catchment
    grid.catchment(data='dir', x=x, y=y, dirmap=dirmap, out_name='catch',
               recursionlimit=15000, xytype='label')

    # Crop and plot the catchment
    # ---------------------------
    # Clip the bounding box to the catchment
    grid.clip_to('catch')
    grid.view('catch')

    # Calculate flow accumulation
    # --------------------------
    grid.accumulation(data='catch', dirmap=dirmap, out_name='acc')
    grid.view('acc')

    # Calculate distance to outlet from each cell
    # -------------------------------------------
    grid.flow_distance(data='catch', x=x, y=y, dirmap=dirmap,
                   out_name='dist', xytype='label')
    grid.view('dist')

    # Extract river network
    # ---------------------
    branches = grid.extract_river_network(fdir='catch', acc='acc',
                                          threshold=50, dirmap=dirmap)

    # Combine with land cover data
    # ---------------------
    grid.read_raster('nlcd_2011_impervious_2011_edition_2014_10_10.img',
                      data_name='terrain', window=grid.bbox, window_crs=grid.crs)
    grid.view('terrain')

    # Convert catchment raster to vector and combine with soils shapefile
    # ---------------------
    # Read soils shapefile
    import geopandas as gpd
    from shapely import geometry, ops
    soils = gpd.read_file('nrcs-soils-tarrant_439.shp')
    # Convert catchment raster to vector geometry and find intersection
    shapes = grid.polygonize()
    catchment_polygon = ops.unary_union([geometry.shape(shape)
                                         for shape, value in shapes])
    soils = soils[soils.intersects(catchment_polygon)]
    catchment_soils = soils.intersection(catchment_polygon)

    # Convert soils polygons to raster
    # ---------------------
    soil_polygons = zip(catchment_soils.geometry.values,
                        catchment_soils['soil_type'].values)
    soil_raster = grid.rasterize(soil_polygons, fill=np.nan)