API refactoring for Rasterizing Reproject

Create RasterRegionReproject type class
locationtech · Dec 12, 2017 · 0451152 · 0451152
1 parent 9ef7c64
commit 0451152
Show file tree

Hide file tree

Showing 15 changed files with 134 additions and 263 deletions.
diff --git a/raster/src/main/scala/geotrellis/raster/package.scala b/raster/src/main/scala/geotrellis/raster/package.scala
@@ -98,13 +98,11 @@ package object raster
 
   implicit class withSinglebandRasterMethods(val self: SinglebandRaster) extends MethodExtensions[SinglebandRaster]
       with reproject.SinglebandRasterReprojectMethods
-      with reproject.SinglebandRasterRasterizeReprojectMethods
       with resample.SinglebandRasterResampleMethods
       with vectorize.SinglebandRasterVectorizeMethods
 
   implicit class withMultibandRasterMethods(val self: MultibandRaster) extends MethodExtensions[MultibandRaster]
       with reproject.MultibandRasterReprojectMethods
-      with reproject.MultibandRasterRasterizeReprojectMethods
       with resample.MultibandRasterResampleMethods
 
   implicit class withTileSeqMethods(val self: Traversable[Tile]) extends MethodExtensions[Traversable[Tile]]

diff --git a/raster/src/main/scala/geotrellis/raster/reproject/Implicits.scala b/raster/src/main/scala/geotrellis/raster/reproject/Implicits.scala
@@ -22,11 +22,4 @@ object Implicits extends Implicits
 
 trait Implicits {
   implicit class withProjectedRasterReprojectMethods[T <: CellGrid](self: ProjectedRaster[T])(implicit ev: Raster[T] => RasterReprojectMethods[Raster[T]])
-    extends ProjectedRasterReprojectMethods[T](self)
-
-  implicit class withProjectedRasterRasterizeReprojectMethods[T <: CellGrid](
-    self: ProjectedRaster[T]
-  )(
-    implicit ev: Raster[T] => RasterRasterizeReprojectMethods[T]
-  ) extends ProjectedRasterRasterizeReprojectMethods[T](self)
 }
diff --git a/raster/src/main/scala/geotrellis/raster/reproject/MultibandRasterReprojectMethods.scala b/raster/src/main/scala/geotrellis/raster/reproject/MultibandRasterReprojectMethods.scala
@@ -18,7 +18,8 @@ package geotrellis.raster.reproject
 
 import geotrellis.raster._
 import geotrellis.raster.resample._
-import geotrellis.vector.Extent
+import geotrellis.raster.rasterize._
+import geotrellis.vector.Polygon
 import geotrellis.proj4._
 
 import spire.syntax.cfor._
@@ -34,11 +35,9 @@ trait MultibandRasterReprojectMethods extends RasterReprojectMethods[MultibandRa
   ): MultibandRaster = {
     val Raster(tile, extent) = self
     val bands =
-      for(bandIndex <- 0 until tile.bandCount ) yield {
-        Raster(tile.band(bandIndex), self.extent)
-          .reproject(targetRasterExtent, transform, inverseTransform, options)
-      }
+      for (bandIndex <- 0 until tile.bandCount ) yield
+        Raster(tile.band(bandIndex), extent).reproject(targetRasterExtent, transform, inverseTransform, options).tile
 
-    Raster(ArrayMultibandTile(bands.map(_.tile)), bands.head.extent)
+    Raster(ArrayMultibandTile(bands.toArray), targetRasterExtent.extent)
   }
 }
diff --git a/raster/src/main/scala/geotrellis/raster/reproject/MultibandRasterizeReprojectMethods.scala b/raster/src/main/scala/geotrellis/raster/reproject/MultibandRasterizeReprojectMethods.scala
diff --git a/raster/src/main/scala/geotrellis/raster/reproject/ProjectedRasterReprojectMethods.scala b/raster/src/main/scala/geotrellis/raster/reproject/ProjectedRasterReprojectMethods.scala
@@ -25,8 +25,9 @@ import geotrellis.util.MethodExtensions
 import spire.syntax.cfor._
 
 
