Skip to content

Commit

Permalink
Match R2 AND, OR Behavior To Other Viewshed Algs
Browse files Browse the repository at this point in the history
  • Loading branch information
James McClain committed Jun 14, 2017
1 parent 7259308 commit 64ba090
Show file tree
Hide file tree
Showing 4 changed files with 144 additions and 20 deletions.
132 changes: 124 additions & 8 deletions raster-test/src/test/scala/geotrellis/raster/viewshed/R2ViewshedSpec.scala
View file
Expand Up @@ -153,19 +153,31 @@ class R2ViewshedSpec extends FunSpec
all should be (20)
}

it("computes the viewshed of a flat int plane") {
it("computes the viewshed of a flat int plane (OR)") {
val r = createTile(Array.fill(7 * 8)(1), 7, 8)
val shed = R2Viewshed(r, 4, 3)
assertEqual(BitConstantTile(true, 7, 8), shed)
}

it("computes the viewshed of a flat double plane") {
it("computes the viewshed of a flat int plane (AND)") {
val r = createTile(Array.fill(7 * 8)(1), 7, 8)
val shed = R2Viewshed(r, 4, 3, true)
assertEqual(BitConstantTile(true, 7, 8), shed)
}

it("computes the viewshed of a flat double plane (OR)") {
val r = createTile(Array.fill(7 * 8)(1.5), 7, 8)
val shed = R2Viewshed(r, 4, 3)
assertEqual(BitConstantTile(true, 7, 8), shed)
}

it("computes the viewshed of a double line") {
it("computes the viewshed of a flat double plane (AND)") {
val r = createTile(Array.fill(7 * 8)(1.5), 7, 8)
val shed = R2Viewshed(r, 4, 3, true)
assertEqual(BitConstantTile(true, 7, 8), shed)
}

it("computes the viewshed of a double line (OR)") {
val rasterData = Array (
300.0, 1.0, 99.0, 0.0, 10.0, 200.0, 137.0
)
Expand All @@ -178,7 +190,20 @@ class R2ViewshedSpec extends FunSpec
assertEqual(viewRaster, shed)
}

it("computes the viewshed of a double plane") {
it("computes the viewshed of a double line (AND)") {
val rasterData = Array (
300.0, 1.0, 99.0, 0.0, 10.0, 200.0, 137.0
)
val viewable = Array (
1, 0, 1, 1, 1, 1, 0
)
val r = createTile(rasterData, 7, 1)
val viewRaster = createTile(viewable, 7, 1).convert(BitCellType)
val shed = R2Viewshed(r, 3, 0, true)
assertEqual(viewRaster, shed)
}

it("computes the viewshed of a double plane (OR)") {
val rasterData = Array (
999.0, 1.0, 1.0, 1.0, 1.0, 1.0, 999.0,
1.0, 1.0, 1.0, 1.0, 1.0, 499.0, 1.0,
Expand All @@ -203,7 +228,32 @@ class R2ViewshedSpec extends FunSpec
assertEqual(viewRaster, shed)
}

it("computes the viewshed of a int plane") {
it("computes the viewshed of a double plane (AND)") {
val rasterData = Array (
999.0, 1.0, 1.0, 1.0, 1.0, 1.0, 999.0,
1.0, 1.0, 1.0, 1.0, 1.0, 499.0, 1.0,
1.0, 1.0, 1.0, 1.0, 99.0, 1.0, 1.0,
1.0, 1.0, 999.0, 1.0, 1.0, 1.0, 1.0,
1.0, 1.0, 1.0, 1.0, 100.0, 1.0, 1.0,
1.0, 1.0, 1.0, 1.0, 1.0, 101.0, 1.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 102.0
)
val viewable = Array (
1, 1, 1, 1, 1, 0, 1,
1, 1, 1, 1, 0, 1, 0,
0, 0, 1, 1, 1, 0, 1,
0, 0, 1, 1, 1, 1, 1,
0, 0, 1, 1, 1, 0, 1,
1, 1, 1, 1, 0, 0, 0,
1, 1, 1, 1, 1, 0, 0
)
val r = createTile(rasterData, 7, 7)
val viewRaster = createTile(viewable, 7, 7).convert(BitCellType)
val shed = R2Viewshed(r, 3, 3, true)
assertEqual(viewRaster, shed)
}

it("computes the viewshed of a int plane (OR)") {
val rasterData = Array (
999, 1, 1, 1, 1, 499, 999,
1, 1, 1, 1, 1, 499, 250,
Expand All @@ -228,7 +278,32 @@ class R2ViewshedSpec extends FunSpec
assertEqual(viewRaster, shed)
}

it("ignores NoData values and indicates they're unviewable"){
it("computes the viewshed of a int plane (AND)") {
val rasterData = Array (
999, 1, 1, 1, 1, 499, 999,
1, 1, 1, 1, 1, 499, 250,
1, 1, 1, 1, 99, 1, 1,
1, 999, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 0, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1
)
val viewable = Array (
1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 0, 1, 0,
1, 1, 1, 1, 1, 0, 1,
0, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 0, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1
)
val r = createTile(rasterData, 7, 7)
val viewRaster = createTile(viewable, 7, 7).convert(BitCellType)
val shed = R2Viewshed(r, 3, 3, true)
assertEqual(viewRaster, shed)
}

it("ignores NoData values and indicates they're unviewable (OR)"){
val rasterData = Array (
300.0, 1.0, 99.0, 0.0, Double.NaN, 200.0, 137.0
)
Expand All @@ -241,7 +316,20 @@ class R2ViewshedSpec extends FunSpec
assertEqual(viewRaster, shed)
}

it("should match veiwshed generated by ArgGIS") {
it("ignores NoData values and indicates they're unviewable (AND)"){
val rasterData = Array (
300.0, 1.0, 99.0, 0.0, Double.NaN, 200.0, 137.0
)
val viewable = Array (
1, 0, 1, 1, 0, 1, 0
)
val r = createTile(rasterData, 7, 1)
val viewRaster = createTile(viewable, 7, 1).convert(BitCellType)
val shed = R2Viewshed(r, 3, 0, true)
assertEqual(viewRaster, shed)
}

it("should match veiwshed generated by ArgGIS (OR)") {
val rs = loadTestArg("data/viewshed-elevation")
val elevation = rs.tile
val rasterExtent = rs.rasterExtent
Expand All @@ -262,7 +350,35 @@ class R2ViewshedSpec extends FunSpec
}

val diff = (countDiff(expected, actual) / (256 * 256).toDouble) * 100
val allowable = 8.75
val allowable = 8.72
// System.out.println(s"${diff} / ${256 * 256} = ${diff / (256 * 256).toDouble}")
withClue(s"Percent difference from ArgGIS raster is more than $allowable%:") {
diff should be < allowable
}
}

it("should match veiwshed generated by ArgGIS (AND)") {
val rs = loadTestArg("data/viewshed-elevation")
val elevation = rs.tile
val rasterExtent = rs.rasterExtent
val expected = loadTestArg("data/viewshed-expected")

val (x, y) = (-93.63300872055451407, 30.54649743277299123) // create overload
val (col, row) = rasterExtent.mapToGrid(x, y)
val actual = R2Viewshed(elevation, col, row, true)

def countDiff(a: Tile, b: Tile): Int = {
var ans = 0
for(col <- 0 until 256) {
for(row <- 0 until 256) {
if (a.get(col, row) != b.get(col, row)) ans += 1;
}
}
ans
}

val diff = (countDiff(expected, actual) / (256 * 256).toDouble) * 100
val allowable = 9.01
// System.out.println(s"${diff} / ${256 * 256} = ${diff / (256 * 256).toDouble}")
withClue(s"Percent difference from ArgGIS raster is more than $allowable%:") {
diff should be < allowable
Expand Down
17 changes: 12 additions & 5 deletions raster/src/main/scala/geotrellis/raster/viewshed/R2Viewshed.scala
View file
Expand Up @@ -66,7 +66,7 @@ object R2Viewshed extends Serializable {
* @param rows The number of rows
*/
def generateEmptyViewshedTile(cols: Int, rows: Int) =
ArrayTile.empty(IntCellType, cols, rows)
ArrayTile.empty(IntConstantNoDataCellType, cols, rows)

/**
* Compute the drop in elevation due to Earth's curvature (please
Expand All @@ -85,11 +85,16 @@ object R2Viewshed extends Serializable {
* @param startCol The x position of the vantage point
* @param startRow The y position of the vantage point
*/
def apply(elevationTile: Tile, startCol: Int, startRow: Int): Tile = {
def apply(
elevationTile: Tile,
startCol: Int, startRow: Int,
and: Boolean = false): Tile = {
val cols = elevationTile.cols
val rows = elevationTile.rows
val viewHeight = elevationTile.getDouble(startCol, startRow)
val viewshedTile = generateEmptyViewshedTile(cols, rows)
val viewshedTile =
if (!and) ArrayTile.empty(IntCellType, cols, rows)
else ArrayTile.empty(IntConstantNoDataCellType, cols, rows)

R2Viewshed.compute(
elevationTile, viewshedTile,
Expand All @@ -98,7 +103,7 @@ object R2Viewshed extends Serializable {
FromInside(),
null,
{ (_, _) => },
false, // OR
and,
false // Ignore curvature
)
viewshedTile
Expand Down Expand Up @@ -248,10 +253,12 @@ object R2Viewshed extends Serializable {
if (distance >= maxDistance) alpha.terminated = true
if (!alpha.terminated) {
val visible = alpha <= angle
val bit = viewshedTile.get(col, row)

if (visible) alpha.alpha = angle
if (!and && visible) viewshedTile.set(col, row, 1)
else if (and && !visible) viewshedTile.set(col, row, 0)
else if (and && visible && isNoData(viewshedTile.get(col, row))) viewshedTile.set(col, row, 1)
else if (and && visible && isNoData(bit)) viewshedTile.set(col, row, 1)
}
}
}
Expand Down
2 changes: 1 addition & 1 deletion spark/src/main/scala/geotrellis/spark/viewshed/IterativeViewshed.scala
View file
Expand Up @@ -231,7 +231,7 @@ object IterativeViewshed {
} while (rays.value.size > 0)

// Return the computed viewshed layer
val metadata = TileLayerMetadata(DoubleCellType, md.layout, md.extent, md.crs, md.bounds)
val metadata = TileLayerMetadata(IntConstantNoDataCellType, md.layout, md.extent, md.crs, md.bounds)
val rdd = sheds.map({ case (k, _, v) => (k, v.asInstanceOf[Tile]) })
ContextRDD(rdd, metadata)
}
Expand Down
13 changes: 7 additions & 6 deletions spark/src/test/scala/geotrellis/spark/viewshed/IterativeViewshedSpec.scala
View file
Expand Up @@ -45,7 +45,7 @@ class IterativeViewshedSpec extends FunSpec
}

val viewshed = IterativeViewshed(rdd, Point(7, 7), -0.0, Double.PositiveInfinity, false, false)
val actual = viewshed.map({ case (_, v) => v.toArray.sum }).reduce(_ + _)
var actual = 0 ; viewshed.collect.foreach({ case (_, v) => v.foreach({ z => if (isData(z)) actual += z }) })
val expected = 15*15

actual should be (expected)
Expand All @@ -64,7 +64,7 @@ class IterativeViewshedSpec extends FunSpec
}

val viewshed = IterativeViewshed(rdd, Point(7, 7), -0.0, Double.PositiveInfinity, false, false)
val actual = viewshed.map({ case (_, v) => v.toArray.sum }).reduce(_ + _)
var actual = 0 ; viewshed.collect.foreach({ case (_, v) => v.foreach({ z => if (isData(z)) actual += z }) })
val expected = 180

actual should be (expected)
Expand All @@ -84,12 +84,13 @@ class IterativeViewshedSpec extends FunSpec
}

val viewshed = IterativeViewshed(rdd, Point(7, 7), -0.0, Double.PositiveInfinity, false, false)
val ND = NODATA
val expected: Array[Int] = Array(
1, 1, 1, 1, 1,
1, 1, 0, 1, 1,
1, 1, 0, 1, 1,
1, 0, 0, 0, 1,
1, 0, 1, 0, 1
1, 1, ND, 1, 1,
1, 1, ND, 1, 1,
1, ND, ND, ND, 1,
1, ND, 1, ND, 1
)
val actual = viewshed
.collect
Expand Down

0 comments on commit 64ba090

Please sign in to comment.