Merge pull request #645 from jbouffard/improvement/hillshade

The hillshade Function Will Now Use the ZFactorCalculator

Author: jbouffard

geopyspark-backend/geotrellis/src/main/scala/geopyspark/geotrellis/SpatialTiledRasterLayer.scala

@@ -17,6 +17,7 @@ import geotrellis.raster.histogram._
 import geotrellis.raster.io.geotiff._
 import geotrellis.raster.io.geotiff.compression._
 import geotrellis.raster.mapalgebra.focal.{Square, Slope}
+import geotrellis.raster.mapalgebra.focal.hillshade._
 import geotrellis.raster.rasterize._
 import geotrellis.raster.render._
 import geotrellis.raster.resample.{ResampleMethod, PointResampleMethod, Resample}
@@ -227,6 +228,35 @@ class SpatialTiledRasterLayer(
     )
   }
 
+  def hillshade(
+    azimuth: Double,
+    altitude: Double,
+    zFactorCalculator: ZFactorCalculator,
+    band: Int
+  ): SpatialTiledRasterLayer = {
+    val mt = rdd.metadata.mapTransform
+    val cellSize = rdd.metadata.cellSize
+    val neighborhood = Square(1)
+    val gridBounds = rdd.metadata.gridBounds
+    val partitioner = rdd.partitioner
+
+    SpatialTiledRasterLayer(
+      zoomLevel,
+      rdd.withContext { rdd =>
+        rdd.bufferTiles(neighborhood.extent, gridBounds, partitioner).mapPartitions[(SpatialKey, MultibandTile)](
+        { iter =>
+          iter.map { case (key, BufferedTile(tile, bounds)) =>
+            val zfactor = zFactorCalculator.deriveZFactor(mt.keyToExtent(key))
+            val hillshadeTile = Hillshade(tile.bands(band), neighborhood, Some(bounds), cellSize, azimuth, altitude, zfactor)
+
+            key -> MultibandTile(hillshadeTile)
+          }
+        }, preservesPartitioning = true
+        )
+      }.mapContext(_.copy(cellType = ShortConstantNoDataCellType))
+    )
+  }
+
   def mask(geometries: Seq[MultiPolygon]): TiledRasterLayer[SpatialKey] =
     SpatialTiledRasterLayer(zoomLevel, Mask(rdd, geometries, Mask.Options.DEFAULT))
 
@@ -324,25 +354,6 @@ class SpatialTiledRasterLayer(
     SpatialTiledRasterLayer(None, multibandRDD)
   }
 
-  def hillshade(
-    sc: SparkContext,
-    azimuth: Double,
-    altitude: Double,
-    zFactor: Double,
-    band: Int
-  ): TiledRasterLayer[SpatialKey] = {
-    val tileLayer = TileLayerRDD(rdd.mapValues(_.band(band)), rdd.metadata)
-
-    implicit val _sc = sc
-
-    val result = tileLayer.hillshade(azimuth, altitude, zFactor)
-
-    val multibandRDD: MultibandTileLayerRDD[SpatialKey] =
-      MultibandTileLayerRDD(result.mapValues{ tile => MultibandTile(tile) }, result.metadata)
-
-    SpatialTiledRasterLayer(None, multibandRDD)
-  }
-
   def reclassify(reclassifiedRDD: RDD[(SpatialKey, MultibandTile)]): TiledRasterLayer[SpatialKey] =
     SpatialTiledRasterLayer(zoomLevel, MultibandTileLayerRDD(reclassifiedRDD, rdd.metadata))
 

geopyspark-backend/geotrellis/src/main/scala/geopyspark/geotrellis/TemporalTiledRasterLayer.scala

@@ -15,6 +15,7 @@ import geotrellis.raster.histogram._
 import geotrellis.raster.io.geotiff._
 import geotrellis.raster.io.geotiff.compression._
 import geotrellis.raster.mapalgebra.focal.{Square, Slope}
