How to extract location from Strided GeoTIFF raster?

Specific Domains Geo
1

Hello,

I’ve been trying to read some geotiff datasets. Reading it with GeoTIFF.jl results in a GeoTIFF.GeoTIFFImageIterator{TiffImages.StridedTaggedImage{...}}}.
This case is not really documented anywhere neither in GeoTIFF nor in TiffImages so I am unsure on how one is supposed to manipulate it.
This is the full output:

GeoTIFF Output 
GeoTIFF.GeoTIFFImageIterator{TiffImages.StridedTaggedImage{UInt64, Matrix{ColorTypes.Gray{Float32}}}, Base.Generator{Vector{TiffImages.IFD{UInt64}}, GeoTIFF.var"#19#20"}}(Any[ColorTypes.Gray{Float32}[Gray{Float32}(-7.5412407f0) Gray{Float32}(-7.70549f0) … Gray{Float32}(-2375.686f0) Gray{Float32}(-2376.507f0); Gray{Float32}(-8.608861f0) Gray{Float32}(-8.773109f0) … Gray{Float32}(-2374.7825f0) Gray{Float32}(-2376.425f0); … ; Gray{Float32}(43.62238f0) Gray{Float32}(45.018497f0) … Gray{Float32}(-1532.1022f0) Gray{Float32}(-1532.1842f0); Gray{Float32}(45.10062f0) Gray{Float32}(46.16824f0) … Gray{Float32}(-1532.1022f0) Gray{Float32}(-1532.1842f0)], ColorTypes.Gray{Float32}[Gray{Float32}(-8.157152f0) Gray{Float32}(-8.732026f0) … Gray{Float32}(-2377.2463f0) Gray{Float32}(-2375.85f0); Gray{Float32}(-10.087005f0) Gray{Float32}(-10.497636f0) … Gray{Float32}(-2377.6362f0) Gray{Float32}(-2375.727f0); … ; Gray{Float32}(43.66345f0) Gray{Float32}(43.499256f0) … Gray{Float32}(-1532.1227f0) Gray{Float32}(-1532.2664f0); Gray{Float32}(44.977463f0) Gray{Float32}(44.669476f0) … Gray{Float32}(-1532.1022f0) Gray{Float32}(-1532.1432f0)], ColorTypes.Gray{Float32}[Gray{Float32}(-9.368346f0) Gray{Float32}(-10.163881f0) … Gray{Float32}(-2377.9905f0) Gray{Float32}(-2376.6147f0); Gray{Float32}(-11.5189495f0) Gray{Float32}(-11.585778f0) … Gray{Float32}(-2375.5735f0) Gray{Float32}(-2376.148f0); … ; Gray{Float32}(43.396812f0) Gray{Float32}(43.288292f0) … Gray{Float32}(-1532.1022f0) Gray{Float32}(-1532.4409f0); Gray{Float32}(44.202477f0) Gray{Float32}(42.852448f0) … Gray{Float32}(-1532.1022f0) Gray{Float32}(-1532.1586f0)], ColorTypes.Gray{Float32}[Gray{Float32}(-10.659187f0) Gray{Float32}(-11.235324f0) … Gray{Float32}(-2376.605f0) Gray{Float32}(-2376.5815f0); Gray{Float32}(-11.747573f0) Gray{Float32}(-11.647497f0) … Gray{Float32}(-2372.8823f0) Gray{Float32}(-2375.5977f0); … ; Gray{Float32}(43.17711f0) Gray{Float32}(43.858208f0) … Gray{Float32}(-1531.9242f0) Gray{Float32}(-1532.7655f0); Gray{Float32}(43.435093f0) Gray{Float32}(44.234425f0) … Gray{Float32}(-1531.7432f0) Gray{Float32}(-1532.2009f0)], ColorTypes.Gray{Float32}[Gray{Float32}(-11.322267f0) Gray{Float32}(-11.448627f0) … Gray{Float32}(-2370.4438f0) Gray{Float32}(-2375.417f0); Gray{Float32}(-11.569621f0) Gray{Float32}(-11.449963f0) … Gray{Float32}(-2366.1345f0) Gray{Float32}(-2371.188f0); … ; Gray{Float32}(44.226574f0) Gray{Float32}(43.909412f0) … Gray{Float32}(-1531.3268f0) Gray{Float32}(-1532.9606f0); Gray{Float32}(43.67626f0) Gray{Float32}(44.4127f0) … Gray{Float32}(-1530.5358f0) Gray{Float32}(-1532.1583f0)], ColorTypes.Gray{Float32}[Gray{Float32}(-11.4475765f0) Gray{Float32}(-11.208084f0) … Gray{Float32}(-2365.8428f0) Gray{Float32}(-2370.7974f0); Gray{Float32}(-11.140415f0) Gray{Float32}(-11.298808f0) … Gray{Float32}(-2354.5442f0) Gray{Float32}(-2360.42f0); … ; Gray{Float32}(43.584164f0) Gray{Float32}(43.737946f0) … Gray{Float32}(-1529.1589f0) Gray{Float32}(-1533.4852f0); Gray{Float32}(44.05623f0) Gray{Float32}(44.68009f0) … Gray{Float32}(-1528.362f0) Gray{Float32}(-1531.7451f0)]], Base.Generator{Vector{TiffImages.IFD{UInt64}}, GeoTIFF.var"#19#20"}(GeoTIFF.var"#19#20"(), TiffImages.IFD{UInt64}[IFD, with tags:
        Tag(IMAGEWIDTH, 16666)
        Tag(IMAGELENGTH, 23333)
        Tag(BITSPERSAMPLE, 32)
        Tag(COMPRESSION, COMPRESSION_ADOBE_DEFLATE)
        Tag(PHOTOMETRIC, 1)
        Tag(SAMPLESPERPIXEL, 1)
        Tag(PLANARCONFIG, 1)
        Tag(PREDICTOR, 1)
        Tag(TILEWIDTH, 512)
        Tag(TILELENGTH, 512)
        Tag(TILEOFFSETS, UInt64[424318014, 424600685, 424898509, 425059230, 425180156, ...])
        Tag(TILEBYTECOUNTS, UInt32[282663, 297816, 160713, 120918, 132183, ...])
        Tag(SAMPLEFORMAT, 3)
        Tag(MODELPIXELSCALE, Float64[0.0002999999999999996, 0.0003000000000000001, 0.0])
        Tag(MODELTIEPOINT, Float64[0.0, 0.0, 0.0, 142.0000611111, -10.000138888905807, ...])
        Tag(GEOKEYDIRECTORY, UInt16[1, 1, 0, 7, 1024, ...])
        Tag(GEODOUBLEPARAMS, Float64[298.257223563, 6.378137e6])
        Tag(GEOASCIIPARAMS, "WGS 84|")
        Tag(GDALNODATA, "-9999"), IFD, with tags:
        Tag(SUBFILETYPE, 1)
        Tag(IMAGEWIDTH, 8333)
        Tag(IMAGELENGTH, 11667)
        Tag(BITSPERSAMPLE, 32)
        Tag(COMPRESSION, COMPRESSION_ADOBE_DEFLATE)
        Tag(PHOTOMETRIC, 1)
        Tag(SAMPLESPERPIXEL, 1)
        Tag(PLANARCONFIG, 1)
        Tag(PREDICTOR, 1)
        Tag(TILEWIDTH, 512)
        Tag(TILELENGTH, 512)
        Tag(TILEOFFSETS, UInt64[110023558, 110944836, 111707108, 112501107, 113298705, ...])
        Tag(TILEBYTECOUNTS, UInt32[921270, 762264, 793991, 797590, 899336, ...])
        Tag(SAMPLEFORMAT, 3)
        Tag(GDALNODATA, "-9999"), IFD, with tags:
        Tag(SUBFILETYPE, 1)
        Tag(IMAGEWIDTH, 4167)
        Tag(IMAGELENGTH, 5834)
        Tag(BITSPERSAMPLE, 32)
        Tag(COMPRESSION, COMPRESSION_ADOBE_DEFLATE)
        Tag(PHOTOMETRIC, 1)
        Tag(SAMPLESPERPIXEL, 1)
        Tag(PLANARCONFIG, 1)
        Tag(PREDICTOR, 1)
        Tag(TILEWIDTH, 512)
        Tag(TILELENGTH, 512)
        Tag(TILEOFFSETS, UInt64[27895940, 28804616, 29687974, 30614470, 31549359, ...])
        Tag(TILEBYTECOUNTS, UInt32[908668, 883350, 926488, 934881, 886679, ...])
        Tag(SAMPLEFORMAT, 3)
        Tag(GDALNODATA, "-9999"), IFD, with tags:
        Tag(SUBFILETYPE, 1)
        Tag(IMAGEWIDTH, 2084)
        Tag(IMAGELENGTH, 2917)
        Tag(BITSPERSAMPLE, 32)
        Tag(COMPRESSION, COMPRESSION_ADOBE_DEFLATE)
        Tag(PHOTOMETRIC, 1)
        Tag(SAMPLESPERPIXEL, 1)
        Tag(PLANARCONFIG, 1)
        Tag(PREDICTOR, 1)
        Tag(TILEWIDTH, 512)
        Tag(TILELENGTH, 512)
        Tag(TILEOFFSETS, UInt64[6818639, 7745514, 8691403, 9563744, 10427788, ...])
        Tag(TILEBYTECOUNTS, UInt32[926867, 945881, 872333, 864036, 64373, ...])
        Tag(SAMPLEFORMAT, 3)
        Tag(GDALNODATA, "-9999"), IFD, with tags:
        Tag(SUBFILETYPE, 1)
        Tag(IMAGEWIDTH, 1042)
        Tag(IMAGELENGTH, 1459)
        Tag(BITSPERSAMPLE, 32)
        Tag(COMPRESSION, COMPRESSION_ADOBE_DEFLATE)
        Tag(PHOTOMETRIC, 1)
        Tag(SAMPLESPERPIXEL, 1)
        Tag(PLANARCONFIG, 1)
        Tag(PREDICTOR, 1)
        Tag(TILEWIDTH, 512)
        Tag(TILELENGTH, 512)
        Tag(TILEOFFSETS, UInt64[1416326, 2378142, 3259160, 3293416, 4254955, ...])
        Tag(TILEBYTECOUNTS, UInt32[961808, 881010, 34248, 961531, 892386, ...])
        Tag(SAMPLEFORMAT, 3)
        Tag(GDALNODATA, "-9999"), IFD, with tags:
        Tag(SUBFILETYPE, 1)
        Tag(IMAGEWIDTH, 521)
        Tag(IMAGELENGTH, 730)
        Tag(BITSPERSAMPLE, 32)
        Tag(COMPRESSION, COMPRESSION_ADOBE_DEFLATE)
        Tag(PHOTOMETRIC, 1)
        Tag(SAMPLESPERPIXEL, 1)
        Tag(PLANARCONFIG, 1)
        Tag(PREDICTOR, 1)
        Tag(TILEWIDTH, 512)
        Tag(TILELENGTH, 512)
        Tag(TILEOFFSETS, UInt64[27052, 977349, 996614, 1407686])
        Tag(TILEBYTECOUNTS, UInt32[950289, 19257, 411064, 8632])
        Tag(SAMPLEFORMAT, 3)
        Tag(GDALNODATA, "-9999")]))

