docs: def rasterToPolygon(raster_file, polygon_file):
          """
          Converts watershed raster to polygon and then dissolves it.
          It dissolves features based on the LINKNO attribute.
          """
          log("Process: Raster to Polygon ...")
          time_start = datetime.utcnow()
          temp_polygon_file = \
              "{0}_temp.shp".format(
                  os.path.splitext(os.path.basename(polygon_file))[0])
  
          GDALGrid(raster_file).to_polygon(out_shapefile=temp_polygon_file,
                                           fieldname="LINKNO",
                                           self_mask=True)
  
          log("Time to convert to polygon: {0}"
              .format(datetime.utcnow()-time_start))
  
          log("Dissolving ...")
          time_start_dissolve = datetime.utcnow()
          ogr_polygin_shapefile = ogr.Open(temp_polygon_file)
          ogr_polygon_shapefile_lyr = ogr_polygin_shapefile.GetLayer()
          number_of_features = ogr_polygon_shapefile_lyr.GetFeatureCount()
          polygon_rivid_list = np.zeros(number_of_features, dtype=np.int32)
          for feature_idx, catchment_feature in \
                  enumerate(ogr_polygon_shapefile_lyr):
              polygon_rivid_list[feature_idx] = \
                  catchment_feature.GetField('LINKNO')
  
          shp_drv = ogr.GetDriverByName('ESRI Shapefile')
          # Remove output shapefile if it already exists
          if os.path.exists(polygon_file):
              shp_drv.DeleteDataSource(polygon_file)
  
          dissolve_shapefile = shp_drv.CreateDataSource(polygon_file)
          dissolve_layer = \
              dissolve_shapefile.CreateLayer(
                  '',
                  ogr_polygon_shapefile_lyr.GetSpatialRef(),
                  ogr.wkbPolygon)
          dissolve_layer.CreateField(ogr.FieldDefn('LINKNO', ogr.OFTInteger))
          dissolve_layer_defn = dissolve_layer.GetLayerDefn()
  
          for unique_rivid in np.unique(polygon_rivid_list):
              # get indices where it is in the polygon
              feature_indices = np.where(polygon_rivid_list == unique_rivid)[0]
              new_feat = ogr.Feature(dissolve_layer_defn)
              new_feat.SetField('LINKNO', int(unique_rivid))
  
              if len(feature_indices) == 1:
                  # write feature to file
                  feature = \
                      ogr_polygon_shapefile_lyr.GetFeature(feature_indices[0])
                  new_feat.SetGeometry(feature.GetGeometryRef())
              else:
                  # dissolve
                  dissolve_poly_list = []
                  for feature_index in feature_indices:
                      feature = \
                          ogr_polygon_shapefile_lyr.GetFeature(feature_index)
                      feat_geom = feature.GetGeometryRef()
                      dissolve_poly_list.append(
                          shapely_loads(feat_geom.ExportToWkb()))
                  dissolve_polygon = cascaded_union(dissolve_poly_list)
                  new_feat.SetGeometry(
                      ogr.CreateGeometryFromWkb(dissolve_polygon.wkb))
              dissolve_layer.CreateFeature(new_feat)
          # clean up
          shp_drv.DeleteDataSource(temp_polygon_file)
          log("Time to dissolve: {0}".format(datetime.utcnow() -
                                             time_start_dissolve))
          log("Total time to convert: {0}".format(datetime.utcnow() -
                                                  time_start))