+import geotrellis.raster.mapalgebra.focal.hillshade._
 import geotrellis.raster.rasterize._
 import geotrellis.raster.render._
 import geotrellis.raster.resample.{ResampleMethod, PointResampleMethod, Resample}
@@ -351,6 +352,35 @@ class TemporalTiledRasterLayer(
     )
   }
 
+  def hillshade(
+    azimuth: Double,
+    altitude: Double,
+    zFactorCalculator: ZFactorCalculator,
+    band: Int
+  ): TemporalTiledRasterLayer = {
+    val mt = rdd.metadata.mapTransform
+    val cellSize = rdd.metadata.cellSize
+    val neighborhood = Square(1)
+    val gridBounds = rdd.metadata.gridBounds
+    val partitioner = rdd.partitioner
+
+    TemporalTiledRasterLayer(
+      zoomLevel,
+      rdd.withContext { rdd =>
+        rdd.bufferTiles(neighborhood.extent, gridBounds, partitioner).mapPartitions[(SpaceTimeKey, MultibandTile)](
+        { iter =>
+          iter.map { case (key, BufferedTile(tile, bounds)) =>
+            val zfactor = zFactorCalculator.deriveZFactor(mt.keyToExtent(key))
+            val hillshadeTile = Hillshade(tile.bands(band), neighborhood, Some(bounds), cellSize, azimuth, altitude, zfactor)
+
+            key -> MultibandTile(hillshadeTile)
+          }
+        }, preservesPartitioning = true
+        )
+      }.mapContext(_.copy(cellType = ShortConstantNoDataCellType))
+    )
+  }
+
   def costDistance(
     sc: SparkContext,
     geometries: Seq[Geometry],
@@ -375,25 +405,6 @@ class TemporalTiledRasterLayer(
     TemporalTiledRasterLayer(None, multibandRDD)
   }
 
-  def hillshade(
-    sc: SparkContext,
-    azimuth: Double,
-    altitude: Double,
-    zFactor: Double,
-    band: Int
-  ): TiledRasterLayer[SpaceTimeKey] = {
-    val tileLayer = TileLayerRDD(rdd.mapValues(_.band(band)), rdd.metadata)
-
-    implicit val _sc = sc
-
-    val result = tileLayer.hillshade(azimuth, altitude, zFactor)
-
-    val multibandRDD: MultibandTileLayerRDD[SpaceTimeKey] =
-      MultibandTileLayerRDD(result.mapValues(MultibandTile(_)), result.metadata)
-
-    TemporalTiledRasterLayer(None, multibandRDD)
-  }
-
   def reclassify(reclassifiedRDD: RDD[(SpaceTimeKey, MultibandTile)]): TiledRasterLayer[SpaceTimeKey] =
     TemporalTiledRasterLayer(zoomLevel, MultibandTileLayerRDD(reclassifiedRDD, rdd.metadata))
 

geopyspark-backend/geotrellis/src/main/scala/geopyspark/geotrellis/TiledRasterLayer.scala

@@ -198,13 +198,6 @@ abstract class TiledRasterLayer[K: SpatialComponent: JsonFormat: ClassTag: Bound
     maxDistance: Double
   ): TiledRasterLayer[K]
 
-  def hillshade(sc: SparkContext,
-    azimuth: Double,
-    altitude: Double,
-    zFactor: Double,
-    band: Int
-  ): TiledRasterLayer[K]
-
   def localMax(i: Int): TiledRasterLayer[K] =
     withRDD(rdd.mapValues { x => MultibandTile(x.bands.map(_.localMax(i))) })
 
@@ -402,6 +395,13 @@ abstract class TiledRasterLayer[K: SpatialComponent: JsonFormat: ClassTag: Bound
     }
   }
 