It seems that the dataset contains the same region but at different resolutions.

I can access the raw elevation data, by iterating manually or defining nth(itr, n) = first(Iterators.drop(itr, n-1)) to access the images, but the extracted images don’t seem to have any metadata.

I would need the spatial information of the dataset, such as resolution, or extents, because I need to extract a cutout, but I can’t seem to find a way to do it with this dataset.

I have also tried with Rasters.jl but to no avail.

Am I using the wrong package? is GeoTIFF.jl too low level?

1 Like
2

What error do you have in Rasters.jl?
I don’t know whether GeoTIFF.jl can handle these overviews but Rasters should be able to handle this format when you load ArchGDAL.jl to load the data using GDAL in the background.

3

@cshen take a look into the test suite of GeoTIFF.jl. We load these iterators all the time in CoordGridTransforms.jl as part of CoordRefSystems.jl. No need for GDAL.

@eliascarv can share more details if you get stuck.

1 Like
4

with

imgs = GeoTIFF.load(file)
img = nth(imgs, 1) # 23333×16666 GeoTIFF.GeoTIFFImage{ColorTypes.Gray{Float32}, Matrix{ColorTypes.Gray{Float32}}}
Raster(img, dims=(23333,16666))

I have

ERROR: MethodError: no method matching _format(::Int64, ::Base.OneTo{Int64})

