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

geotrellis.feature.op.geometry.IDWInterpolate.scala Maven / Gradle / Ivy

The newest version!
package geotrellis.feature.op.geometry

import geotrellis._
import geotrellis.feature._
import geotrellis.raster._

import spire.syntax.cfor._

/**
 * IDW Interpolation
 */
case class IDWInterpolate(points:Op[Seq[Point[Int]]],re:Op[RasterExtent],radius:Op[Option[Int]]=None)
    extends Op3(points,re,radius)({
  (points,re,radius) =>
    val cols = re.cols
    val rows = re.rows
    val data = RasterData.emptyByType(TypeInt, cols, rows)
    if(points.isEmpty) {
      Result(Raster(data,re))
    } else {



      val r = radius match {
        case Some(r) =>
          val rr = r*r
          val index:SpatialIndex[Point[Int]] = SpatialIndex(points)(p => (p.x,p.y))
          cfor(0)(_ < cols, _ + 1) { col =>
            cfor(0)(_ < rows, _ + 1) { row =>
              val destX = re.gridColToMap(col)
              val destY = re.gridRowToMap(row)
              val pts = index.pointsInExtent(Extent(destX - r, destY - r, destX + r, destY + r))

              if (pts.isEmpty) {
                data.set(col, row, NODATA)
              } else {
                var s = 0.0
                var c = 0
                var ws = 0.0
                val length = pts.length

                cfor(0)(_ < length, _ + 1) { i =>
                  val point = pts(i)
                  val dX = (destX - point.x)
                  val dY = (destY - point.y)
                  val d = dX * dX + dY * dY
                  if (d < rr) {
                    val w = 1 / d
                    s += point.data * w
                    ws += w
                    c += 1
                  }
                }

                if (c == 0) {
                  data.set(col, row, NODATA)
                } else {
                  val mean = s / ws
                  data.set(col, row, mean.toInt)
                }
              }
            }
          }
        case None =>
          val length = points.length
          cfor(0)(_ < cols, _ + 1) { col =>
            cfor(0)(_ < rows, _ + 1) { row =>
              val destX = re.gridColToMap(col)
              val destY = re.gridRowToMap(row)
              var s = 0.0
              var c = 0
              var ws = 0.0

              cfor(0)(_ < length, _ + 1) { i =>
                val point = points(i)
                val dX = (destX - point.x)
                val dY = (destY - point.y)
                val d = dX * dX + dY * dY
                val w = 1 / d
                s += point.data * w
                ws += w
                c += 1
              }

              if (c == 0) {
                data.set(col, row, NODATA)
              } else {
                val mean = s / ws
                data.set(col, row, mean.toInt)
              }
            }
          }
      }
      Result(Raster(data,re))
    }
})




© 2015 - 2024 Weber Informatics LLC | Privacy Policy