+  def hillshade(
+    azimuth: Double,
+    altitude: Double,
+    zFactorCalculator: ZFactorCalculator,
+    band: Int
+  ): TiledRasterLayer[K]
+
   def aggregateByCell(operation: String): TiledRasterLayer[K] = {
     val bands: RDD[(K, Array[ArrayBuffer[Tile]])] =
       rdd.combineByKey(

geopyspark/geotrellis/hillshade.py

@@ -4,13 +4,19 @@
 __all__ = ['hillshade']
 
 
-def hillshade(tiled_raster_layer, band=0, azimuth=315.0, altitude=45.0, z_factor=1.0):
+def hillshade(tiled_raster_layer, zfactor_calculator, band=0, azimuth=315.0, altitude=45.0):
     """Computes Hillshade (shaded relief) from a raster.
 
     The resulting raster will be a shaded relief map (a hill shading) based
-    on the sun altitude, azimuth, and the z factor. The z factor is a
+    on the sun altitude, azimuth, and the ``zfactor``. The ``zfactor`` is a
     conversion factor from map units to elevation units.
 
+    The ``hillshade``` operation will be carried out in a ``SQUARE`` neighborhood with with an
+    ``extent`` of 1.  The ``zfactor`` will be derived from the ``zfactor_calculator``
+    for each ``Tile`` in the Layer. The resulting Layer will have a ``cell_type``
+    of ``INT16`` regardless of the input Layer's ``cell_type``; as well as
+    have a single band, that represents the calculated ``hillshade``.
+
     Returns a raster of ShortConstantNoDataCellType.
 
     For descriptions of parameters, please see Esri Desktop's
@@ -19,17 +25,19 @@ def hillshade(tiled_raster_layer, band=0, azimuth=315.0, altitude=45.0, z_factor
     Args:
         tiled_raster_layer (:class:`~geopyspark.geotrellis.layer.TiledRasterLayer`): The base layer
             that contains the rasters used to compute the hillshade.
+        zfactor_calculator (py4j.JavaObject): A ``JavaObject`` that represents the
+            Scala ``ZFactorCalculator`` class. This can be created using either the
+            :meth:`~geopyspark.geotrellis.zfactor_lat_lng_calculator` or the
+            :meth:`~geopyspark.geotrellis.zfactor_calculator` methods.
         band (int, optional): The band of the raster to base the hillshade calculation on. Default is 0.
         azimuth (float, optional): The azimuth angle of the source of light. Default value is 315.0.
         altitude (float, optional): The angle of the altitude of the light above the horizon. Default is
             45.0.
-        z_factor (float, optional): How many x and y units in a single z unit. Default value is 1.0.
 
     Returns:
         :class:`~geopyspark.geotrellis.layer.TiledRasterLayer`
     """
 
-    srdd = tiled_raster_layer.srdd.hillshade(tiled_raster_layer.pysc._jsc.sc(), azimuth,
-                                             altitude, z_factor, band)
+    srdd = tiled_raster_layer.srdd.hillshade(azimuth, altitude, zfactor_calculator, band)
 
     return TiledRasterLayer(tiled_raster_layer.layer_type, srdd)

geopyspark/tests/geotrellis/tiled_layer_tests/hillshade_test.py

@@ -5,11 +5,11 @@
 import pytest
 
 from shapely.geometry import Point
-from geopyspark.geotrellis import SpatialKey, Tile
+from geopyspark.geotrellis import SpatialKey, Tile, zfactor_lat_lng_calculator
 from geopyspark.tests.base_test_class import BaseTestClass
 from geopyspark.geotrellis import hillshade
 from geopyspark.geotrellis.layer import TiledRasterLayer
-from geopyspark.geotrellis.constants import LayerType
+from geopyspark.geotrellis.constants import LayerType, Unit
 
 
 class HillshadeTest(BaseTestClass):
@@ -45,7 +45,8 @@ def tearDown(self):
         BaseTestClass.pysc._gateway.close()
 
     def test_hillshade(self):
-        result = hillshade(self.raster_rdd, band=0, azimuth=99.0, altitude=33.0, z_factor=0.0)
+        calc = zfactor_lat_lng_calculator(Unit.METERS)
+        result = hillshade(self.raster_rdd, calc, band=0, azimuth=99.0, altitude=33.0)
 
         data = result.to_numpy_rdd().first()[1].cells[0][0][0]
         self.assertEqual(data, 63)