Ability to read/write color tables for GeoTIFFs encoded with palette photometric interpretation

raster-test/src/test/scala/geotrellis/raster/io/geotiff/reader/GeoTiffReaderSpec.scala

@@ -16,28 +16,18 @@
 
 package geotrellis.raster.io.geotiff.reader
 
+import geotrellis.proj4._
 import geotrellis.raster._
 import geotrellis.raster.io.arg._
 import geotrellis.raster.io.geotiff._
-import geotrellis.raster.io.geotiff.util._
 import geotrellis.raster.io.geotiff.tags._
-import geotrellis.raster.summary.polygonal._
-
-import geotrellis.vector.{Point, Extent}
+import geotrellis.raster.io.geotiff.tags.codes.ColorSpace
+import geotrellis.raster.render.RGB
 import geotrellis.raster.testkit._
-import geotrellis.proj4._
-
+import geotrellis.vector.{Extent, Point}
 import monocle.syntax.apply._
-import org.scalactic.Tolerance
-
-import scala.io.{Source, Codec}
-import scala.collection.immutable.HashMap
-
-import java.util.BitSet
-import java.nio.ByteBuffer
-
-import spire.syntax.cfor._
 import org.scalatest._
+import spire.syntax.cfor._
 
 class GeoTiffReaderSpec extends FunSpec
     with Matchers
@@ -81,8 +71,6 @@ class GeoTiffReaderSpec extends FunSpec
     val tile = compressed.tile
     val bounds = tile.gridBounds
     bounds.width should be (1121)
-
-    import Tolerance._
     compressed.crs should be (CRS.fromName("EPSG:4326"))
     if(compressed.extent.min.distance(Point(59.9955397,  30.0044603))>0.0001) {
       compressed.extent.min should be (Point(59.9955397,  30.0044603))
@@ -473,6 +461,36 @@ class GeoTiffReaderSpec extends FunSpec
       val geoTiff = SinglebandGeoTiff.compressed(s"$baseDataPath/sbn/SBN_inc_percap-nodata-clip.tif")
       geoTiff.tile.cellType should be (ByteConstantNoDataCellType)
     }
