Fixes #68

Added better natural color rendering to slippy map export.

Signed-off-by: Simeon H.K. Fitch <[email protected]>

Author: metasim

core/src/main/scala/astraea/spark/rasterframes/util/MultibandRender.scala

@@ -21,7 +21,9 @@
 package astraea.spark.rasterframes.util
 
 import geotrellis.raster._
-import geotrellis.raster.render.Png
+import geotrellis.raster.render.{ColorRamp, Png}
+
+import scala.util.Try
 
 /**
  * Rework of process courtesy of @lossyrob for creating natural color RGB images in GeoTrellis.
@@ -55,6 +57,11 @@ object MultibandRender {
   }
   import CellTransforms._
 
+  /** Create an RGB composite PNG file from the given MultibandTile and color profile. */
+  @deprecated("Use Profile.render instead", "0.7.0")
+  def rgbComposite(tile: MultibandTile, profile: Profile): Png = profile.render(tile)
+
+  /** Base type for Rendering profiles. */
   trait Profile {
     /** Value from -255 to 255 */
     def brightness: Int = 0
@@ -73,8 +80,10 @@ object MultibandRender {
     /** Convert the tile to an Int-based cell type. */
     def normalizeCellType(tile: Tile): Tile = tile.convert(IntCellType)
 
-    /** Convert tile such that cells values fall between 0 and 255. */
-    def compressRange(tile: Tile): Tile = tile
+    /** Convert tile such that cells values fall between 0 and 255, if desired. */
+    def compressRange(tile: Tile): Tile =
+      // `Try` below is due to https://github.com/locationtech/geotrellis/issues/2621
+      Try(tile.rescale(0, 255)).getOrElse(tile)
 
     /** Apply color correction so it "looks nice". */
     def colorAdjust(tile: Tile): Tile = {
@@ -89,7 +98,15 @@ object MultibandRender {
       normalizeCellType(tile).map(pipeline)
     }
 
-    val applyAdjustment = compressRange _ andThen colorAdjust
+    val applyAdjustment: Tile ⇒ Tile =
+      compressRange _ andThen colorAdjust
+
+    def render(tile: MultibandTile) = {
+      val r = applyAdjustment(red(tile))
+      val g = applyAdjustment(green(tile))
+      val b = applyAdjustment(blue(tile))
+      ArrayMultibandTile(r, g, b).renderPng
+    }
   }
   case object Default extends Profile
 
@@ -108,13 +125,11 @@ object MultibandRender {
 
   case object NAIPNaturalColor extends Profile {
     override val gamma = 1.4
-    override def compressRange(tile: Tile): Tile = tile.rescale(0, 255)
   }
 
-  def rgbComposite(tile: MultibandTile, profile: Profile): Png = {
-    val red = profile.applyAdjustment(profile.red(tile))
-    val green = profile.applyAdjustment(profile.green(tile))
-    val blue = profile.applyAdjustment(profile.blue(tile))
-    ArrayMultibandTile(red, green, blue).renderPng
+  case class ColorRampProfile(ramp: ColorRamp) extends Profile {
+    // Are there other ways to use the other bands?
+    override def render(tile: MultibandTile): Png =
+      colorAdjust(tile.band(0)).renderPng(ramp)
   }
 }

core/src/test/scala/astraea/spark/rasterframes/RasterFrameSpec.scala

@@ -301,7 +301,7 @@ class RasterFrameSpec extends TestEnvironment with MetadataKeys
 
         noException shouldBe thrownBy {
           val raster = joined.toMultibandRaster(Seq($"red", $"green", $"blue"), 128, 128)
-          val png = MultibandRender.rgbComposite(raster.tile, MultibandRender.Landsat8NaturalColor)
+          val png = MultibandRender.Landsat8NaturalColor.render(raster.tile)
           //png.write(s"target/${getClass.getSimpleName}.png")
         }
       }
@@ -313,7 +313,7 @@ class RasterFrameSpec extends TestEnvironment with MetadataKeys
 
         noException shouldBe thrownBy {
           val raster = joined.toMultibandRaster(Seq($"red", $"green", $"blue"), 256, 256)
-          val png = MultibandRender.rgbComposite(raster.tile, MultibandRender.NAIPNaturalColor)
+          val png = MultibandRender.NAIPNaturalColor.render(raster.tile)
           png.write(s"target/${getClass.getSimpleName}.png")
         }
       }

core/src/test/scala/astraea/spark/rasterframes/TestData.scala

@@ -25,7 +25,7 @@ import com.vividsolutions.jts.geom.{Coordinate, GeometryFactory}
 import geotrellis.proj4.LatLng
 import geotrellis.raster
 import geotrellis.raster._
-import geotrellis.raster.io.geotiff.SinglebandGeoTiff
+import geotrellis.raster.io.geotiff.{MultibandGeoTiff, SinglebandGeoTiff}
 import geotrellis.spark._
 import geotrellis.spark.testkit.TileLayerRDDBuilders
 import geotrellis.spark.tiling.LayoutDefinition
@@ -89,6 +89,7 @@ trait TestData {
   }
 
   def readSingleband(name: String) = SinglebandGeoTiff(IOUtils.toByteArray(getClass.getResourceAsStream("/" + name)))
+  def readMultiband(name: String) = MultibandGeoTiff(IOUtils.toByteArray(getClass.getResourceAsStream("/" + name)))
 
   /** 774 x 500 GeoTiff */
   def sampleGeoTiff = readSingleband("L8-B8-Robinson-IL.tiff")
@@ -106,6 +107,8 @@ trait TestData {
     readSingleband(s"NAIP-VA-b$band.tiff")
   }
 
+  def rgbCogSample = readMultiband("LC08_RGB_Norfolk_COG.tiff")
+
   def sampleTileLayerRDD(implicit spark: SparkSession): TileLayerRDD[SpatialKey] = {
     val rf = sampleGeoTiff.projectedRaster.toRF(128, 128)
     rf.toTileLayerRDD(rf.tileColumns.head).left.get

datasource/src/main/scala/astraea/spark/rasterframes/datasource/geotiff/DefaultSource.scala

@@ -19,14 +19,16 @@
 
 package astraea.spark.rasterframes.datasource.geotiff
 
-import _root_.geotrellis.raster.io.geotiff.GeoTiff
 import astraea.spark.rasterframes._
 import astraea.spark.rasterframes.datasource._
 import astraea.spark.rasterframes.util._
 import com.typesafe.scalalogging.LazyLogging
 import org.apache.spark.sql.sources.{BaseRelation, CreatableRelationProvider, DataSourceRegister, RelationProvider}
 import org.apache.spark.sql.types.LongType
 import org.apache.spark.sql.{DataFrame, SQLContext, SaveMode, functions ⇒ F}
+import _root_.geotrellis.raster.io.geotiff.{GeoTiffOptions, MultibandGeoTiff, Tags, Tiled}
+import _root_.geotrellis.raster.io.geotiff.compression._
+import _root_.geotrellis.raster.io.geotiff.tags.codes.ColorSpace
 
 /**
  * Spark SQL data source over GeoTIFF files.
@@ -88,10 +90,25 @@ class DefaultSource extends DataSourceRegister with RelationProvider with Creata
     if(cols.toDouble * rows * 64.0 > Runtime.getRuntime.totalMemory() * 0.5)
       logger.warn(s"You've asked for the construction of a very large image ($cols x $rows), destined for ${pathO.get}. Out of memory error likely.")
 
-    logger.debug(s"Writing DataFrame to GeoTIFF ($cols by $rows) at ${pathO.get}")
-    val raster = rf.toMultibandRaster(rf.tileColumns, cols.toInt, rows.toInt)
+    val tcols = rf.tileColumns
+    val raster = rf.toMultibandRaster(tcols, cols.toInt, rows.toInt)
 
-    GeoTiff(raster).write(pathO.get.getPath)
+    // We make some assumptions here.... eventually have column metadata encode this.
+    val colorSpace = tcols.size match {
+      case 3 | 4 ⇒  ColorSpace.RGB
+      case _ ⇒ ColorSpace.BlackIsZero
+    }
+
+    val compress = parameters.get(DefaultSource.COMPRESS).map(_.toBoolean).getOrElse(false)
+    val options = GeoTiffOptions(Tiled, if (compress) DeflateCompression else NoCompression, colorSpace)
+    val tags = Tags(
+      RFBuildInfo.toMap.filter(_._1.startsWith("rf")).mapValues(_.toString),
+      tcols.map(c ⇒ Map("RF_COL" -> c.columnName)).toList
+    )
+    val geotiff = new MultibandGeoTiff(raster.tile, raster.extent, raster.crs, tags, options)
+
+    logger.debug(s"Writing DataFrame to GeoTIFF ($cols x $rows) at ${pathO.get}")
+    geotiff.write(pathO.get.getPath)
     GeoTiffRelation(sqlContext, pathO.get)
   }
 }
@@ -101,4 +118,5 @@ object DefaultSource {
   final val PATH_PARAM = "path"
   final val IMAGE_WIDTH_PARAM = "imageWidth"
   final val IMAGE_HEIGHT_PARAM = "imageWidth"
+  final val COMPRESS = "compress"
 }

datasource/src/main/scala/astraea/spark/rasterframes/datasource/geotiff/GeoTiffInfoSupport.scala

@@ -59,7 +59,6 @@ trait GeoTiffInfoSupport {
     // * `info.rasterExtent.{cellwidth|cellheight}` is the per-pixel spatial resolution
     // * `info.extent` and `info.rasterExtent.extent` are the same thing
 
-
     val tlm = {
       val tileLayout = if(info.segmentLayout.isTiled) {
         info.segmentLayout.tileLayout
@@ -77,17 +76,9 @@ trait GeoTiffInfoSupport {
         SpatialKey(tileLayout.layoutCols - 1, tileLayout.layoutRows - 1)
       )
 
-      val layoutExtentWidth = tileLayout.totalCols * info.rasterExtent.cellwidth
-      val layoutExtentHeight = tileLayout.totalRows * info.rasterExtent.cellheight
-
-      val layoutExtent = Extent(
-        extent.xmin,
-        extent.ymin,
-        extent.xmin + layoutExtentWidth,
-        extent.ymin + layoutExtentHeight
-      )
-
-      TileLayerMetadata(cellType, LayoutDefinition(layoutExtent, tileLayout), extent, crs, bounds)
+      TileLayerMetadata(cellType,
+        LayoutDefinition(info.rasterExtent, tileLayout.tileCols, tileLayout.tileRows),
+        extent, crs, bounds)
     }
 
     (info, tlm)

datasource/src/main/scala/astraea/spark/rasterframes/datasource/geotiff/GeoTiffRelation.scala

@@ -23,12 +23,9 @@ import java.net.URI
 
 import astraea.spark.rasterframes._
 import astraea.spark.rasterframes.util._
-import geotrellis.raster.{MultibandTile, Tile, TileLayout}
-import geotrellis.raster.io.geotiff.reader.GeoTiffReader
 import geotrellis.spark._
 import geotrellis.spark.io._
 import geotrellis.spark.io.hadoop._
-import geotrellis.spark.tiling.LayoutDefinition
 import geotrellis.util._
 import org.apache.hadoop.fs.Path
 import org.apache.spark.rdd.RDD

datasource/src/test/scala/astraea/spark/rasterframes/datasource/geotiff/GeoTiffDataSourceSpec.scala

@@ -21,7 +21,10 @@ package astraea.spark.rasterframes.datasource.geotiff
 import java.nio.file.Paths
 
 import astraea.spark.rasterframes._
+import geotrellis.proj4.LatLng
 import org.apache.spark.sql.functions._
+import geotrellis.vector._
+import geotrellis.vector.io.json.JsonFeatureCollection
 
 /**
  * @since 1/14/18
@@ -103,33 +106,17 @@ class GeoTiffDataSourceSpec
     }
 
     it("should write GeoTIFF") {
-
       val rf = spark.read
         .geotiff
         .loadRF(cogPath)
 
       logger.info(s"Read extent: ${rf.tileLayerMetadata.merge.extent}")
 
-      val cellsize = rf.groupBy(
-        tileDimensions(col("tile_1"))("rows") as "tileRows",
-        tileDimensions(col("tile_1"))("cols") as "tileCols"
-      ).count()
-
-      cellsize.show(false)
-
-      rf.withBounds().select(
-      tileDimensions(col("tile_1"))("rows") as "tileRows",
-      tileDimensions(col("tile_1"))("cols") as "tileCols",
-        BOUNDS_COLUMN)
-        .show(false)
-
-      val out = Paths.get(outputLocalPath, "example-geotiff.tiff")
+      val out = Paths.get("target", "example-geotiff.tiff")
       logger.info(s"Writing to $out")
-      //val out = Paths.get("target", "example-geotiff.tiff")
       noException shouldBe thrownBy {
         rf.write.geotiff.save(out.toString)
       }
-      logger.info("Done!")
     }
   }
 }

experimental/src/it/scala/astraea/spark/rasterframes/experimental/slippy/SlippyExportDriver.scala

@@ -27,6 +27,7 @@ import geotrellis.raster._
 import geotrellis.raster.io.geotiff.SinglebandGeoTiff
 import org.apache.spark.sql.SparkSession
 import SlippyExport._
+import astraea.spark.rasterframes.util.MultibandRender
 
 object SlippyExportDriver {
   def main(args: Array[String]): Unit = {
@@ -38,21 +39,34 @@ object SlippyExportDriver {
       .getOrCreate()
       .withRasterFrames
 
-    val bands: Seq[SinglebandGeoTiff] = for (i ← 1 to 3) yield {
-      TestData.readSingleband(s"NAIP-VA-b$i.tiff")
-      //s"L8-B$i-Elkton-VA.tiff")
-    }
+    def mergeBands = {
+      val bands: Seq[SinglebandGeoTiff] = for (i ← 1 to 3) yield {
+        TestData.l8Sample(i)
+      }
+
+      val mtile = MultibandTile(bands.map(_.tile))
 
-    val mtile = MultibandTile(bands.map(_.tile))
+      val pr = ProjectedRaster(mtile, bands.head.extent, bands.head.crs)
 
-    val pr = ProjectedRaster(mtile, bands.head.extent, bands.head.crs)
+      implicit val bandCount = PairRDDConverter.forSpatialMultiband(bands.length)
 
-    implicit val bandCount = PairRDDConverter.forSpatialMultiband(bands.length)
+      val rf = pr.toRF(64, 64)
 
-    val rf = pr.toRF(64, 64)
+      //rf.exportGeoTiffTiles(new File("target/slippy-tiff").toURI)
 
-    //rf.exportGeoTiffTiles(new File("target/slippy-tiff").toURI)
+      rf.exportSlippyMap(new File("target/slippy-1/").toURI)
+    }
+
+   def multiband = {
+     implicit val bandCount = PairRDDConverter.forSpatialMultiband(3)
+     val rf = TestData.rgbCogSample.projectedRaster.toRF(256, 256)
+     import astraea.spark.rasterframes.datasource.geotiff._
+     //val rf = spark.read.geotiff.loadRF(getClass.getResource("/LC08_RGB_Norfolk_COG.tiff").toURI)
+     println(rf.tileLayerMetadata.merge.toString)
+     rf.exportSlippyMap(new File("target/slippy-2/").toURI, MultibandRender.Landsat8NaturalColor)
+    }
 
-    rf.exportSlippyMap(new File("target/slippy-png/").toURI)
+    multiband
+    //mergeBands
   }
 }

experimental/src/main/resources/slippy.html

@@ -21,7 +21,6 @@
 
 <html lang="en">
 <head>
-    <meta charset="UTF-8">
     <title>RasterFrames</title>
     <meta charset="utf-8" />
    	<meta name="viewport" content="width=device-width, initial-scale=1.0">

experimental/src/main/scala/astraea/spark/rasterframes/experimental/slippy/SlippyExport.scala

@@ -20,15 +20,16 @@
 
 package astraea.spark.rasterframes.experimental.slippy
 
+import java.io.PrintStream
 import java.net.URI
 
 import astraea.spark.rasterframes._
+import astraea.spark.rasterframes.util.MultibandRender.ColorRampProfile
 import astraea.spark.rasterframes.util._
 import geotrellis.proj4.{LatLng, WebMercator}
 import geotrellis.raster._
 import geotrellis.raster.io.geotiff._
 import geotrellis.raster.io.geotiff.tags.codes.ColorSpace
-import geotrellis.raster.render.{ColorMap, ColorRamps}
 import geotrellis.raster.resample.Bilinear
 import geotrellis.spark._
 import geotrellis.spark.io.slippy.HadoopSlippyTileWriter
@@ -41,7 +42,6 @@ import org.apache.hadoop.fs.{FileSystem, Path}
 import org.apache.spark.annotation.Experimental
 
 import scala.io.Source
-import scala.util.Try
 
 /**
  * Experimental support for exporting a RasterFrame into Slippy map format.
@@ -94,7 +94,7 @@ trait SlippyExport extends MethodExtensions[RasterFrame]{
    * @param colorMap Optional color map to use for rendering tiles in non-RGB RasterFrames.
    */
   @Experimental
-  def exportSlippyMap(dest: URI, colorMap: Option[ColorMap] = None): Unit = {
+  def exportSlippyMap(dest: URI, renderer: MultibandRender.Profile = MultibandRender.Default): Unit = {
     val spark = self.sparkSession
     implicit val sc = spark.sparkContext
 
@@ -116,15 +116,9 @@ trait SlippyExport extends MethodExtensions[RasterFrame]{
 
     val (zoom, reprojected) = inputRDD.reproject(WebMercator, layoutScheme, Bilinear)
     val writer = new HadoopSlippyTileWriter[MultibandTile](dest.toASCIIString + "/" + tileDirName, "png")({ (_, tile) =>
-      val png = if(colorMap.isEmpty && tile.bands.lengthCompare(3) == 0) {
-        // `Try` below is due to https://github.com/locationtech/geotrellis/issues/2621
-        tile.mapBands((_, t) ⇒ Try(t.rescale(0, 255)).getOrElse(t)).renderPng()
-      }
-      else {
-        // Are there other ways to use the other bands?
-        val selected = tile.bands.head
-        colorMap.map(m ⇒ selected.renderPng(m)).getOrElse(selected.renderPng(ColorRamps.greyscale(256)))
-      }
+      require(tile.bandCount >= 3 || renderer.isInstanceOf[ColorRampProfile],
+        "Single-band and dual-band RasterFrames require a ColorRampProfile for rendering")
+      val png = renderer.render(tile)
       png.bytes
     })
 
@@ -153,10 +147,10 @@ object SlippyExport {
     val subst = new StrSubstitutor(subs.asJava)
 
     val fs = FileSystem.get(dest, conf)
-    withResource(fs.create(new Path(new Path(dest), "index.html"), true)) { out ⇒
+    withResource(fs.create(new Path(new Path(dest), "index.html"), true)) { hout ⇒
+      val out = new PrintStream(hout, true ,"UTF-8")
       for(line ← rawLines) {
-        out.writeBytes(subst.replace(line))
-        out.writeChar('\n')
+        out.println(subst.replace(line))
       }
     }
   }