docs: def addLengthMeters(stream_network):
          """
          Adds length field in meters to network
          (The added field name will be 'LENGTH_M').
  
          .. note:: This may be needed for generating the kfac file
                    depending on the units of your raster. See: :doc:`gis_tools`.
  
          Parameters
          ----------
          stream_network: str
              Path to stream network file.
  
  
          Here is an example of how to use this:
  
          .. code:: python
  
              import os
              from RAPIDpy.gis.taudem import TauDEM
  
              output_directory = '/path/to/output/files'
              TauDEM.addLengthMeters(os.path.join(output_directory,
                                                  "stream_reach_file.shp"))
  
          """
          network_shapefile = ogr.Open(stream_network, 1)
          network_layer = network_shapefile.GetLayer()
          network_layer_defn = network_layer.GetLayerDefn()
  
          # make sure projection EPSG:4326
          network_layer_proj = network_layer.GetSpatialRef()
          geographic_proj = osr.SpatialReference()
          geographic_proj.ImportFromEPSG(4326)
          proj_transform = None
          if network_layer_proj != geographic_proj:
              proj_transform = osr.CoordinateTransformation(network_layer_proj,
                                                            geographic_proj)
  
          # check for field
          create_field = True
          for i in xrange(network_layer_defn.GetFieldCount()):
              field_name = network_layer_defn.GetFieldDefn(i).GetName()
              if field_name == 'LENGTH_M':
                  create_field = False
                  break
  
          if create_field:
              network_layer.CreateField(ogr.FieldDefn('LENGTH_M', ogr.OFTReal))
  
          geo_manager = Geod(ellps="WGS84")
          for network_feature in network_layer:
              feat_geom = network_feature.GetGeometryRef()
              # make sure coordinates are geographic
              if proj_transform:
                  feat_geom.Transform(proj_transform)
  
              line = shapely_loads(feat_geom.ExportToWkb())
              lon_list, lat_list = line.xy
              dist = geo_manager.inv(lon_list[:-1], lat_list[:-1],
                                     lon_list[1:], lat_list[1:])[2]
              network_feature.SetField('LENGTH_M', sum(dist))
              network_layer.SetFeature(network_feature)