+
+    it("should read photometric interpretation code") {
+      val expected = Map(
+        "colormap.tif" -> ColorSpace.Palette,
+        "multi-tag.tif" -> ColorSpace.RGB,
+        "alaska-polar-3572.tif" -> ColorSpace.BlackIsZero,
+        "3bands/bit/3bands-striped-band.tif" -> ColorSpace.RGB
+      )
+
+      Inspectors.forEvery(expected) {
+        case (file, space) ⇒
+          MultibandGeoTiff(geoTiffPath(file)).options.colorSpace should be (space)
+      }
+    }
+
+
+    it("should read and convert color table") {
+      val geoTiff = SinglebandGeoTiff.compressed(geoTiffPath("colormap.tif"))
+
+      geoTiff.options.colorMap should be ('defined)
+
+      val cmap = geoTiff.options.colorMap.get
+      cmap.colors.size should be (256)
+      // These was determined by inspecting the color table
+      cmap.map(0) should be (RGB(0, 0, 0))
+      cmap.map(1) should be (RGB(0, 249, 0))
+      cmap.map(12) should be (RGB(209, 221, 249))
+      cmap.map(95) should be (RGB(112, 163, 186))
+      cmap.map(255) should be (RGB(0, 0, 0))
+    }
   }
 
   describe("Reading and writing special metadata tags ") {

raster-test/src/test/scala/geotrellis/raster/io/geotiff/writer/GeoTiffWriterSpec.scala

@@ -16,16 +16,15 @@
 
 package geotrellis.raster.io.geotiff.writer
 
-import geotrellis.raster._
-import geotrellis.raster.io.geotiff._
-import geotrellis.raster.io.geotiff.reader._
-import geotrellis.vector.Extent
-import geotrellis.proj4.CRS
-import geotrellis.proj4.LatLng
-import geotrellis.raster.testkit._
 import java.io._
 
+import geotrellis.proj4.{CRS, LatLng}
+import geotrellis.raster._
+import geotrellis.raster.io.geotiff._
 import geotrellis.raster.io.geotiff.tags.codes.ColorSpace
+import geotrellis.raster.render.{ColorRamps, IndexedColorMap}
+import geotrellis.raster.testkit._
+import geotrellis.vector.Extent
 import org.scalatest._
 
 class GeoTiffWriterSpec extends FunSpec
@@ -38,10 +37,11 @@ class GeoTiffWriterSpec extends FunSpec
   override def afterAll = purge
 
   private val testCRS = CRS.fromName("EPSG:3857")
+  private val testExtent = Extent(100.0, 400.0, 120.0, 420.0)
 
   describe ("writing GeoTiffs without errors and with correct tiles, crs and extent") {
     val temp = File.createTempFile("geotiff-writer", ".tif")
-    val path = temp.getPath()
+    val path = temp.getPath
 
     it("should write GeoTiff with tags") {
       val geoTiff = MultibandGeoTiff(geoTiffPath("multi-tag.tif"))
@@ -105,19 +105,18 @@ class GeoTiffWriterSpec extends FunSpec
       actualCRS.epsgCode should be (geoTiff.crs.epsgCode)
     }
 
-    it ("should write floating point rasters correctly") {
-      val e = Extent(100.0, 400.0, 120.0, 420.0)
+    it("should write floating point rasters correctly") {
       val t = DoubleArrayTile(Array(11.0, 22.0, 33.0, 44.0), 2, 2)
 
-      val geoTiff = SinglebandGeoTiff(t, e, testCRS, Tags.empty, GeoTiffOptions.DEFAULT)
+      val geoTiff = SinglebandGeoTiff(t, testExtent, testCRS, Tags.empty, GeoTiffOptions.DEFAULT)
 
       GeoTiffWriter.write(geoTiff, path)
 
       addToPurge(path)
 
       val SinglebandGeoTiff(tile, extent, crs, _, _) = SinglebandGeoTiff(path)
 
-      extent should equal (e)
+      extent should equal (testExtent)
       crs should equal (testCRS)
       assertEqual(tile, t)
     }
@@ -210,20 +209,12 @@ class GeoTiffWriterSpec extends FunSpec
         assertEqual(actualBand, expectedBand)
       }
     }
+  }
+  describe ("writing GeoTiffs with correct color handling") {
+    val temp = File.createTempFile("geotiff-writer", ".tif")
+    val path = temp.getPath
+
 
-    it("should read photometric interpretation code") {
-      val expected = Map(
-        "colormap.tif" -> ColorSpace.Palette,
-        "multi-tag.tif" -> ColorSpace.RGB,
-        "alaska-polar-3572.tif" -> ColorSpace.BlackIsZero,
-        "3bands/bit/3bands-striped-band.tif" -> ColorSpace.RGB
-      )
-
-      Inspectors.forEvery(expected) {
-        case (file, space) ⇒
-          MultibandGeoTiff(geoTiffPath(file)).options.colorSpace should be (space)
-      }
-    }
 
     it("should write photometric interpretation code") {
       // Read in a 4-band file interpreted as RGB(A)
@@ -239,7 +230,64 @@ class GeoTiffWriterSpec extends FunSpec
 
       val reread = MultibandGeoTiff(path)
 
+      addToPurge(path)
+
       reread.options.colorSpace should be (ColorSpace.CMYK)
     }
+
+    it("should write color map when photometric interpretation is 'Palette'") {
+      val hundreds = createConsecutiveTile(10).map(_ - 1).convert(ByteCellType)
+
+      val colorMap = ColorRamps.HeatmapBlueToYellowToRedSpectrum
+        .stops(100)
+        .toColorMap(Array.tabulate[Int](100)(identity))
+      val indexedColorMap = IndexedColorMap.fromColorMap(colorMap)
+
+      val indexed = SinglebandGeoTiff(hundreds, testExtent, testCRS, Tags.empty, GeoTiffOptions(indexedColorMap))
+
+      GeoTiffWriter.write(indexed, path)
+
+      val reread = MultibandGeoTiff(path)
+
+      addToPurge(path)
+
+      reread.options.colorSpace should be (ColorSpace.Palette)
+      reread.options.colorMap should be('defined)
+
+      val p1 = reread.options.colorMap.get.colors
+      val p2 = indexed.options.colorMap.get.colors
+
+      Inspectors.forEvery(p1.zip(p2)) { case (c1, c2) ⇒
+        c1 should equal (c2)
+      }
+    }
+
+    it("should preserve color map in existing file") {
+      val base = SinglebandGeoTiff(geoTiffPath("colormap.tif"))
+      GeoTiffWriter.write(base, path)
+
+      val reread = MultibandGeoTiff(path)
+
+      addToPurge(path)
+
+      val p1 = reread.options.colorMap.get.colors
+      val p2 = base.options.colorMap.get.colors
+
+      Inspectors.forEvery(p1.zip(p2)) { case (c1, c2) ⇒
+        c1 should equal (c2)
+      }
+    }
+
+    it("should inhibit writing unsupported 'Palette' color map configuration") {
+      val base = SinglebandGeoTiff(geoTiffPath("colormap.tif"))
+      val illegal = Seq(FloatCellType, DoubleCellType, IntCellType)
+
+      Inspectors.forEvery(illegal) { cellType ⇒
+        val naughty = base.copy(tile = base.tile.convert(cellType))
+        intercept[IncompatibleGeoTiffOptionsException] {
+          GeoTiffWriter.write(naughty, path)
+        }
+      }
+    }
   }
 }

