docs: def serialize_streamnet(streamnet_file, output_reach_file):
          """Eliminate reach with zero length and return the reach ID map.
          Args:
              streamnet_file: original stream net ESRI shapefile
              output_reach_file: serialized stream net, ESRI shapefile
  
          Returns:
              id pairs {origin: newly assigned}
          """
          FileClass.copy_files(streamnet_file, output_reach_file)
          ds_reach = ogr_Open(output_reach_file, update=True)
          layer_reach = ds_reach.GetLayer(0)
          layer_def = layer_reach.GetLayerDefn()
          i_link = layer_def.GetFieldIndex(FLD_LINKNO)
          i_link_downslope = layer_def.GetFieldIndex(FLD_DSLINKNO)
          i_len = layer_def.GetFieldIndex(REACH_LENGTH)
  
          old_id_list = []
          # there are some reaches with zero length.
          # this program will remove these zero-length reaches
          # output_dic is used to store the downstream reaches of these zero-length
          # reaches
          output_dic = {}
          ft = layer_reach.GetNextFeature()
          while ft is not None:
              link_id = ft.GetFieldAsInteger(i_link)
              reach_len = ft.GetFieldAsDouble(i_len)
              if link_id not in old_id_list:
                  if reach_len < DELTA:
                      downstream_id = ft.GetFieldAsInteger(i_link_downslope)
                      output_dic[link_id] = downstream_id
                  else:
                      old_id_list.append(link_id)
  
              ft = layer_reach.GetNextFeature()
          old_id_list.sort()
  
          id_map = {}
          for i, old_id in enumerate(old_id_list):
              id_map[old_id] = i + 1
          # print(id_map)
          # change old ID to new ID
          layer_reach.ResetReading()
          ft = layer_reach.GetNextFeature()
          while ft is not None:
              link_id = ft.GetFieldAsInteger(i_link)
              if link_id not in id_map:
                  layer_reach.DeleteFeature(ft.GetFID())
                  ft = layer_reach.GetNextFeature()
                  continue
  
              ds_id = ft.GetFieldAsInteger(i_link_downslope)
              ds_id = output_dic.get(ds_id, ds_id)
              ds_id = output_dic.get(ds_id, ds_id)
  
              ft.SetField(FLD_LINKNO, id_map[link_id])
              if ds_id in id_map:
                  ft.SetField(FLD_DSLINKNO, id_map[ds_id])
              else:
                  # print(ds_id)
                  ft.SetField(FLD_DSLINKNO, -1)
              layer_reach.SetFeature(ft)
              ft = layer_reach.GetNextFeature()
          ds_reach.ExecuteSQL(str('REPACK reach'))
          layer_reach.SyncToDisk()
          ds_reach.Destroy()
          del ds_reach
          return id_map