[WIP] Add support for GeoTIFFs with multiple grids

src/CoordGridTransforms.jl

@@ -11,7 +11,7 @@ using Unitful: unit, numtype
 using Unitful: uconvert, ustrip
 using StaticArrays: SVector, SA
 using MappedArrays: mappedarray
-using Interpolations: interpolate, Gridded, Linear
+using Interpolations: interpolate, Gridded, Linear, bounds
 using DataDeps: @datadep_str, register, DataDep
 
 import GeoTIFF

src/grids.jl

@@ -12,13 +12,26 @@ function geotiff end
 # cache interpolator objects to avoid interpolating the same grid twice
 const INTERPOLATOR = IdDict()
 
+function interpolatepoint(Datumₛ, Datumₜ, lat, lon)
+  lat′ = ustrip(lat)
+  lon′ = ustrip(lon)
+  interps = interpolators(Datumₛ, Datumₜ)
+  for interp in interps
+    (lonmin, lonmax), (latmin, latmax) = bounds(interp)
+    if lonmin < lon′ < lonmax && latmin < lat′ < latmax
+      return interp(lon′, lat′)
+    end
+  end
+  nothing
+end
+
 """
     interpolator(Datumₛ, Datumₜ)
 
 Linear interpolation of GeoTIFF grid that converts `Datumₛ` to `Datumₜ`.
 All of the GeoTIFF channels are combined into the interpolated grid as a vector.
 """
-function interpolator(Datumₛ, Datumₜ)
+function interpolators(Datumₛ, Datumₜ)
   if haskey(INTERPOLATOR, (Datumₛ, Datumₜ))
     return INTERPOLATOR[(Datumₛ, Datumₜ)]
   end
@@ -27,7 +40,20 @@ function interpolator(Datumₛ, Datumₜ)
   file = downloadgeotiff(Datumₛ, Datumₜ)
   geotiff = GeoTIFF.load(file)
 
-  # construct the interpolation grid
+  # interpolator of each grid
+  interps = if geotiff isa GeoTIFF.GeoTIFFImageIterator
+    map(interpolator, geotiff)
+  else
+    [interpolator(geotiff)]
+  end
+
+  # store interpolators in cache
+  INTERPOLATOR[(Datumₛ, Datumₜ)] = interps
+
+  interps
+end
+
+function interpolator(geotiff)
   img = GeoTIFF.image(geotiff)
   grid = mappedarray(img) do color
     N = GeoTIFF.nchannels(color)
@@ -59,12 +85,8 @@ function interpolator(Datumₛ, Datumₜ)
     grid = @view grid[:, n:-1:1]
   end
 
-  # create the interpolation
-  interp = interpolate((lonrange, latrange), grid, Gridded(Linear()))
-
-  INTERPOLATOR[(Datumₛ, Datumₜ)] = interp
-
-  interp
+  # create the interpolator
+  interpolate((lonrange, latrange), grid, Gridded(Linear()))
 end
 
 """

src/transforms/geocgridtranslation.jl

@@ -48,8 +48,7 @@ function geocgridtranslationparams(Datumₛ, Datumₜ, xyz)
   T = numtype(eltype(xyz))
   # target latlon coordinates for interpolation
   latlon = convert(LatLon, Cartesian{Datumₜ}(xyz...))
-  interp = interpolator(Datumₛ, Datumₜ)
-  itp = interp(ustrip(latlon.lon), ustrip(latlon.lat))
+  itp = interpolatepoint(Datumₛ, Datumₜ, latlon.lat, latlon.lon)
   # type assertion is necessary for type stability
   δx::T = T(itp[1])
   δy::T = T(itp[2])

src/transforms/hgridshift.jl

@@ -24,14 +24,13 @@ macro hgridshift(Datumₛ, Datumₜ)
 end
 
 function hgridshiftfwd(Datumₛ, Datumₜ, (lat, lon))
-  latshift, lonshift = hgridshiftparams(Datumₛ, Datumₜ, (lat, lon))
+  latshift, lonshift = hgridshiftparams(Datumₛ, Datumₜ, lat, lon)
   lat + latshift, lon + lonshift
 end
 
-function hgridshiftparams(Datumₛ, Datumₜ, (lat, lon))
+function hgridshiftparams(Datumₛ, Datumₜ, lat, lon)
   T = numtype(lon)
-  interp = interpolator(Datumₛ, Datumₜ)
-  itp = interp(ustrip(lon), ustrip(lat))
+  itp = interpolatepoint(Datumₛ, Datumₜ, lat, lon)
   # type assertion is necessary for type stability
   latshift::T = T(itp[1])
   lonshift::T = T(itp[2])

src/transforms/pointmotion.jl

@@ -27,7 +27,7 @@ function pointmotionfwd(Datumₛ, Datumₜ, (lat, lon, alt))
   ϕ = ustrip(deg2rad(lat))
   h = ustrip(m, alt)
 
-  λₛ, ϕₛ, hₛ = pointmotionparams(Datumₛ, Datumₜ, lon, lat, ϕ, h)
+  λₛ, ϕₛ, hₛ = pointmotionparams(Datumₛ, Datumₜ, lat, lon, ϕ, h)
   λ′ = λ + λₛ
   ϕ′ = ϕ + ϕₛ
   h′ = h + hₛ
@@ -40,14 +40,13 @@ function pointmotionfwd(Datumₛ, Datumₜ, (lat, lon, alt))
   lat′, lon′, alt′
 end
 
-function pointmotionparams(Datumₛ, Datumₜ, lon, lat, ϕ, h)
+function pointmotionparams(Datumₛ, Datumₜ, lat, lon, ϕ, h)
   🌎 = ellipsoid(Datumₛ)
   T = numtype(lon)
   a = T(ustrip(m, majoraxis(🌎)))
   e² = T(eccentricity²(🌎))
 
-  interp = interpolator(Datumₛ, Datumₜ)
-  itp = interp(ustrip(lon), ustrip(lat))
+  itp = interpolatepoint(Datumₛ, Datumₜ, lat, lon)
   # type assertion is necessary for type stability
   # convert millimeters to meters
   eᵥ::T = T(itp[1]) / 1000