raster/src/main/scala/geotrellis/raster/io/geotiff/GeoTiffOptions.scala

@@ -2,6 +2,7 @@ package geotrellis.raster.io.geotiff
 
 import geotrellis.raster.io.geotiff.compression._
 import geotrellis.raster.io.geotiff.tags.codes.ColorSpace
+import geotrellis.raster.render.IndexedColorMap
 
 /**
   * This case class holds information about how the data is stored in
@@ -11,26 +12,33 @@ import geotrellis.raster.io.geotiff.tags.codes.ColorSpace
 case class GeoTiffOptions(
   storageMethod: StorageMethod = GeoTiffOptions.DEFAULT.storageMethod,
   compression: Compression = GeoTiffOptions.DEFAULT.compression,
-  colorSpace: Int = GeoTiffOptions.DEFAULT.colorSpace
+  colorSpace: Int = GeoTiffOptions.DEFAULT.colorSpace,
+  colorMap: Option[IndexedColorMap] = GeoTiffOptions.DEFAULT.colorMap
 )
 
 /**
-  * The companion object to [[GeoTiffOptions]]
-  */
+ * The companion object to [[GeoTiffOptions]]
+ */
 object GeoTiffOptions {
-  val DEFAULT = GeoTiffOptions(Striped, NoCompression, ColorSpace.BlackIsZero)
+  val DEFAULT = GeoTiffOptions(Striped, NoCompression, ColorSpace.BlackIsZero, None)
 
   /**
-    * Creates a new instance of [[GeoTiffOptions]] with the given
-    * StorageMethod and the default compression value
-    */
+   * Creates a new instance of [[GeoTiffOptions]] with the given
+   * StorageMethod and the default compression value
+   */
   def apply(storageMethod: StorageMethod): GeoTiffOptions =
-    GeoTiffOptions(storageMethod, DEFAULT.compression, DEFAULT.colorSpace)
+    DEFAULT.copy(storageMethod = storageMethod)
 
   /**
-    * Creates a new instance of [[GeoTiffOptions]] with the given
-    * Compression and the default [[StorageMethod]] value
-    */
+   * Creates a new instance of [[GeoTiffOptions]] with the given
+   * Compression and the default [[StorageMethod]] value
+   */
   def apply(compression: Compression): GeoTiffOptions =
-    GeoTiffOptions(DEFAULT.storageMethod, compression, DEFAULT.colorSpace)
+    DEFAULT.copy(compression = compression)
+
+  /**
+   * Creates a new instance of [[GeoTiffOptions]] with the given color map.
+   */
+  def apply(colorMap: IndexedColorMap): GeoTiffOptions =
+    DEFAULT.copy(colorSpace = ColorSpace.Palette, colorMap = Some(colorMap))
 }

raster/src/main/scala/geotrellis/raster/io/geotiff/reader/GeoTiffReader.scala

