All Downloads are FREE. Search and download functionalities are using the official Maven repository.

geotrellis.raster.interpolation.InverseDistanceWeighting.scala Maven / Gradle / Ivy

package geotrellis.raster.interpolation

import geotrellis.raster._
import geotrellis.vector._
import spire.syntax.cfor._

object InverseDistanceWeighting {
  /**
    * Compute an Inverse Distance Weighting raster over the given
    * extent from the given set known-points.  Please see
    * https://en.wikipedia.org/wiki/Inverse_distance_weighting for
    * more details.
    *
    * @param   points            A collection of known-points
    * @param   rasterExtent      The study area
    * @param   radius            Interpolation radius in coordinate unit
    * @param   equalWeightRadius If any points lie at that distance from a cell, the interpolation result for it would
    *                            be the mean of those points
    * @param   cellType          Interpolation radius in coordinate unit
    * @param   onSet             Function to be applied before storing interpolation values
    * @return The data interpolated across the study area
    */
  def apply[D](
                points: Seq[PointFeature[D]],
                rasterExtent: RasterExtent,
                radius: Double = Double.PositiveInfinity,
                equalWeightRadius: Double = 0,
                cellType: CellType = IntConstantNoDataCellType,
                onSet: Double => Double = x => x
              )(implicit ev: D => Double): Tile = {
    val cols = rasterExtent.cols
    val rows = rasterExtent.rows
    val tile = ArrayTile.empty(cellType, cols, rows)

    if (points.isEmpty) {
      tile
    } else {
      val ewr2 = equalWeightRadius * equalWeightRadius

      def idw(points: Seq[PointFeature[D]], x: Double, y: Double)(filter: Double => Boolean) = {
        var sum = 0.0
        var count = 0
        var weightSum = 0.0
        var sampleSum = 0.0
        var sampleCount = 0

        for (point <- points) {
          val dX = x - point.geom.x
          val dY = y - point.geom.y
          val d2 = dX * dX + dY * dY

          if (filter(d2)) {
            if (d2 <= ewr2) {
              sampleSum += point.data
              sampleCount += 1
            } else {
              val w = 1 / d2
              sum += point.data * w
              weightSum += w
              count += 1
            }
          }
        }

        if (sampleCount == 0) {
          if (count == 0) {
            doubleNODATA
          } else {
            onSet(sum / weightSum)
          }
        } else {
          onSet(sampleSum / sampleCount)
        }
      }

      radius match {
        case Double.PositiveInfinity =>
          cfor(0)(_ < rows, _ + 1) { row =>
            cfor(0)(_ < cols, _ + 1) { col =>
              val x = rasterExtent.gridColToMap(col)
              val y = rasterExtent.gridRowToMap(row)
              val value = idw(points, x, y)(d2 => true)

              tile.setDouble(col, row, value)
            }
          }
        case _ =>
          val r2 = radius * radius
          val index: SpatialIndex[PointFeature[D]] = SpatialIndex(points)(p => (p.geom.x, p.geom.y))

          cfor(0)(_ < rows, _ + 1) { row =>
            cfor(0)(_ < cols, _ + 1) { col =>
              val x = rasterExtent.gridColToMap(col)
              val y = rasterExtent.gridRowToMap(row)
              val rPoints = index.pointsInExtent(Extent(x - radius, y - radius, x + radius, y + radius))
              val value = if (rPoints.isEmpty) {
                doubleNODATA
              } else {
                idw(rPoints, x, y)(d2 => d2 <= r2)
              }

              tile.setDouble(col, row, value)
            }
          }
      }
      tile
    }
  }
}




© 2015 - 2025 Weber Informatics LLC | Privacy Policy