Merge pull request #231 from s22s/fix/gtiff-distort-224

Refactor construction of LayoutDefinition in GeoTiff writing process

Author: vpipkt

core/src/main/scala/org/locationtech/rasterframes/expressions/aggregates/ProjectedLayerMetadataAggregate.scala

@@ -28,7 +28,7 @@ import org.locationtech.rasterframes.model.TileDimensions
 import geotrellis.proj4.{CRS, Transform}
 import geotrellis.raster._
 import geotrellis.raster.reproject.{Reproject, ReprojectRasterExtent}
-import geotrellis.spark.tiling.{FloatingLayoutScheme, LayoutLevel}
+import geotrellis.spark.tiling.LayoutDefinition
 import geotrellis.spark.{KeyBounds, SpatialKey, TileLayerMetadata}
 import geotrellis.vector.Extent
 import org.apache.spark.sql.expressions.{MutableAggregationBuffer, UserDefinedAggregateFunction}
@@ -76,7 +76,10 @@ class ProjectedLayerMetadataAggregate(destCRS: CRS, destDims: TileDimensions) ex
   override def evaluate(buffer: Row): Any = {
     import org.locationtech.rasterframes.encoders.CatalystSerializer._
     val buf = buffer.to[BufferRecord]
-    val LayoutLevel(_, layout) = FloatingLayoutScheme(destDims.cols, destDims.rows).levelFor(buf.extent, buf.cellSize)
+
+    val re = RasterExtent(buf.extent, buf.cellSize)
+    val layout = LayoutDefinition(re, destDims.cols, destDims.rows)
+
     val kb = KeyBounds(layout.mapTransform(buf.extent))
     TileLayerMetadata(buf.cellType, layout, buf.extent, destCRS, kb).toRow
   }

core/src/test/scala/org/locationtech/rasterframes/ExtensionMethodSpec.scala

@@ -25,7 +25,7 @@ import geotrellis.proj4.LatLng
 import geotrellis.raster.{ByteCellType, GridBounds, TileLayout}
 import geotrellis.spark.tiling.{CRSWorldExtent, LayoutDefinition}
 import geotrellis.spark.{KeyBounds, SpatialKey, TileLayerMetadata}
