Handle TIFFs with more than one image

src/GeoTIFF.jl

@@ -4,7 +4,7 @@
 
 module GeoTIFF
 
-using TiffImages: AbstractTIFF, DenseTaggedImage, WidePixel
+using TiffImages: AbstractTIFF, DenseTaggedImage, StridedTaggedImage, IFD, WidePixel
 using ColorTypes: Colorant, Gray
 using MappedArrays: mappedarray
 using StaticArrays: SVector, SMatrix, SA

src/image.jl

@@ -15,7 +15,7 @@ to get the metadata and the image with corrected axes, respectively.
 * The [`GeoTIFF.image`](@ref) function is necessary because 
   the GeoTIFF format swaps the order of the image axes;
 """
-struct GeoTIFFImage{T,N,I<:AbstractTIFF{T,N}} <: AbstractArray{T,N}
+struct GeoTIFFImage{T,I<:AbstractMatrix{T}} <: AbstractMatrix{T}
   tiff::I
   metadata::Metadata
 end
@@ -59,4 +59,34 @@ channel(geotiff::GeoTIFFImage, i) = mappedarray(c -> channel(c, i), image(geotif
 Base.size(geotiff::GeoTIFFImage) = size(geotiff.tiff)
 Base.getindex(geotiff::GeoTIFFImage, i...) = getindex(geotiff.tiff, i...)
 Base.setindex!(geotiff::GeoTIFFImage, v, i...) = setindex!(geotiff.tiff, v, i...)
-Base.IndexStyle(::Type{GeoTIFFImage{T,N,I}}) where {T,N,I} = IndexStyle(I)
+Base.IndexStyle(::Type{GeoTIFFImage{T,I}}) where {T,I} = IndexStyle(I)
+
+"""
+    GeoTIFF.GeoTIFFImageIterator 
+
+Iterator of [`GeoTIFF.GeoTIFFImage`](@ref) returned by the [`GeoTIFF.load`](@ref) function.
+"""
+struct GeoTIFFImageIterator{I,M}
+  tiffs::I
+  metadata::M
+end
+
+Base.length(geotiffs::GeoTIFFImageIterator) = length(geotiffs.tiffs)
+
+function Base.iterate(geotiffs::GeoTIFFImageIterator)
+  tiff, stateₜ = iterate(geotiffs.tiffs)
+  metadata, stateₘ = iterate(geotiffs.metadata)
+  GeoTIFFImage(tiff, metadata), (stateₜ, stateₘ)
+end
+
+function Base.iterate(geotiffs::GeoTIFFImageIterator, (stateₜ, stateₘ))
+  valueₜ = iterate(geotiffs.tiffs, stateₜ)
+  valueₘ = iterate(geotiffs.metadata, stateₘ)
+  if !isnothing(valueₜ) && !isnothing(valueₘ)
+    tiff, newstateₜ = valueₜ
+    metadata, newstateₘ = valueₘ
+    GeoTIFFImage(tiff, metadata), (newstateₜ, newstateₘ)
+  else
+    nothing
+  end
+end

src/load.jl

@@ -6,8 +6,11 @@
     GeoTIFF.load(fname; kwargs...)
 
 Load a GeoTIFF file returning an image with the processed GeoTIFF metadata.
+If the file contains more than one image, an iterator of images will be returned.
 The keyword arguments are forward to the `TiffImages.load` function.
 
+See also [`GeoTIFF.GeoTIFFImage`](@ref), [`GeoTIFF.GeoTIFFImageIterator`](@ref).
+
 ### Notes
 
 * TIFF files without GeoTIFF metadata are load with default and mandatory metadata:
@@ -16,16 +19,23 @@ The keyword arguments are forward to the `TiffImages.load` function.
 """
 function load(fname; kwargs...)
   tiff = TiffImages.load(fname; kwargs...)
-  metadata = _getmetadata(tiff)
-  GeoTIFFImage(tiff, metadata)
+  ifds = TiffImages.ifds(tiff)
+  metadata = ifds isa IFD ? _getmetadata(ifds) : (_getmetadata(ifd) for ifd in ifds)
+  if tiff isa StridedTaggedImage
+    GeoTIFFImageIterator(tiff, metadata)
+  elseif ndims(tiff) == 3
+    imgs = eachslice(tiff, dims=3)
+    GeoTIFFImageIterator(imgs, metadata)
+  else
+    GeoTIFFImage(tiff, metadata)
+  end
 end
 
 # -----------------
 # HELPER FUNCTIONS
 # -----------------
 