-class ProjectedRasterReprojectMethods[T <: CellGrid](val self: ProjectedRaster[T])(implicit ev: Raster[T] => RasterReprojectMethods[Raster[T]])
-    extends MethodExtensions[ProjectedRaster[T]] {
+class ProjectedRasterReprojectMethods[T <: CellGrid](val self: ProjectedRaster[T])(
+  implicit ev: Raster[T] => RasterReprojectMethods[Raster[T]]
+) extends MethodExtensions[ProjectedRaster[T]] {
   import Reproject.Options
 
   def reproject(dest: CRS, options: Options): ProjectedRaster[T] =
@@ -41,4 +42,15 @@ class ProjectedRasterReprojectMethods[T <: CellGrid](val self: ProjectedRaster[T
 
   def reproject(gridBounds: GridBounds, dest: CRS): ProjectedRaster[T] =
     reproject(gridBounds, dest, Options.DEFAULT)
+
+  def reproject(dest: CRS, rasterExtent: RasterExtent, resampleMethod: ResampleMethod)
+               (implicit ev: RasterRegionReproject[T]): ProjectedRaster[T] = {
+    ProjectedRaster(
+      raster = ev.regionReproject(
+        self.raster, self.crs, dest, rasterExtent,
+        self.projectedExtent.reprojectAsPolygon(dest, 0.005),
+        resampleMethod),
+      crs = dest
+    )
+  }
 }
diff --git a/raster/src/main/scala/geotrellis/raster/reproject/RasterRegionReproject.scala b/raster/src/main/scala/geotrellis/raster/reproject/RasterRegionReproject.scala
@@ -0,0 +1,85 @@
+package geotrellis.raster.reproject
+
+import geotrellis.raster._
+import geotrellis.raster.resample._
+import geotrellis.raster.rasterize._
+import geotrellis.vector.Polygon
+import geotrellis.proj4._
+
+import spire.syntax.cfor._
+
+trait RasterRegionReproject[T <: CellGrid]{
+  /** Reproject raster to a region that may partially intersect target raster extent.
+    *
+    * Back-projects only the cells overlapping the destination region before sampling their value from the source raster.
+    * This can be used to avoid assigning destination cell values where there is insufficient data for resample kernel.
+    *
+    * @param src source raster CRS
+    * @param dest target raster CRS
+    * @param rasterExtent extent and resolution of target raster
+    * @param region polygon boundry of source raster in target CRS
+    * @param resampleMethod cell value resample method
+    */
+  def regionReproject(raster: Raster[T], src: CRS, dest: CRS, rasterExtent: RasterExtent, region: Polygon, resampleMethod: ResampleMethod): Raster[T]
+}
+
+object RasterRegionReproject {
+  implicit val singlebandInstance = new RasterRegionReproject[Tile] {
+    def regionReproject(raster: Raster[Tile], src: CRS, dest: CRS, rasterExtent: RasterExtent, region: Polygon, resampleMethod: ResampleMethod): Raster[Tile] = {
+      val buffer = raster.tile.prototype(rasterExtent.cols, rasterExtent.rows)
+      val trans = Proj4Transform(dest, src)
+      val resampler = Resample.apply(resampleMethod, raster.tile, raster.extent, CellSize(raster.rasterExtent.cellwidth, raster.rasterExtent.cellheight))
+
+      if (raster.cellType.isFloatingPoint) {
+        Rasterizer.foreachCellByPolygon(region, rasterExtent) { (px, py) =>
+          val (x, y) = rasterExtent.gridToMap(px, py)
+          val (tx, ty) = trans(x, y)
+          buffer.setDouble(px, py, resampler.resampleDouble(tx, ty))
+        }
+      } else {
+        Rasterizer.foreachCellByPolygon(region, rasterExtent) { (px, py) =>
+          val (x, y) = rasterExtent.gridToMap(px, py)
+          val (tx, ty) = trans(x, y)
+          buffer.set(px, py, resampler.resample(tx, ty))
+        }
+      }
+
+      Raster(buffer, rasterExtent.extent)
+    }
+  }
+
+  implicit val multibandInstance = new RasterRegionReproject[MultibandTile] {
+    def regionReproject(raster: Raster[MultibandTile], src: CRS, dest: CRS, rasterExtent: RasterExtent, region: Polygon, resampleMethod: ResampleMethod): Raster[MultibandTile] = {
+      val trans = Proj4Transform(dest, src)
+      val bands = Array.ofDim[MutableArrayTile](raster.tile.bandCount)
+
+      cfor(0)(_ < bands.length, _ + 1) { i =>
+        bands(i) = raster.band(i).prototype(rasterExtent.cols, rasterExtent.rows).mutable
+      }
+
+      val resampler = (0 until raster.bandCount).map { i =>
+        Resample(resampleMethod, raster.band(i), raster.extent, raster.rasterExtent.cellSize)
+      }
+
+      if (raster.cellType.isFloatingPoint) {
+        Rasterizer.foreachCellByPolygon(region, rasterExtent) { (px, py) =>
+          val (x, y) = rasterExtent.gridToMap(px, py)
+          val (tx, ty) = trans(x, y)
+          cfor(0)(_ < bands.length, _ + 1) { i =>
+            bands(i).setDouble(px, py, resampler(i).resampleDouble(tx, ty))
+          }
+        }
+      } else {
+        Rasterizer.foreachCellByPolygon(region, rasterExtent) { (px, py) =>
+          val (x, y) = rasterExtent.gridToMap(px, py)
+          val (tx, ty) = trans(x, y)
+          cfor(0)(_ < bands.length, _ + 1) { i =>
+            bands(i).set(px, py, resampler(i).resample(tx, ty))
+          }
+        }
+      }
+
+      Raster(MultibandTile(bands), rasterExtent.extent)
+    }
+  }
+}
diff --git a/raster/src/main/scala/geotrellis/raster/reproject/RasterReprojectMethods.scala b/raster/src/main/scala/geotrellis/raster/reproject/RasterReprojectMethods.scala
@@ -67,5 +67,4 @@ trait RasterReprojectMethods[+T <: Raster[_]] extends MethodExtensions[T] {
 
   def reproject(gridBounds: GridBounds, transform: Transform, inverseTransform: Transform): T =
     reproject(gridBounds, transform, inverseTransform, Options.DEFAULT)
-
 }
diff --git a/raster/src/main/scala/geotrellis/raster/reproject/RasterizeReprojectMethods.scala b/raster/src/main/scala/geotrellis/raster/reproject/RasterizeReprojectMethods.scala
diff --git a/raster/src/main/scala/geotrellis/raster/reproject/Reproject.scala b/raster/src/main/scala/geotrellis/raster/reproject/Reproject.scala
@@ -13,7 +13,6 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-
 package geotrellis.raster.reproject
 
 import geotrellis.raster._

diff --git a/raster/src/main/scala/geotrellis/raster/reproject/SinglebandRasterReprojectMethods.scala b/raster/src/main/scala/geotrellis/raster/reproject/SinglebandRasterReprojectMethods.scala
@@ -18,8 +18,9 @@ package geotrellis.raster.reproject
 
 import geotrellis.raster._
 import geotrellis.raster.resample._
-import geotrellis.vector.Extent
+import geotrellis.raster.rasterize._
 import geotrellis.proj4._
+import geotrellis.vector.Polygon
 
 import spire.syntax.cfor._
 
@@ -33,7 +34,6 @@ trait SinglebandRasterReprojectMethods extends RasterReprojectMethods[Singleband
     options: Options
   ): SinglebandRaster = {
     val Raster(tile, extent) = self
-    val RasterExtent(_, cellwidth, cellheight, _, _) = self.rasterExtent
     val RasterExtent(newExtent, newCellWidth, newCellHeight, newCols, newRows) = targetRasterExtent
 
     val newTile = ArrayTile.empty(tile.cellType, newCols, newRows)
@@ -63,28 +63,27 @@ trait SinglebandRasterReprojectMethods extends RasterReprojectMethods[Singleband
     val resampler = Resample(options.method, tile, extent, CellSize(newCellWidth, newCellHeight))
 
     if(tile.cellType.isFloatingPoint) {
-      val resample = resampler.resampleDouble _
       cfor(0)(_ < newRows, _ + 1) { row =>
         // Reproject this whole row.
         rowTransform(destX, destY, srcX, srcY)
         cfor(0)(_ < newCols, _ + 1) { col =>
-          val v = resample(srcX(col), srcY(col))
+          val x = srcX(col)
+          val y = srcY(col)
+          val v = resampler.resampleDouble(x, y)
           newTile.setDouble(col, row, v)
 
           // Add row height for next iteration
           destY(col) -= newCellHeight
         }
       }
     } else {
-      val resample = resampler.resample _
       cfor(0)(_ < newRows, _ + 1) { row =>
         // Reproject this whole row.
         rowTransform(destX, destY, srcX, srcY)
         cfor(0)(_ < newCols, _ + 1) { col =>
           val x = srcX(col)
           val y = srcY(col)
-
-          val v = resample(x, y)
+          val v = resampler.resample(x, y)
           newTile.set(col, row, v)
 
           // Add row height for next iteration