@@ -19,17 +19,16 @@ package geotrellis.raster.io.geotiff.reader
 import geotrellis.raster._
 import geotrellis.raster.io.geotiff._
 import geotrellis.raster.io.geotiff.compression._
-import geotrellis.raster.io.geotiff.util._
 import geotrellis.raster.io.geotiff.tags._
 import geotrellis.vector.Extent
 import geotrellis.proj4.CRS
-import geotrellis.util.{Filesystem, ByteReader}
-
+import geotrellis.util.{ByteReader, Filesystem}
 import monocle.syntax.apply._
-
-import scala.io._
 import java.nio.{ByteBuffer, ByteOrder}
 
+import geotrellis.raster.io.geotiff.tags.codes.ColorSpace
+import geotrellis.raster.render.{ColorMap, IndexedColorMap, RGB}
+
 class MalformedGeoTiffException(msg: String) extends RuntimeException(msg)
 
 class GeoTiffReaderLimitationException(msg: String)
@@ -343,11 +342,15 @@ object GeoTiffReader {
 
     val colorSpace = tiffTags.basicTags.photometricInterp
 
+    val colorMap = if (colorSpace == ColorSpace.Palette && tiffTags.basicTags.colorMap.nonEmpty) {
+      Option(IndexedColorMap.fromTiffPalette(tiffTags.basicTags.colorMap))
+    } else None
+
     GeoTiffInfo(
       tiffTags.extent,
       tiffTags.crs,
       tiffTags.tags,
-      GeoTiffOptions(storageMethod, compression, colorSpace),
+      GeoTiffOptions(storageMethod, compression, colorSpace, colorMap),
       bandType,
       segmentBytes,
       decompressor,
@@ -358,4 +361,6 @@ object GeoTiffReader {
       noDataValue
     )
   }
+
+
 }

raster/src/main/scala/geotrellis/raster/io/geotiff/reader/TiffTagsReader.scala

@@ -412,7 +412,9 @@ object TiffTagsReader {
       if ((tiffTags &|->
         TiffTags._basicTags ^|->
         BasicTags._photometricInterp get) == 3) {
-        val divider = shorts.size / 3
+        // In GDAL world, `divider` ends up being the same as `bitsPerSample`
+        // but theoretically it's valid to have color tables that are smaller
+        val divider = shorts.length / 3
 
         val arr = Array.ofDim[(Short, Short, Short)](divider)
         cfor(0)(_ < divider, _ + 1) { i =>

raster/src/main/scala/geotrellis/raster/io/geotiff/writer/GeoTiffWriter.scala

@@ -16,23 +16,16 @@
 
 package geotrellis.raster.io.geotiff.writer
 
-import geotrellis.raster._
-import geotrellis.raster.io._
 import geotrellis.raster.io.geotiff._
-import geotrellis.vector.Extent
-import geotrellis.proj4.CRS
 
-import geotrellis.raster.io.geotiff.tags.codes._
-import scala.collection.mutable
+import spire.syntax.cfor._
 
 import java.io.ByteArrayOutputStream
 import java.io.DataOutputStream
 import java.io.File
 import java.io.FileOutputStream
 import java.nio.ByteOrder
 
-import spire.syntax.cfor._
-
 object GeoTiffWriter {
   def write(geoTiff: GeoTiffData, path: String): Unit = {
     val fos = new FileOutputStream(new File(path))
@@ -202,3 +195,14 @@ class GeoTiffWriter(geoTiff: GeoTiffData, dos: DataOutputStream) {
     dos.flush()
   }
 }
+
+/**
+ * This exception may be thrown by [[GeoTiffWriter]] in the case where a combination of color space,
+ * color map, and sample depth are not supported by the GeoTiff specification. A specific case is
+ * when [[GeoTiffOptions.colorSpace]] is set to [[geotrellis.raster.io.geotiff.tags.codes.ColorSpace.Palette]]
+ * and [[GeoTiffOptions.colorMap]] is `None` and/or the raster's [[geotrellis.raster.CellType]] is not an 8-bit
+ * or 16-bit integral value.
+ */
+class IncompatibleGeoTiffOptionsException(msg: String, cause: Throwable) extends RuntimeException(msg, cause) {
+  def this(msg: String) = this(msg, null)
+}