Merge pull request #70 from locationtech/fix/68

Partial fix for 68

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/main/scala/astraea/spark/rasterframes/util/debug/package.scala

@@ -21,6 +21,10 @@
 package astraea.spark.rasterframes.util
 
 import astraea.spark.rasterframes._
+import geotrellis.proj4.LatLng
+import geotrellis.vector.{Feature, Geometry}
+import geotrellis.vector.io.json.JsonFeatureCollection
+import spray.json.JsValue
 
 /**
  * Additional debugging routines. No guarantees these are or will remain stable.
@@ -34,5 +38,19 @@ package object debug {
     def describeFullSchema: String = {
       self.schema.prettyJson
     }
+
+    /** Renders all the extents in this RasterFrame as GeoJSON in EPSG:4326. This does a full
+     * table scan and collects **all** the geometry into the driver, and then converts it into a
+     * Spray JSON data structure. Not performant, and for debugging only. */
+    def geoJsonExtents: JsValue = {
+      import spray.json.DefaultJsonProtocol._
+
+      val features = self
+        .select(BOUNDS_COLUMN, SPATIAL_KEY_COLUMN)
+        .collect()
+        .map{ case (p, s) ⇒ Feature(Geometry(p).reproject(self.crs, LatLng), Map("col" -> s.col, "row" -> s.row)) }
+
+      JsonFeatureCollection(features).toJson
+    }
   }
 }

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,13 +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.{DataFrame, SQLContext, SaveMode}
+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.
@@ -50,25 +53,62 @@ class DefaultSource extends DataSourceRegister with RelationProvider with Creata
     require(pathO.get.getScheme == "file" || pathO.get.getScheme == null, "Currently only 'file://' destinations are supported")
     sqlContext.withRasterFrames
 
-
     require(data.isRF, "GeoTIFF can only be constructed from a RasterFrame")
     val rf = data.certify
 
-    val tl = rf.tileLayerMetadata.merge.layout.tileLayout
+    // If no desired image size is given, write at full size.
+    lazy val (fullResCols, fullResRows) = {
+      // get the layout size given that the tiles may be heterogenously sized
+      // first get any valid row and column in the spatial key structure
+      val sk = rf.select(SPATIAL_KEY_COLUMN).first()
+
+      val tc = rf.tileColumns.head
+
+      val c = rf
+        .where(SPATIAL_KEY_COLUMN("row") === sk.row)
+        .agg(
+          F.sum(tileDimensions(tc)("cols") cast(LongType))
+        ).first()
+        .getLong(0)
+
+      val r = rf
+        .where(SPATIAL_KEY_COLUMN("col") === sk.col)
+        .agg(
+          F.sum(tileDimensions(tc)("rows") cast(LongType))
+        ).first()
+        .getLong(0)
 
-    val cols = numParam(DefaultSource.IMAGE_WIDTH_PARAM, parameters).getOrElse(tl.totalCols)
-    val rows = numParam(DefaultSource.IMAGE_HEIGHT_PARAM, parameters).getOrElse(tl.totalRows)
+      (c, r)
+    }
+
+    val cols = numParam(DefaultSource.IMAGE_WIDTH_PARAM, parameters).getOrElse(fullResCols)
+    val rows = numParam(DefaultSource.IMAGE_HEIGHT_PARAM, parameters).getOrElse(fullResRows)
 
     require(cols <= Int.MaxValue && rows <= Int.MaxValue, s"Can't construct a GeoTIFF of size $cols x $rows. (Too big!)")
 
     // Should we really play traffic cop here?
     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.")
 
-    println()
-    val raster = rf.toMultibandRaster(rf.tileColumns, cols.toInt, rows.toInt)
+    val tcols = rf.tileColumns
+    val raster = rf.toMultibandRaster(tcols, cols.toInt, rows.toInt)
+
+    // 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)
 
-    GeoTiff(raster).write(pathO.get.getPath)
+    logger.debug(s"Writing DataFrame to GeoTIFF ($cols x $rows) at ${pathO.get}")
+    geotiff.write(pathO.get.getPath)
     GeoTiffRelation(sqlContext, pathO.get)
   }
 }
@@ -78,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

@@ -27,9 +27,10 @@ import geotrellis.raster.io.geotiff.reader.GeoTiffReader.GeoTiffInfo
 import geotrellis.spark.{KeyBounds, SpatialKey, TileLayerMetadata}
 import geotrellis.spark.tiling.LayoutDefinition
 import geotrellis.util.ByteReader
+import geotrellis.vector.Extent
 
 /**
- * Utility mix in for generating a tlm from GeoTiff headers.
+ * Utility mix-in for generating a tlm from GeoTiff headers.
  *
  * @since 5/4/18
  */
@@ -49,23 +50,35 @@ trait GeoTiffInfoSupport {
 
   def extractGeoTiffLayout(reader: ByteReader): (GeoTiffReader.GeoTiffInfo, TileLayerMetadata[SpatialKey]) = {
     val info: GeoTiffInfo = Shims.readGeoTiffInfo(reader, false, true)
+    // Some notes on GeoTiffInfo properties:
+    // * `info.extent` is the actual geotiff extent
+    // * `info.segmentLayout.tileLayout` contains the internal, regularized gridding of a tiled GeoTIFF
+    // * `info.segmentLayout.tileLayout.{totalCols|totalRows}` is the largest number of possible cells in the internal gridding
+    // * `info.segmentLayout.{totalCols|totalRows}` are the real dimensions of the GeoTIFF. This is likely smaller than
+    //     the total size of `info.segmentLayout.tileLayout.{totalCols|totalRows}`
+    // * `info.rasterExtent.{cellwidth|cellheight}` is the per-pixel spatial resolution
+    // * `info.extent` and `info.rasterExtent.extent` are the same thing
+
     val tlm = {
-      val layout = if(!info.segmentLayout.isTiled) {
+      val tileLayout = if(info.segmentLayout.isTiled) {
+        info.segmentLayout.tileLayout
+      }
+      else {
         val width = info.segmentLayout.totalCols
         val height = info.segmentLayout.totalRows
         defaultLayout(width, height)
       }
-      else {
-        info.segmentLayout.tileLayout
-      }
       val extent = info.extent
       val crs = info.crs
       val cellType = info.cellType
       val bounds = KeyBounds(
         SpatialKey(0, 0),
-        SpatialKey(layout.layoutCols - 1, layout.layoutRows - 1)
+        SpatialKey(tileLayout.layoutCols - 1, tileLayout.layoutRows - 1)
       )
-      TileLayerMetadata(cellType, LayoutDefinition(extent, layout), 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,13 +106,14 @@ class GeoTiffDataSourceSpec
     }
 
     it("should write GeoTIFF") {
-
       val rf = spark.read
         .geotiff
         .loadRF(cogPath)
 
-      val out = Paths.get(outputLocalPath, "example-geotiff.tiff")
-      //val out = Paths.get("target", "example-geotiff.tiff")
+      logger.info(s"Read extent: ${rf.tileLayerMetadata.merge.extent}")
+
+      val out = Paths.get("target", "example-geotiff.tiff")
+      logger.info(s"Writing to $out")
       noException shouldBe thrownBy {
         rf.write.geotiff.save(out.toString)
       }

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">