But I am most likely not using it correctly, so I don’t know how relevant this error is.

That’s seems an excellent idea! Thanks, I’ll see what I can do!

5

Raster(file) if file is the path to a tiff

Trying random arguments is unlikely to succeed :wink:

6

Ok, I think GeoTIFF.affineparams2D was the missing piece I needed.

But interestingly, this works only for the first image.

julia> geoimg = GeoTIFF.load(file)
julia> GeoTIFF.modeltype(GeoTIFF.metadata(nth(geoimg, 1)))
0x0002

julia> GeoTIFF.modeltype(GeoTIFF.metadata(nth(geoimg, 2)))

julia>

Now that I look back at the output when loading the file, there are a lot of repeating tags, that seem associated with the various images at different resolutions.
Is the dataset badly formatted?

7

More options. I just picked one GeoTIFF file at random among those that I have in a work dir.

using GMT

julia> I = gmtread("KR_S2B_MSIL2A_20210614_20.tiff")
A GMTimage object with 10 bands of type UInt16
Pixel node registration used
x_min: 445000.0 x_max :463000.0 x_inc :20.0     n_columns :900
y_min: 5.106e6  y_max :5.124e6  y_inc :20.0     n_rows :900
z_min: 0.0      z_max :17936.0
Mem layout:     BRBa
PROJ: +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs
8

