Issue reading Extent from GeoTiff - offset encountered

[lossyrob]

See the failing test on this branch: https://github.com/lossyrob/geotrellis/tree/bugreport/offset-geotiff-extent

[fosskers]

Is this still relevant?

[lossyrob]

Potentially, can you add the test to master and see if it fails? lossyrob@d6d82ee

[fosskers]

Sure.

[fosskers]

On master with your commit cherry-picked:

[info] - Fails to read correct extent *** FAILED ***
[info]   Extent(1.0, 0.0, 11.0, 10.0) was not equal to Extent(0.5, 0.5, 10.5, 10.5) (FailingTestCase.scala:36)

[fosskers]

From gdalinfo:

Corner Coordinates:
Upper Left  (   0.5000000,  10.5000000) (  0d 0' 0.02"E,  0d 0' 0.34"N)
Lower Left  (   0.5000000,   0.5000000) (  0d 0' 0.02"E,  0d 0' 0.02"N)
Upper Right (  10.5000000,  10.5000000) (  0d 0' 0.34"E,  0d 0' 0.34"N)
Lower Right (  10.5000000,   0.5000000) (  0d 0' 0.34"E,  0d 0' 0.02"N)
Center      (   5.5000000,   5.5000000) (  0d 0' 0.18"E,  0d 0' 0.18"N)
Band 1 Block=10x1 Type=Float32, ColorInterp=Gray

[fosskers]

Yellow brick road:

• Tests call this: https://github.com/geotrellis/geotrellis/blob/master/raster/src/main/scala/geotrellis/raster/io/geotiff/SinglebandGeoTiff.scala#L77
• Eventually gets here: https://github.com/geotrellis/geotrellis/blob/master/raster/src/main/scala/geotrellis/raster/io/geotiff/reader/GeoTiffReader.scala#L83
• Then here: https://github.com/geotrellis/geotrellis/blob/master/raster/src/main/scala/geotrellis/raster/io/geotiff/reader/GeoTiffReader.scala#L240
• Extents are pulled right from tifftags: https://github.com/geotrellis/geotrellis/blob/master/raster/src/main/scala/geotrellis/raster/io/geotiff/reader/GeoTiffReader.scala#L255

[fosskers]

The major Haskell image IO library couldn't parse the GeoTIFF. tiffinfo yields:

TIFFReadDirectory: Warning, Unknown field with tag 33550 (0x830e) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 33922 (0x8482) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 34735 (0x87af) encountered.
TIFFReadDirectory: Warning, Unknown field with tag 34736 (0x87b0) encountered.
TIFF Directory at offset 0x198 (408)
  Image Width: 10 Image Length: 10
  Bits/Sample: 32
  Sample Format: IEEE floating point
  Compression Scheme: None
  Photometric Interpretation: min-is-black
  Orientation: row 0 top, col 0 lhs
  Samples/Pixel: 1
  Rows/Strip: 1
  SMin Sample Value: 0
  SMax Sample Value: 10
  Planar Configuration: single image plane
  Tag 33550: 1.000000,1.000000,1.000000
  Tag 33922: 0.000000,0.000000,0.000000,1.000000,10.000000,0.000000
  Tag 34735: 1,1,0,26,1024,0,1,1,1025,0,1,2,2054,0,1,9102,2057,34736,1,0,2058,34736,1,1,3072,0,1,3857,3075,0,1,7,3076,0,1,9001,3078,34736,1,2,3079,34736,1,3,3080,34736,1,4,3081,34736,1,5,3082,34736,1,6,3083,34736,1,7,3084,34736,1,8,3085,34736,1,9,3086,34736,1,10,3087,34736,1,11,3088,34736,1,12,3089,34736,1,13,3090,34736,1,14,3091,34736,1,15,3092,34736,1,16,3094,34736,1,17,3095,34736,1,18,4099,0,1,9001
  Tag 34736: 6378137.000000,6378137.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,1.000000,0.000000,0.000000

[fosskers]

Debugged from TiffTags.tiePointsModelSpace:

PIXEL SCALES (1.0,1.0,1.0)
PIXEL3D: Pixel3D(0.0,10.0,0.0)
 --> FIRST: Pixel3D(0.0,0.0,0.0), SECOND: Pixel3D(1.0,10.0,0.0)
 --> TRANSLATED PIXEL: Pixel3D(1.0,0.0,0.0)
PIXEL3D: Pixel3D(10.0,0.0,0.0)
 --> FIRST: Pixel3D(0.0,0.0,0.0), SECOND: Pixel3D(1.0,10.0,0.0)
 --> TRANSLATED PIXEL: Pixel3D(11.0,10.0,0.0)

[fosskers]

More importantly:

TIE POINTS FROM TAGS: (Pixel3D(0.0,0.0,0.0),Pixel3D(1.0,10.0,0.0))

Whatever these "tie points" are, they're coming right out of the tags as whole numbers.

That's consistent with what tiffinfo showed:

Tag 33922: 0.000000,0.000000,0.000000,1.000000,10.000000,0.000000

[fosskers]

From what I can see of the code that's here, no one is misbehaving.

  private def tiePointsModelSpace(tiePoints: Array[(Pixel3D, Pixel3D)],
    pixelScaleOption: Option[(Double, Double, Double)]) =
    pixelScaleOption match {
      case Some(pixelScales) => {
        def modelFunc(pixel: Pixel3D) = {
          val (first, second) = tiePoints.head

          val scaleX = (pixel.x - first.x) * pixelScales._1
          val scaleY = (pixel.y - first.y) * pixelScales._2
          val scaleZ = (pixel.z - first.z) * pixelScales._3

          Pixel3D(scaleX + second.x, second.y - scaleY, scaleZ + second.z)
        }

        getExtentFromModelFunction(modelFunc)
      }
...
}

getExtentFromModelFunction calls getRasterBoundaries which asks for image bounds. The length and width are both stored as 10 (although this might be wrong. The tiff info says min value of 0, max value of 10. Technically that's length 11), and the tie points are clearly visible as (0,0,0) and (1,10,0) right in the tiff. So the translation seen above is doing what it's told.

What I'm curious about is how gdalinfo claims the corner coords that it does.

[fosskers]

Notice gdalinfo claims the origin is at Origin = (0.500000000000000,10.500000000000000), but those numbers aren't visible in the tiffinfo output.

[fosskers]

I think we could get the bottom of this together fairly quickly on Monday, @lossyrob .

[lossyrob]

There's model tie points in that raster, and those should be accounted for - somehow this logic here is not accounting for it properly:

geotrellis/raster/src/main/scala/geotrellis/raster/io/geotiff/tags/TiffTags.scala

Line 178 in 9fc2980

def extent: Extent =

[lossyrob]

See https://github.com/OSGeo/gdal/blob/7480b720aa74c87d63d1fbf563373dc7c24fa191/gdal/frmts/gtiff/gt_wkt_srs.cpp#L2722