diff --git a/raster-test/src/test/scala/geotrellis/raster/viewshed/R2ViewshedSpec.scala b/raster-test/src/test/scala/geotrellis/raster/viewshed/R2ViewshedSpec.scala index b8fb80689a..457aa947ec 100644 --- a/raster-test/src/test/scala/geotrellis/raster/viewshed/R2ViewshedSpec.scala +++ b/raster-test/src/test/scala/geotrellis/raster/viewshed/R2ViewshedSpec.scala @@ -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 ) @@ -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, @@ -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, @@ -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 ) @@ -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 @@ -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 diff --git a/raster/src/main/scala/geotrellis/raster/viewshed/R2Viewshed.scala b/raster/src/main/scala/geotrellis/raster/viewshed/R2Viewshed.scala index a6201183fb..3b8e71a62c 100644 --- a/raster/src/main/scala/geotrellis/raster/viewshed/R2Viewshed.scala +++ b/raster/src/main/scala/geotrellis/raster/viewshed/R2Viewshed.scala @@ -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 @@ -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, @@ -98,7 +103,7 @@ object R2Viewshed extends Serializable { FromInside(), null, { (_, _) => }, - false, // OR + and, false // Ignore curvature ) viewshedTile @@ -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) } } } diff --git a/spark/src/main/scala/geotrellis/spark/viewshed/IterativeViewshed.scala b/spark/src/main/scala/geotrellis/spark/viewshed/IterativeViewshed.scala index 6485d956b9..1e08a291b7 100644 --- a/spark/src/main/scala/geotrellis/spark/viewshed/IterativeViewshed.scala +++ b/spark/src/main/scala/geotrellis/spark/viewshed/IterativeViewshed.scala @@ -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) } diff --git a/spark/src/test/scala/geotrellis/spark/viewshed/IterativeViewshedSpec.scala b/spark/src/test/scala/geotrellis/spark/viewshed/IterativeViewshedSpec.scala index 1573628dc6..5ce32737c5 100644 --- a/spark/src/test/scala/geotrellis/spark/viewshed/IterativeViewshedSpec.scala +++ b/spark/src/test/scala/geotrellis/spark/viewshed/IterativeViewshedSpec.scala @@ -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) @@ -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) @@ -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