Hi @cshen, I don’t remember the convention out of my head. But the metadata is stored at the top level, including these datasets with iterators. If you encounter any issue with the GeoTIFF.jl implementation, please feel free to suggest improvements.

9

It wasn’t completely random, just a bit.
I saw that Raster(filepath) was available but I didn’t want to bring in the whole GDAL, felt overkill for what I needed to do :slight_smile:

I don’t know the standard, nor the conventions, so I don’t know if this is something that whoever made the file decided or if it can be relied upon, but It would have been awesome, in this case, to have spatial information automatically inferred for all images in the iterator, based on the top level image metadata, and rescaled using the difference in image size.

for example:

julia> GeoTIFF.affineparams2D(GeoTIFF.metadata(nth(geotiffs, 1)))
([0.0002999999999999996 0.0; 0.0 -0.0003000000000000001], [142.0000611111, -10.000138888905807])

julia> GeoTIFF.affineparams2D(GeoTIFF.metadata(nth(geotiffs, 2)))
([1.0 0.0; 0.0 1.0], [0.0, 0.0])

julia> size(nth(geotiffs, 1))
(23333, 16666)

julia> size(nth(geotiffs, 2))
(11667, 8333)

while it could be inferred to be:

modelpixelscale = [0.0006 0.0; 0.0 -0.0006] # rescaled from the size ratio
modeltiepoints = [142.0000611111, -10.000138888905807] #inherited from main image

or something like that.

EDIT: I did a little digging, It is actually possible to know that this is indeed the case using the SUBFILETYPE tag, in the TIFF Spec [0][1]:

A general indication of the kind of data contained in this subfile.

Replaces the old SubfileType field, due to limitations in the definition of that field.

NewSubfileType is mainly useful when there are multiple subfiles in a single TIFF file.

This field is made up of a set of 32 flag bits. Unused bits are expected to be 0. Bit 0 is the low-order bit. Currently defined values are:

Bit 0 is 1 if the image is a reduced-resolution version of another image in this TIFF file; else the bit is 0.
Bit 1 is 1 if the image is a single page of a multi-page image (see the PageNumber field description); else the bit is 0.
Bit 2 is 1 if the image defines a transparency mask for another image in this TIFF file. The PhotometricInterpretation value must be 4, designating a transparency mask.

Which for my file is indeed a chain of IFDs with Tag(SUBFILETYPE, 1)

[0]: Baseline TIFF Tag NewSubfileType (TIFFTAG_SUBFILETYPE), code 254 (0x00FE)
[1]: PDF of TIFF SPEC

1 Like
10

Thanks for digging into the spec @cshen. Feel free to submit a PR with any improvement that is compliant. GeoTIFF.jl should support any spec-compatible file.

11

Sure, I’ll make an issue and setup a draft PR in the next days!

1 Like
12

As a side note, you may be interested in GeoIO.jl, which uses GeoTIFF.jl as the backend to load these datasets as geotables. All explained in the GDSJL book.

If you need a lightweight dependency, then GeoTIFF.jl is the way to go. Thanks for considering a contribution.

13

By the way, it would be awesome to have GeoTIFF.jl instead of GDAL in Rasters.jl :smiley: (lets julia all the things)