Permalink
Browse files

Fix reprojection bug in pyramid

Signed-off-by: jpolchlo <jpolchlopek@azavea.com>
  • Loading branch information...
jpolchlo authored and echeipesh committed Apr 24, 2018
1 parent e3b4fca commit 2eefb6c55bc6e6cf23ee0077bce98f2611da1fd0
Showing with 31 additions and 21 deletions.
  1. +31 −21 spark/src/main/scala/geotrellis/spark/pyramid/Pyramid.scala
@@ -44,9 +44,13 @@ object Pyramid extends LazyLogging {
/** Resample base layer to generate the next level "up" in the pyramid.
*
* When building the next pyramid level each tile is resampled individually, without sampling pixels from neighboring tiles.
* In the common case of [[ZoomedLayoutScheme]], used when building a power of two pyramid, this means that
* each resampled pixel is going to be equidistant to four pixels in the base layer.
* Builds the pyramid level with a cell size twice that of the input level---the "next level up" in the pyramid.
* Each tile is resampled individually, without sampling pixels from neighboring tiles to speed up the process.
* The algorithm proceeds by reducing the input tiles by half using a resampling method over 2x2 pixel neighborhoods.
* We support all [[AggregateResampleMethod]]s as well as NearestNeighbor and Bilinear resampling. Nearest neighbor
* resampling is, strictly speaking, non-deterministic in this setting, but is included to support categorical layers
* (e.g., NLCD). Given this method, input tile layers are obviously expected to comprise tiles with even pixel
* dimensions.
*
* @param rdd the base layer to be resampled
* @param layoutScheme the scheme used to generate next pyramid level
@@ -58,7 +62,7 @@ object Pyramid extends LazyLogging {
*/
def up[
K: SpatialComponent: ClassTag,
V <: CellGrid: ClassTag: ? => TileMergeMethods[V]: ? => TilePrototypeMethods[V],
V <: CellGrid: ClassTag: ? => TilePrototypeMethods[V]: ? => TileMergeMethods[V],
M: Component[?, LayoutDefinition]: Component[?, Bounds[K]]
](rdd: RDD[(K, V)] with Metadata[M],
layoutScheme: LayoutScheme,
@@ -67,10 +71,16 @@ object Pyramid extends LazyLogging {
): (Int, RDD[(K, V)] with Metadata[M]) = {
val Options(resampleMethod, partitioner) = options
assert(! Seq(CubicConvolution, CubicSpline, Lanczos).map(_ canEqual resampleMethod).reduce(_||_),
s"${resampleMethod} resample method is not supported for pyramid construction")
val sourceLayout = rdd.metadata.getComponent[LayoutDefinition]
val sourceBounds = rdd.metadata.getComponent[Bounds[K]]
val LayoutLevel(nextZoom, nextLayout) = layoutScheme.zoomOut(LayoutLevel(zoom, sourceLayout))
assert(sourceLayout.tileCols % 2 == 0 && sourceLayout.tileRows % 2 == 0,
s"Pyramid operation requires tiles with even dimensions, got ${sourceLayout.tileCols} x ${sourceLayout.tileRows}")
val nextKeyBounds =
sourceBounds match {
case EmptyBounds => EmptyBounds
@@ -103,35 +113,35 @@ object Pyramid extends LazyLogging {
.setComponent(nextKeyBounds)
// Functions for combine step
def createTiles(tile: (K, V)): Seq[(K, V)] = Seq(tile)
def mergeTiles1(tiles: Seq[(K, V)], tile: (K, V)): Seq[(K, V)] = tiles :+ tile
def mergeTiles2(tiles1: Seq[(K, V)], tiles2: Seq[(K, V)]): Seq[(K, V)] = tiles1 ++ tiles2
def createTiles(tile: Raster[V]): Seq[Raster[V]] = Seq(tile)
def mergeTiles1(tiles: Seq[Raster[V]], tile: Raster[V]): Seq[Raster[V]] = tiles :+ tile
def mergeTiles2(tiles1: Seq[Raster[V]], tiles2: Seq[Raster[V]]): Seq[Raster[V]] = tiles1 ++ tiles2
val nextRdd = {
val transformedRdd = rdd
val transformedRdd = rdd
.map { case (key, tile) =>
val extent: Extent = key.getComponent[SpatialKey].extent(sourceLayout)
val newSpatialKey = nextLayout.mapTransform(extent.center)
// Resample the tile on the map side of the pyramid step.
// This helps with shuffle size.
val resampled = tile.prototype(nextLayout.tileLayout.tileCols, nextLayout.tileLayout.tileRows)
val resampled = tile.prototype(sourceLayout.tileCols / 2, sourceLayout.tileRows / 2)
resampled.merge(extent, extent, tile, resampleMethod)
(key.setComponent(newSpatialKey), (key, resampled))
(key.setComponent(newSpatialKey), Raster(resampled, extent))
}
partitioner
.fold(transformedRdd.combineByKey(createTiles, mergeTiles1, mergeTiles2))(transformedRdd.combineByKey(createTiles _, mergeTiles1 _, mergeTiles2 _, _))
.mapPartitions ( partition => partition.map { case (newKey: K, seq: Seq[(K, V)]) =>
val newExtent = newKey.getComponent[SpatialKey].extent(nextLayout)
val newTile = seq.head._2.prototype(nextLayout.tileLayout.tileCols, nextLayout.tileLayout.tileRows)
partitioner
.fold(transformedRdd.combineByKey(createTiles, mergeTiles1, mergeTiles2))(transformedRdd.combineByKey(createTiles _, mergeTiles1 _, mergeTiles2 _, _))
.mapPartitions ( partition => partition.map { case (newKey: K, seq: Seq[Raster[V]]) =>
val newExtent = newKey.getComponent[SpatialKey].extent(nextLayout)
val newTile = seq.head.tile.prototype(nextLayout.tileLayout.tileCols, nextLayout.tileLayout.tileRows)
for ((oldKey, tile) <- seq) {
val oldExtent = oldKey.getComponent[SpatialKey].extent(sourceLayout)
newTile.merge(newExtent, oldExtent, tile, NearestNeighbor)
}
(newKey, newTile: V)
}, preservesPartitioning = true)
for (raster <- seq) {
newTile.merge(newExtent, raster.extent, raster.tile, NearestNeighbor)
}
(newKey, newTile: V)
}, preservesPartitioning = true)
}
nextZoom -> new ContextRDD(nextRdd, nextMetadata)

0 comments on commit 2eefb6c

Please sign in to comment.