code: def extract_line(geom, dem, **kwargs):
      
      kwargs.setdefault('masked', True)
  
      coords_in = coords_array(geom)
  
      # Transform geometry into pixels
      f = lambda *x: ~dem.transform * x
      px = transform(f,geom)
  
      # Subdivide geometry segments if option is given
      interval = kwargs.pop('subdivide', 1)
      if interval is not None:
          px = subdivide(px, interval=interval)
  
      # Transform pixels back to geometry
      # to capture subdivisions
      f = lambda *x: dem.transform * (x[0],x[1])
      geom = transform(f,px)
  
      # Get min and max coords for windowing
      # Does not deal with edge cases where points
      # are outside of footprint of DEM
  
      coords_px = coords_array(px)
      mins = N.floor(coords_px.min(axis=0))
      maxs = N.ceil(coords_px.max(axis=0))
  
      window = tuple((int(mn),int(mx))
          for mn,mx in zip(mins[::-1],maxs[::-1]))
  
      aff = Affine.translation(*(-mins))
  
      f = lambda *x: aff * x
      px_to_extract = transform(f,px)
  
      band = dem.read(1, window=window, **kwargs)
      extracted = bilinear(band, px_to_extract)
      coords = coords_array(extracted)
  
      coords[:,:2] = coords_array(geom)
      return coords
docs: Extract a linear feature from a `rasterio` geospatial dataset.