-function _getmetadata(tiff)
-  ifd = TiffImages.ifds(tiff)
+function _getmetadata(ifd)
   geokeydirectory = _gettag(ifd, GeoKeyDirectoryTag, GeoKeyDirectory)
   geodoubleparams = _gettag(ifd, GeoDoubleParamsTag, GeoDoubleParams)
   geoasciiparams = _gettag(ifd, GeoAsciiParamsTag, GeoAsciiParams)

test/runtests.jl

@@ -121,6 +121,67 @@ savedir = mktempdir()
     @test GeoTIFF.channel(geotiff, 4) == channel4
   end
 
+  @testset "GeoTIFFImageIterator" begin
+    # the "test3D.tif" is a file with random data created using TiffImages.jl
+    geotiffitr = GeoTIFF.load(joinpath(datadir, "test3D.tif"))
+    @test geotiffitr isa GeoTIFF.GeoTIFFImageIterator
+    @test length(geotiffitr) == 2
+    geotiff1, state1 = iterate(geotiffitr)
+    geotiff2, state2 = iterate(geotiffitr, state1)
+    @test isnothing(iterate(geotiffitr, state2))
+
+    metadata = GeoTIFF.metadata(geotiff1)
+    @test eltype(geotiff1) <: Gray
+    @test size(geotiff1) == (50, 50)
+    @test GeoTIFF.rastertype(metadata) == GeoTIFF.PixelIsArea
+    @test GeoTIFF.modeltype(metadata) == GeoTIFF.Geographic2D
+    @test GeoTIFF.epsgcode(metadata) == 4326
+    A = [7.2 0.0; 0.0 -3.6]
+    b = [-180.0, 90.0]
+    @test GeoTIFF.affineparams2D(metadata) == (A, b)
+
+    metadata = GeoTIFF.metadata(geotiff2)
+    @test eltype(geotiff2) <: Gray
+    @test size(geotiff2) == (50, 50)
+    @test GeoTIFF.rastertype(metadata) == GeoTIFF.PixelIsArea
+    @test GeoTIFF.modeltype(metadata) == GeoTIFF.Geographic2D
+    @test GeoTIFF.epsgcode(metadata) == 4326
+    A = [7.2 0.0; 0.0 -3.6]
+    b = [-180.0, 90.0]
+    @test GeoTIFF.affineparams2D(metadata) == (A, b)
+
+    # the "ca_nrc_CRD27_00.tif" file is from PROJ CDN
+    # link: https://cdn.proj.org/ca_nrc_CRD27_00.tif
+    geotiffitr = GeoTIFF.load(joinpath(datadir, "ca_nrc_CRD27_00.tif"))
+    @test geotiffitr isa GeoTIFF.GeoTIFFImageIterator
+    @test length(geotiffitr) == 2
+    geotiff1, state1 = iterate(geotiffitr)
+    geotiff2, state2 = iterate(geotiffitr, state1)
+    @test isnothing(iterate(geotiffitr, state2))
+
+    metadata = GeoTIFF.metadata(geotiff1)
+    @test eltype(geotiff1) <: TiffImages.WidePixel
+    @test size(geotiff1) == (46, 91)
+    @test size(GeoTIFF.image(geotiff1)) == (91, 46)
+    @test GeoTIFF.rastertype(metadata) == GeoTIFF.PixelIsPoint
+    @test GeoTIFF.modeltype(metadata) == GeoTIFF.Geographic2D
+    @test GeoTIFF.epsgcode(metadata) == 4267
+    A = [0.016666666666666666 0.0; 0.0 -0.016666666666666666]
+    b = [-124.5, 49.0]
+    @test GeoTIFF.affineparams2D(metadata) == (A, b)
+
+    metadata = GeoTIFF.metadata(geotiff2)
+    @test eltype(geotiff2) <: TiffImages.WidePixel
+    @test size(geotiff2) == (73, 69)
+    @test size(GeoTIFF.image(geotiff2)) == (69, 73)
+    @test GeoTIFF.rastertype(metadata) == GeoTIFF.PixelIsPoint
+    @test GeoTIFF.modeltype(metadata) == GeoTIFF.Geographic2D
+    @test GeoTIFF.epsgcode(metadata) == 4267
+    A = [0.004166666666666667 0.0; 0.0 -0.004166666666666667]
+    b = [-123.53333333333333, 48.7]
+    @test GeoTIFF.affineparams2D(metadata) == (A, b)
+  end
+
   @testset "save" begin
     # saving geotiffs
     file1 = joinpath(datadir, "test.tif")