docs: def CreateNetworkConnectivityNHDPlus(in_drainage_line,
                                       out_connectivity_file,
                                       file_geodatabase=None):
      """
      Creates Network Connectivity input CSV file for RAPID
      based on the NHDPlus drainage lines with
      COMID, FROMNODE, TONODE, and DIVERGENCE fields.
  
      Parameters
      ----------
      in_drainage_line: str
          Path to the stream network (i.e. Drainage Line) shapefile.
      out_connectivity_file: str
          The path to the output connectivity file.
      file_geodatabase: str, optional
          Path to the file geodatabase. If you use this option,
          in_drainage_line is the name of the stream network feature class
          (WARNING: Not always stable with GDAL).
  
  
      Example::
  
          from RAPIDpy.gis.network import CreateNetworkConnectivityNHDPlus
  
          CreateNetworkConnectivityNHDPlus(
              in_drainage_line='/path/to/drainageline.shp',
              out_connectivity_file='/path/to/rapid_connect.csv')
  
      """
      ogr_drainage_line_shapefile_lyr, ogr_drainage_line_shapefile = \
          open_shapefile(in_drainage_line, file_geodatabase)
  
      ogr_drainage_line_definition = \
          ogr_drainage_line_shapefile_lyr.GetLayerDefn()
  
      orig_field_names = []
      for idx in xrange(ogr_drainage_line_definition.GetFieldCount()):
          orig_field_names.append(
              ogr_drainage_line_definition.GetFieldDefn(idx).GetName())
  
      upper_field_names = [field.upper() for field in orig_field_names]
  
      def get_field_name_index(upper_field_name, _upper_field_names):
          """
          returns index of field name
          """
          try:
              return _upper_field_names.index(upper_field_name)
          except ValueError:
              raise IndexError("{0} not found in shapefile .."
                               .format(_upper_field_names))
  
      rivid_field = \
          orig_field_names[get_field_name_index('COMID', upper_field_names)]
      fromnode_field = \
          orig_field_names[get_field_name_index('FROMNODE', upper_field_names)]
      tonode_field = \
          orig_field_names[get_field_name_index('TONODE', upper_field_names)]
      divergence_field =\
          orig_field_names[get_field_name_index('DIVERGENCE', upper_field_names)]
  
      number_of_features = ogr_drainage_line_shapefile_lyr.GetFeatureCount()
      rivid_list = np.zeros(number_of_features, dtype=np.int32)
      fromnode_list = np.zeros(number_of_features, dtype=np.int32)
      tonode_list = np.zeros(number_of_features, dtype=np.int32)
      divergence_list = np.zeros(number_of_features, dtype=np.int32)
      for feature_idx, catchment_feature in \
              enumerate(ogr_drainage_line_shapefile_lyr):
          rivid_list[feature_idx] = catchment_feature.GetField(rivid_field)
          fromnode_list[feature_idx] = catchment_feature.GetField(fromnode_field)
          tonode_list[feature_idx] = catchment_feature.GetField(tonode_field)
          divergence_list[feature_idx] = \
              catchment_feature.GetField(divergence_field)
  
      del ogr_drainage_line_shapefile
      # -------------------------------------------------------------------------
      # Compute connectivity, based on:
      # https://github.com/c-h-david/rrr/blob/master/src/rrr_riv_tot_gen_all_nhdplus.py
      # -------------------------------------------------------------------------
      fromnode_list[fromnode_list == 0] = -9999
      # Some NHDPlus v1 reaches have FLOWDIR='With Digitized'
      # but no info in VAA table
  
      fromnode_list[divergence_list == 2] = -9999
      # Virtually disconnect the upstream node of all minor divergences
      del divergence_list  # delete information in list
  
      next_down_id_list = np.zeros(number_of_features, dtype=np.int32)
      for rivid_index in xrange(len(rivid_list)):
          try:
              next_down_id_list[rivid_index] = \
                  rivid_list[
                      np.where(fromnode_list == tonode_list[rivid_index])[0][0]]
          except IndexError:
              # this is an outlet
              next_down_id_list[rivid_index] = -1
  
      # determine the downstream reach for each reach
  
      # empty unecessary lists
      del fromnode_list
      del tonode_list
  
      StreamIDNextDownIDToConnectivity(rivid_list,
                                       next_down_id_list,
                                       out_connectivity_file)