Skip to content

Commit

Permalink
Add union to merge methods (#3482)
Browse files Browse the repository at this point in the history
* Add union to merge methods

* Change unionF to unionFunc

* Remove code duplication in merge logic

* Clean up semigroup usage

* Update changelog and docstring
  • Loading branch information
moradology authored Oct 6, 2024
1 parent a43cc60 commit 24f6271
Show file tree
Hide file tree
Showing 8 changed files with 227 additions and 12 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
### Added
- Add RasterSourceRDD.tiledLayerRDD within the geometry intersection [#3474](https://github.com/locationtech/geotrellis/pull/3474)
- Expose AWS_REQUEST_PAYER environment variable [#3479](https://github.com/locationtech/geotrellis/pull/3479)
- Add union/outer join merger of tiles [#3482](https://github.com/locationtech/geotrellis/pull/3482)

### Changed
- Migration to CE3 and other major dependencies upgrade [#3389](https://github.com/locationtech/geotrellis/pull/3389)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -63,4 +63,35 @@ trait MultibandTileMergeMethods extends TileMergeMethods[MultibandTile] {

ArrayMultibandTile(bands)
}

/**
* Union this [[MultibandTile]] with the other one. The output tile's
* extent will be the minimal extent which encompasses both input extents.
* A new MutlibandTile is returned.
*
* @param extent The extent of this MultiBandTile
* @param otherExtent The extent of the other MultiBandTile
* @param other The other MultiBandTile
* @param method The resampling method
* @param unionFunc The function which decides how values from rasters being combined will be transformed
* @return A new MultiBandTile, the result of the merge
*/
def union(
extent: Extent,
otherExtent: Extent,
other: MultibandTile,
method: ResampleMethod,
unionFunc: (Option[Double], Option[Double]) => Double
): MultibandTile = {
val bands: Seq[Tile] =
for {
bandIndex <- 0 until self.bandCount
} yield {
val thisBand = self.band(bandIndex)
val thatBand = other.band(bandIndex)
thisBand.union(extent, otherExtent, thatBand, method, unionFunc)
}

ArrayMultibandTile(bands)
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -50,4 +50,16 @@ abstract class RasterMergeMethods[
*/
def merge(other: Raster[T]): Raster[T] =
merge(other, NearestNeighbor)

/**
* Union this [[Raster]] with the other one. All places in the
* present raster that contain NODATA are filled-in with data from
* the other raster. A new Raster is returned.
*
* @param other The other Raster
* @param method The resampling method
* @return A new Raster, the result of the merge
*/
def union(other: Raster[T], method: ResampleMethod, unionFunc: (Option[Double], Option[Double]) => Double): Raster[T] =
Raster(self.tile.union(self.extent, other.extent, other.tile, method, unionFunc), self.extent.combine(other.extent))
}
Original file line number Diff line number Diff line change
Expand Up @@ -32,4 +32,7 @@ abstract class RasterTileFeatureMergeMethods[

def merge(other: TileFeature[Raster[T], D], method: ResampleMethod): TileFeature[Raster[T], D] =
TileFeature(self.tile.merge(other.tile, method), Semigroup[D].combine(self.data, other.data))

def union(other: TileFeature[Raster[T], D], method: ResampleMethod, unionFunc: (Option[Double], Option[Double]) => Double): TileFeature[Raster[T], D] =
TileFeature(self.tile.union(other.tile, method, unionFunc), Semigroup[D].combine(self.data, other.data))
}
Original file line number Diff line number Diff line change
Expand Up @@ -99,52 +99,48 @@ trait SinglebandTileMergeMethods extends TileMergeMethods[Tile] {
val gridBounds = re.gridBoundsFor(sharedExtent)
val targetCS = CellSize(sharedExtent, gridBounds.width, gridBounds.height)

val interpolate = Resample(method, other, otherExtent, targetCS)

self.cellType match {
case BitCellType | ByteCellType | UByteCellType | ShortCellType | UShortCellType | IntCellType =>
val interpolate: (Double, Double) => Int = Resample(method, other, otherExtent, targetCS).resample _
// Assume 0 as the transparent value
cfor(0)(_ < self.rows, _ + 1) { row =>
cfor(0)(_ < self.cols, _ + 1) { col =>
if (self.get(col, row) == 0) {
val (x, y) = re.gridToMap(col, row)
val v = interpolate(x, y)
val v = interpolate.resample(x, y)
if(isData(v)) {
mutableTile.set(col, row, v)
}
}
}
}
case FloatCellType | DoubleCellType =>
val interpolate: (Double, Double) => Double = Resample(method, other, otherExtent, targetCS).resampleDouble _
// Assume 0.0 as the transparent value
cfor(0)(_ < self.rows, _ + 1) { row =>
cfor(0)(_ < self.cols, _ + 1) { col =>
if (self.getDouble(col, row) == 0.0) {
val (x, y) = re.gridToMap(col, row)
val v = interpolate(x, y)
val v = interpolate.resampleDouble(x, y)
if(isData(v)) {
mutableTile.setDouble(col, row, v)
}
}
}
}
case x if x.isFloatingPoint =>
val interpolate: (Double, Double) => Double = Resample(method, other, otherExtent, targetCS).resampleDouble _
cfor(0)(_ < self.rows, _ + 1) { row =>
cfor(0)(_ < self.cols, _ + 1) { col =>
if (isNoData(self.getDouble(col, row))) {
val (x, y) = re.gridToMap(col, row)
mutableTile.setDouble(col, row, interpolate(x, y))
mutableTile.setDouble(col, row, interpolate.resampleDouble(x, y))
}
}
}
case _ =>
val interpolate: (Double, Double) => Int = Resample(method, other, otherExtent, targetCS).resample _
cfor(0)(_ < self.rows, _ + 1) { row =>
cfor(0)(_ < self.cols, _ + 1) { col =>
if (isNoData(self.get(col, row))) {
val (x, y) = re.gridToMap(col, row)
mutableTile.set(col, row, interpolate(x, y))
mutableTile.set(col, row, interpolate.resample(x, y))
}
}
}
Expand All @@ -154,4 +150,80 @@ trait SinglebandTileMergeMethods extends TileMergeMethods[Tile] {
case _ =>
self
}

/** Unions this tile with another tile, preserving the non-overlapping
* portions of each tile extent.
*
* This method requires a union function to be provided which decides
* how values will be added to the output raster. A `None` represents
* one of the two input rasters not having a value at the resampled location.
* Two `None` values mean that both tiles are missing values at the
* cell in question.
*/
def union(extent: Extent, otherExtent: Extent, other: Tile, method: ResampleMethod, unionFunc: (Option[Double], Option[Double]) => Double): Tile = {
val unionInt: (Option[Int], Option[Int]) => Int =
(l: Option[Int], r: Option[Int]) => {
unionFunc(l.map(_.toDouble), r.map(_.toDouble)).toInt
}

val combinedExtent = otherExtent combine extent
val re = RasterExtent(extent, self.cols, self.rows)
val gridBounds = re.gridBoundsFor(combinedExtent, false)
val targetCS = CellSize(combinedExtent, gridBounds.width, gridBounds.height)
val targetRE = RasterExtent(combinedExtent, targetCS)
val mutableTile = ArrayTile.empty(self.cellType, targetRE.cols, targetRE.rows)

val interpolateLeft = Resample(method, self, extent, targetCS)
val interpolateRight = Resample(method, other, otherExtent, targetCS)

self.cellType match {
case BitCellType | ByteCellType | UByteCellType | ShortCellType | UShortCellType | IntCellType =>
cfor(0)(_ < targetRE.rows, _ + 1) { row =>
cfor(0)(_ < targetRE.cols, _ + 1) { col =>
val (x,y) = targetRE.gridToMap(col, row)
val l = interpolateLeft.resample(x, y)
val r = interpolateRight.resample(x, y)
val maybeL = Some(l)
val maybeR = Some(r)
val v = unionInt(maybeL, maybeR)
mutableTile.set(col, row, v)
}
}
case FloatCellType | DoubleCellType =>
cfor(0)(_ < targetRE.rows, _ + 1) { row =>
cfor(0)(_ < targetRE.cols, _ + 1) { col =>
val (x,y) = targetRE.gridToMap(col, row)
val l = interpolateLeft.resampleDouble(x, y)
val r = interpolateRight.resampleDouble(x, y)
val maybeL = Some(l)
val maybeR = Some(r)
val v = unionFunc(maybeL, maybeR)
mutableTile.setDouble(col, row, v)
}
}
case x if x.isFloatingPoint =>
cfor(0)(_ < targetRE.rows, _ + 1) { row =>
cfor(0)(_ < targetRE.cols, _ + 1) { col =>
val (x,y) = targetRE.gridToMap(col, row)
val l = interpolateLeft.resampleDouble(x, y)
val r = interpolateRight.resampleDouble(x, y)
val maybeL = if (isNoData(l)) None else Some(l)
val maybeR = if (isNoData(r)) None else Some(r)
mutableTile.setDouble(col, row, unionFunc(maybeL, maybeR))
}
}
case _ =>
cfor(0)(_ < targetRE.rows, _ + 1) { row =>
cfor(0)(_ < targetRE.cols, _ + 1) { col =>
val (x,y) = targetRE.gridToMap(col, row)
val l = interpolateLeft.resample(x, y)
val r = interpolateRight.resample(x, y)
val maybeL = if (isNoData(l)) None else Some(l)
val maybeR = if (isNoData(r)) None else Some(r)
mutableTile.set(col, row, unionInt(maybeL, maybeR))
}
}
}
mutableTile
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -19,16 +19,27 @@ package geotrellis.raster.merge
import geotrellis.raster._
import geotrellis.raster.resample._
import geotrellis.vector._

import cats.Semigroup
import cats.syntax.semigroup._


abstract class TileFeatureMergeMethods[
T <: CellGrid[Int]: * => TileMergeMethods[T],
D: Semigroup
](val self: TileFeature[T, D]) extends TileMergeMethods[TileFeature[T, D]] {
def merge(other: TileFeature[T, D], baseCol: Int, baseRow: Int): TileFeature[T, D] =
TileFeature(self.tile.merge(other.tile, baseCol, baseRow), Semigroup[D].combine(self.data, other.data))
TileFeature(self.tile.merge(other.tile, baseCol, baseRow), self.data combine other.data)

def merge(extent: Extent, otherExtent: Extent, other: TileFeature[T, D], method: ResampleMethod): TileFeature[T, D] =
TileFeature(self.tile.merge(extent, otherExtent, other.tile, method), Semigroup[D].combine(self.data, other.data))
TileFeature(self.tile.merge(extent, otherExtent, other.tile, method), self.data combine other.data)

def union(
extent: Extent,
otherExtent: Extent,
other: TileFeature[T, D],
method: ResampleMethod,
unionFunc: (Option[Double], Option[Double]) => Double
): TileFeature[T, D] =
TileFeature(self.tile.union(extent, otherExtent, other.tile, method, unionFunc), self.data combine other.data)
}
Original file line number Diff line number Diff line change
Expand Up @@ -78,4 +78,9 @@ trait TileMergeMethods[T] extends MethodExtensions[T] {
def merge(extent: Extent, otherExtent: Extent, other: T): T =
merge(extent, otherExtent, other, NearestNeighbor)


def union(extent: Extent, otherExtent: Extent, other: T, method: ResampleMethod, unionFunc: (Option[Double], Option[Double]) => Double): T

def union(extent: Extent, otherExtent: Extent, other: T, unionFunc: (Option[Double], Option[Double]) => Double): T =
union(extent, otherExtent, other, NearestNeighbor, unionFunc)
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
/*
* Copyright 2016 Azavea
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/

package geotrellis.raster.merge

import geotrellis.raster._
import geotrellis.raster.testkit._
import geotrellis.raster.resample.NearestNeighbor
import geotrellis.vector.Extent

import org.scalatest.matchers.should.Matchers
import org.scalatest.funspec.AnyFunSpec

class TileUnionMethodsSpec extends AnyFunSpec
with Matchers
with TileBuilders
with RasterMatchers {
describe("SinglebandTileMergeMethods") {

it("should union two tiles such that the extent of the output is equal to the minimum extent which covers both") {
val cellTypes: Seq[CellType] =
Seq(
BitCellType,
ByteCellType,
ByteConstantNoDataCellType,
ByteUserDefinedNoDataCellType(1.toByte),
UByteCellType,
UByteConstantNoDataCellType,
UByteUserDefinedNoDataCellType(1.toByte),
ShortCellType,
ShortConstantNoDataCellType,
ShortUserDefinedNoDataCellType(1.toShort),
UShortCellType,
UShortConstantNoDataCellType,
UShortUserDefinedNoDataCellType(1.toShort),
IntCellType,
IntConstantNoDataCellType,
IntUserDefinedNoDataCellType(1),
FloatCellType,
FloatConstantNoDataCellType,
FloatUserDefinedNoDataCellType(1.0f),
DoubleCellType,
DoubleConstantNoDataCellType,
DoubleUserDefinedNoDataCellType(1.0)
)

for(ct <- cellTypes) {
val arr = Array.ofDim[Double](100).fill(5.0)
arr(50) = 1.0
arr(55) = 0.0
arr(60) = Double.NaN

val tile1 =
DoubleArrayTile(arr, 10, 10).convert(ct)
val e1 = Extent(0, 0, 1, 1)
val tile2 =
tile1.prototype(ct, tile1.cols, tile1.rows)
val e2 = Extent(1, 1, 2, 2)
val unioned = tile1.union(e1, e2, tile2, NearestNeighbor, (d1, d2) => d1.getOrElse(4))
withClue(s"Failing on cell type $ct: ") {
unioned.rows shouldBe (20)
unioned.cols shouldBe (20)
}
}
}
}
}

0 comments on commit 24f6271

Please sign in to comment.