-import org.apache.spark.sql.{Encoder, Encoders}
+import org.apache.spark.sql.Encoders
 import org.locationtech.rasterframes.util.SubdivideSupport
 
 /**

core/src/test/scala/org/locationtech/rasterframes/TestData.scala

@@ -146,6 +146,7 @@ trait TestData {
 
   lazy val localSentinel: URI = getClass.getResource("/B01.jp2").toURI
   lazy val cogPath: URI = getClass.getResource("/LC08_RGB_Norfolk_COG.tiff").toURI
+  lazy val singlebandCogPath: URI = getClass.getResource("/LC08_B7_Memphis_COG.tiff").toURI
   lazy val nonCogPath: URI = getClass.getResource("/L8-B8-Robinson-IL.tiff").toURI
 
   lazy val l8B1SamplePath: URI = l8SamplePath(1)

datasource/src/main/scala/org/locationtech/rasterframes/datasource/geotiff/GeoTiffDataSource.scala

@@ -121,6 +121,7 @@ class GeoTiffDataSource
             rf_dimensions(tileCol)
           ))
         .first()
+      logger.debug(s"Contructed TileLayerMetadata: ${tlm.toString}")
 
       val c = ProjectedRasterDefinition(tlm)
 

datasource/src/test/scala/org/locationtech/rasterframes/datasource/geotiff/GeoTiffDataSourceSpec.scala

@@ -23,6 +23,9 @@ package org.locationtech.rasterframes.datasource.geotiff
 import java.nio.file.Paths
 
 import geotrellis.proj4._
+import geotrellis.raster.ProjectedRaster
+import geotrellis.raster.io.geotiff.{MultibandGeoTiff, SinglebandGeoTiff}
+import geotrellis.vector.Extent
 import org.locationtech.rasterframes._
 import org.apache.spark.sql.functions._
 import org.locationtech.rasterframes.TestEnvironment
@@ -41,7 +44,7 @@ class GeoTiffDataSourceSpec
       assert(rf.count() > 10)
     }
 
-    it("should lay out tiles correctly"){
+    it("should lay out tiles correctly") {
 
       val rf = spark.read.format("geotiff").load(cogPath.toASCIIString).asLayer
 
@@ -83,6 +86,9 @@ class GeoTiffDataSourceSpec
 
       assert(result === expected)
     }
+  }
+
+  describe("GeoTiff writing") {
 
     it("should write GeoTIFF RF to parquet") {
       val rf = spark.read.format("geotiff").load(cogPath.toASCIIString).asLayer
@@ -101,6 +107,81 @@ class GeoTiffDataSourceSpec
       }
     }
 
+    it("should write unstructured raster") {
+      import spark.implicits._
+      val df = spark.read.format("raster")
+        .option("tileDimensions", "32,32")  // oddball
+        .load(nonCogPath.toASCIIString) // core L8-B8-Robinson-IL.tiff
+
+      df.count() should be > 0L
+
+      val crs = df.select(rf_crs($"proj_raster")).first()
+
+      noException shouldBe thrownBy {
+        df.write.geotiff.withCRS(crs).save("unstructured.tif")
+      }
+
+      val (inCols, inRows) = {
+        val id = sampleGeoTiff.imageData // inshallah same as nonCogPath
+        (id.cols, id.rows)
+      }
+      inCols should be (774)
+      inRows should be (500) //from gdalinfo
+
+      val outputTif = SinglebandGeoTiff("unstructured.tif")
+      outputTif.imageData.cols should be (inCols)
+      outputTif.imageData.rows should be (inRows)
+
+      // TODO check datatype, extent.
+    }
+
+    it("should round trip unstructured raster from COG"){
+      import spark.implicits._
+      import org.locationtech.rasterframes.datasource.raster._
+
+      val df = spark.read.raster.withTileDimensions(64, 64).load(singlebandCogPath.toASCIIString)
+
+      val resourceCols = 963 // from gdalinfo
+      val resourceRows = 754
+      val resourceExtent = Extent(752325.0, 3872685.0, 781215.0, 3895305.0)
+
+      df.count() should be > 0L
+
+      val crs = df.select(rf_crs(col("proj_raster"))).first()
+
+      val totalExtentRow = df.select(rf_extent($"proj_raster").alias("ext"))
+        .agg(
+          min($"ext.xmin").alias("xmin"),
+          min($"ext.ymin").alias("ymin"),
+          max($"ext.xmax").alias("xmax"),
+          max($"ext.ymax").alias("ymax")
+        ).first()
+      val dfExtent = Extent(totalExtentRow.getDouble(0), totalExtentRow.getDouble(1), totalExtentRow.getDouble(2), totalExtentRow.getDouble(3))
+      logger.info(s"Dataframe extent: ${dfExtent.toString()}")
+
+      dfExtent shouldBe (resourceExtent)
+
+      noException shouldBe thrownBy {
+        df.write.geotiff.withCRS(crs).save("target/unstructured_cog.tif")
+      }
+
+      val (inCols, inRows, inExtent, inCellType) = {
+        val tif = readSingleband("LC08_B7_Memphis_COG.tiff")
+        val id = tif.imageData
+        (id.cols, id.rows, tif.extent, tif.cellType)
+      }
+      inCols should be (963)
+      inRows should be (754) //from gdalinfo
+      inExtent should be (resourceExtent)
+
+      val outputTif = SinglebandGeoTiff("target/unstructured_cog.tif")
+      outputTif.imageData.cols should be (inCols)
+      outputTif.imageData.rows should be (inRows)
+      outputTif.extent should be (resourceExtent)
+      outputTif.cellType should be (inCellType)
+
+    }
+
     it("should write GeoTIFF without layer") {
       import org.locationtech.rasterframes.datasource.raster._
       val pr = col("proj_raster_b0")
@@ -152,6 +233,8 @@ ${s(1)},${s(4)},${s(3)}
         .withDimensions(256, 256)
         .save("geotiff-overview.tif")
 
+      val outTif = MultibandGeoTiff("geotiff-overview.tif")
+      outTif.bandCount should be (3)
     }
   }
 }

datasource/src/test/scala/org/locationtech/rasterframes/datasource/raster/RasterSourceDataSourceSpec.scala

@@ -245,4 +245,41 @@ class RasterSourceDataSourceSpec extends TestEnvironment with TestData {
 
     }
   }
+
+  describe("RasterSource breaks up scenes into tiles") {
+    val modis_df = spark.read.raster
+      .withTileDimensions(128, 128)
+      .withLazyTiles(true)
+      .load(remoteMODIS.toASCIIString())
+
+    val l8_df = spark.read.raster
+      .withTileDimensions(32, 33)
+      .withLazyTiles(true)
+      .load(remoteL8.toASCIIString())
+
+    ignore("should have at most four tile dimensions reading MODIS; ignore until fix #242") {
+      val dims = modis_df.select(rf_dimensions($"proj_raster")).distinct().collect()
+      dims.length should be > (0)
+      dims.length should be <= (4)
+    }
+
+    it("should have at most four tile dimensions reading landsat") {
+      val dims = l8_df.select(rf_dimensions($"proj_raster")).distinct().collect()
+      dims.length should be > (0)
+      dims.length should be <= (4)
+    }
+
+    ignore("should have consistent tile resolution reading MODIS; ignore until fix #242") {
+      val res = modis_df.select((rf_extent($"proj_raster").getField("xmax") - rf_extent($"proj_raster").getField("xmin")) /
+        rf_dimensions($"proj_raster").getField("cols")).distinct().collect()
+      res.length should be (1)
+    }
+
+    it("should have consistent tile resolution reading Landsat") {
+      val res = l8_df.select((rf_extent($"proj_raster").getField("xmax") - rf_extent($"proj_raster").getField("xmin")) /
+        rf_dimensions($"proj_raster").getField("cols")).distinct().collect()
+      res.length should be (1)
+    }
+
+  }
 }

pyrasterframes/src/main/python/setup.py

@@ -153,7 +153,7 @@ def initialize_options(self):
         'shapely',
         'pandas',
         'rasterio',
-        'boto3'
+        'boto3',
     ],
     packages=[
         'pyrasterframes',

pyrasterframes/src/main/python/tests/GeoTiffWriterTests.py

@@ -29,7 +29,7 @@ class GeoTiffWriter(TestEnvironment):
 
     @staticmethod
     def _tmpfile():
-        return os.path.join(tempfile.tempdir, "pyrf-test.tif")
+        return os.path.join(tempfile.gettempdir(), "pyrf-test.tif")
 
     def test_identity_write(self):
         rf = self.spark.read.geotiff(self.img_uri)
@@ -44,13 +44,19 @@ def test_identity_write(self):
 
     def test_unstructured_write(self):
         rf = self.spark.read.raster(self.img_uri)
-        dest = self._tmpfile()
-        rf.write.geotiff(dest, crs='EPSG:32616')
+        dest_file = self._tmpfile()
+        rf.write.geotiff(dest_file, crs='EPSG:32616')
 
-        rf2 = self.spark.read.raster('file://' + dest)
+        rf2 = self.spark.read.raster('file://' + dest_file)
         self.assertEqual(rf2.count(), rf.count())
 
-        os.remove(dest)
+        with rasterio.open(self.img_uri) as source:
+            with rasterio.open(dest_file) as dest:
+                self.assertEqual((dest.width, dest.height), (source.width, source.height))
+                self.assertEqual(dest.bounds, source.bounds)
+                self.assertEqual(dest.crs, source.crs)
+
+        os.remove(dest_file)
 
     def test_downsampled_write(self):
         rf = self.spark.read.raster(self.img_uri)