updated computeResolution in IterativeCostDistance to match Iterative…

… Viewshed (#3106) * updated computeResolution in IterativeCostDistance to match IterativeViewshed The compute resolution method in Iterative Cost Distance was calcultating the cell size based on the length of a degree at the equator without scaling for the current latitude. This introduces an error as you move away from the equator. The same method in Iterative Viewshed compensates for this and maintains greater spatial accuracy. Signed-off-by: jmtaysom <[email protected]> Signed-off-by: Grigory Pomadchin <[email protected]>
locationtech · Oct 3, 2019 · 82878cb · 82878cb
1 parent dbbeedc
commit 82878cb
Show file tree

Hide file tree

Showing 4 changed files with 27 additions and 45 deletions.
diff --git a/docs/CHANGELOG.rst b/docs/CHANGELOG.rst
@@ -137,6 +137,7 @@ Fixes & Updates
 - AmazonS3URI.getKey() now returns an empty string instead of null if no key was provided (`#3096 https://github.com/locationtech/geotrellis/pull/3096`_).
 - Improved handling of null prefix for S3AttributeStore. Prefer use of `S3AttributeStore.apply` to take advantage of this improved handling (`#3096 https://github.com/locationtech/geotrellis/pull/3096`_).
 - Fix Tiled TIFF BitCellType Segments conversion (`#3102 https://github.com/locationtech/geotrellis/pull/3102`_)
+- Updated computeResolution in IterativeCostDistance to match Iterative Viewshed (`#3106 https://github.com/locationtech/geotrellis/pull/3106`_)
 
 2.3.0
 -----

diff --git a/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala b/spark/src/main/scala/geotrellis/spark/costdistance/IterativeCostDistance.scala
@@ -22,17 +22,15 @@ import geotrellis.raster.costdistance.SimpleCostDistance
 import geotrellis.raster.rasterize.Rasterizer
 import geotrellis.layer._
 import geotrellis.spark._
-import geotrellis.util._
 import geotrellis.vector._
+
 import org.apache.log4j.Logger
 import org.apache.spark.rdd.RDD
-import org.apache.spark.SparkContext
 import org.apache.spark.storage.StorageLevel
 import org.apache.spark.util.AccumulatorV2
 
 import scala.collection.mutable
 
-
 /**
   * This Spark-enabled implementation of the standard cost-distance
   * algorithm mentioned in the "previous work" section of [1] is
@@ -52,6 +50,27 @@ object IterativeCostDistance {
 
   val logger = Logger.getLogger(IterativeCostDistance.getClass)
 
+  /**
+    * Compute the resolution (in meters per pixel) of a layer.
+    */
+  private [spark] def computeResolution[K: (* => SpatialKey), V: (* => Tile)](
+    friction: RDD[(K, V)] with Metadata[TileLayerMetadata[K]]
+  ) = {
+    val md = friction.metadata
+    val mt = md.mapTransform
+    val key: SpatialKey = md.bounds.get.minKey
+    val extent = mt(key).reproject(md.crs, LatLng)
+
+    val latitude = (extent.ymax + extent.ymin) / 2.0
+    val degreesPerKey = extent.xmax - extent.xmin
+    // https://gis.stackexchange.com/questions/2951/algorithm-for-offsetting-a-latitude-longitude-by-some-amount-of-meters
+    val metersPerDegree = 111111 * math.cos(math.toRadians(latitude))
+    val metersPerKey = metersPerDegree * degreesPerKey
+    val keysPerPixel = 1.0 / md.layout.tileCols
+
+    metersPerKey * keysPerPixel
+  }
+
   /**
     * An accumulator to hold lists of edge changes.
     */
@@ -73,19 +92,6 @@ object IterativeCostDistance {
     def value: Changes = list
   }
 
-  def computeResolution[K: (* => SpatialKey), V: (* => Tile)](
-    friction: RDD[(K, V)] with Metadata[TileLayerMetadata[K]]
-  ) = {
-    val md = friction.metadata
-    val mt = md.mapTransform
-    val key: SpatialKey = md.bounds.get.minKey
-    val extent = mt(key).reproject(md.crs, LatLng)
-    val degrees = extent.xmax - extent.xmin
-    val meters = degrees * (6378137 * 2.0 * math.Pi) / 360.0
-    val pixels = md.layout.tileCols
-    math.abs(meters / pixels)
-  }
-
   private def geometryToKeys[K: (* => SpatialKey)](
     md: TileLayerMetadata[K],
     g: Geometry

diff --git a/spark/src/main/scala/geotrellis/spark/viewshed/IterativeViewshed.scala b/spark/src/main/scala/geotrellis/spark/viewshed/IterativeViewshed.scala
@@ -16,21 +16,17 @@
 
 package geotrellis.spark.viewshed
 
-import geotrellis.proj4.LatLng
 import geotrellis.raster._
-import geotrellis.raster.rasterize.Rasterizer
 import geotrellis.raster.viewshed.R2Viewshed
 import geotrellis.raster.viewshed.R2Viewshed._
 import geotrellis.layer._
 import geotrellis.spark._
-import geotrellis.spark.tiling._
-import geotrellis.util._
+import geotrellis.spark.costdistance.IterativeCostDistance._
 import geotrellis.vector._
 
-import org.locationtech.jts.{ geom => jts }
+import org.locationtech.jts.{geom => jts}
 import org.apache.log4j.Logger
 import org.apache.spark.rdd.RDD
-import org.apache.spark.SparkContext
 import org.apache.spark.storage.StorageLevel
 import org.apache.spark.util.AccumulatorV2
 
@@ -125,27 +121,6 @@ object IterativeViewshed {
     def value: Messages = messages.toMap
   }
 
-  /**
-    * Compute the resolution (in meters per pixel) of a layer.
-    */
-  private def computeResolution[K: (* => SpatialKey), V: (* => Tile)](
-    elevation: RDD[(K, V)] with Metadata[TileLayerMetadata[K]]
-  ) = {
-    val md = elevation.metadata
-    val mt = md.mapTransform
-    val key: SpatialKey = md.bounds.get.minKey
-    val extent = mt(key).reproject(md.crs, LatLng)
-
-    val latitude = (extent.ymax + extent.ymin) / 2.0
-    val degrees_per_key = extent.xmax - extent.xmin
-    // https://gis.stackexchange.com/questions/2951/algorithm-for-offsetting-a-latitude-longitude-by-some-amount-of-meters
-    val meters_per_degree = 111111 * math.cos(math.toRadians(latitude))
-    val meters_per_key = meters_per_degree * degrees_per_key
-    val keys_per_pixel = 1.0 / md.layout.tileCols
-
-    meters_per_key * keys_per_pixel
-  }
-
   private case class PointInfo(
     index: Int,
     key: SpatialKey,

diff --git a/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistanceSpec.scala b/spark/src/test/scala/geotrellis/spark/costdistance/IterativeCostDistanceSpec.scala
@@ -104,7 +104,7 @@ class IterativeCostDistanceSpec extends FunSpec
       val resolution = IterativeCostDistance.computeResolution(rdd3)
       val hops = (up.getDouble(2,3) - down.getDouble(2,3)) / resolution
 
-      hops should be (5.0)
+      hops should be (5.0 +- 1e-10)
     }
 
     it("Should propogate down") {
@@ -114,7 +114,7 @@ class IterativeCostDistanceSpec extends FunSpec
       val resolution = IterativeCostDistance.computeResolution(rdd3)
       val hops = (up.getDouble(2,1) - down.getDouble(2,1)) / resolution
 
-      hops should be (5.0)
+      hops should be (5.0 +- 1e-10)
     }
 
   }