georaster icon indicating copy to clipboard operation
georaster copied to clipboard

Current image after opening geotiff

Open michaelkirk opened this issue 4 months ago • 3 comments

It seems like this crate doesn't currently support compression. (Correct me if I'm wrong!)

$ cargo run --example pixel data/tiff/N265E425.tif 2550 3050  
<NoData>

Compare this to gdallocationinfo

$ gdallocationinfo data/tiff/N265E425.tif 2550 3050   
Report:
  Location: (2550P,3050L)
  Band 1:
    Value: 2089

If I manually uncompress the geotiff first, I get the expected answer:

$ gdal_translate -co COMPRESS=NONE data/tiff/N265E425.tif data/tiff/N265E425-uncompressed.tif

$ cargo run --example pixel data/tiff/N265E425-uncompressed.tif 2550 3050  
2089

michaelkirk avatar Jul 28 '25 19:07 michaelkirk

I would have expected support for compressed tiffs, since it depends on the tiff crate for reading. tiff has deflate, jpeg and lzw compresseion as default features. But maybe there is something to do on georaster side as well?

pka avatar Jul 29 '25 13:07 pka

I think I figured out what was going on. That particular file includes multiple overviews, and the GeoTiffReader::open command leaves the GeoTiffReader at the last image. In this case, thats the 312x312 one, which is why I'm seeing NODATA.

$ gdalinfo data/tiff/N265E425.tif

Driver: GTiff/GeoTIFF
Files: data/tiff/N265E425.tif
Size is 5000, 5000
Coordinate System is:
PROJCRS["ETRS89_ETRS_LAEA",
    BASEGEOGCRS["ETRS89",
        DATUM["European Terrestrial Reference System 1989",
            ELLIPSOID["GRS 1980",6378137,298.257222101004,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["Lambert Azimuthal Equal Area",
        METHOD["Lambert Azimuthal Equal Area",
            ID["EPSG",9820]],
        PARAMETER["Latitude of natural origin",52,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",10,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",4321000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",3210000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]]]
Data axis to CRS axis mapping: 1,2
Origin = (4250000.000000000000000,2700000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  LAYOUT=COG
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 4250000.000, 2700000.000) (  9d 3'35.67"E, 47d24'34.38"N)
Lower Left  ( 4250000.000, 2650000.000) (  9d 4' 4.68"E, 46d57'34.17"N)
Upper Right ( 4300000.000, 2700000.000) (  9d43'18.96"E, 47d24'47.70"N)
Lower Right ( 4300000.000, 2650000.000) (  9d43'27.54"E, 46d57'47.36"N)
Center      ( 4275000.000, 2675000.000) (  9d23'36.73"E, 47d11'12.75"N)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Gray
  NoData Value=0
  Overviews: 2500x2500, 1250x1250, 625x625, 312x312

I was expecting my newly opened reader to be at the first image. Do you agree that it would be a better default?

In the meanwhile, the work around for me is to add a line manually rewinding the index:

  let mut decoder = GeoTiffReader::open(file).unwrap();
+ decoder.seek_to_image(0).unwrap();
  decoder.read_pixel(2550, 3050)

michaelkirk avatar Jul 30 '25 23:07 michaelkirk

To be clear, my confusion about compression was a red herring. It just so happened that the file that had compression also had multiple layers.

And my test to convert the compressed tif to an uncompressed tif had the side effect of removing the overview images, so it "fixed" things for me.

If I explicitly include the overview images into my uncompressed tif, I hit the exact same error:

$ gdal_translate -co COMPRESS=NONE -co COPY_SRC_OVERVIEWS=YES data/tiff/N265E425.tif data/tiff/N265E425-uncompressed-overviews.tif

$ cargo run --example pixel data/tiff/N265E425-uncompressed-overviews.tif 2550 3050 
<NoData>

So my problem isn't about compression - it's about where the 'current' image should be set after opening a geotiff.

michaelkirk avatar Jul 30 '25 23:07 michaelkirk