Add utility functions to easily extract values from metadata

src/GeoTIFF.jl

@@ -13,6 +13,7 @@ import TiffImages
 
 include("geokeys.jl")
 include("metadata.jl")
+include("userutils.jl")
 include("image.jl")
 include("load.jl")
 include("save.jl")

src/metadata.jl

@@ -31,10 +31,10 @@ GeoKeyDirectory(; version=1, revision=1, minor=1, geokeys=GeoKey[]) =
 function GeoKeyDirectory(params::Vector{UInt16})
   nkeys = params[4]
   geokeys = map((1:nkeys) * 4) do i
-    id = GeoKeyID(geokey[1 + i])
-    tag = geokey[2 + i]
-    count = geokey[3 + i]
-    value = geokey[4 + i]
+    id = GeoKeyID(params[1 + i])
+    tag = params[2 + i]
+    count = params[3 + i]
+    value = params[4 + i]
     GeoKey(id, tag, count, value)
   end
 

src/userutils.jl

@@ -0,0 +1,80 @@
+# -----------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# -----------------------------------------------------------------
+
+function geokey(metadata::Metadata, id::GeoKeyID)
+  geokeys = metadata.geokeydirectory.geokeys
+  i = findfirst(geokey -> geokey.id == id, geokeys)
+  isnothing(i) ? nothing : geokeys[i]
+end
+
+function geokeyvalue(metadata::Metadata, id::GeoKeyID)
+  gk = geokey(metadata, id)
+  isnothing(gk) ? nothing : gk.value
+end
+
+rastertype(metadata::Metadata) = geokeyvalue(metadata, GTRasterTypeGeoKey)
+
+modeltype(metadata::Metadata) = geokeyvalue(metadata, GTModelTypeGeoKey)
+
+function epsgcode(metadata::Metadata)
+  mt = modeltype(metadata)
+  isnothing(mt) && return nothing
+  if mt == Projected2D
+    geokeyvalue(metadata, ProjectedCRSGeoKey)
+  elseif mt == Geographic2D || mt == Geocentric3D
+    geokeyvalue(metadata, GeodeticCRSGeoKey)
+  else
+    # Undefined or UserDefined
+    nothing
+  end
+end
+
+function affineparams2D(metadata::Metadata)
+  pixelscale = metadata.modelpixelscale
+  tiepoint = metadata.modeltiepoint
+  transformation = metadata.modeltransformation
+  if !isnothing(pixelscale) && !isnothing(tiepoint)
+    sx = pixelscale.x
+    sy = pixelscale.y
+    (; i, j, x, y) = tiepoint
+    tx = x - i / sx
+    ty = y + j / sy
+    A = [
+      sx 0
+      0 -sy
+    ]
+    b = [tx, ty]
+    A, b
+  elseif !isnothing(transformation)
+    transformation.A[1:2, 1:2], transformation.b[1:2]
+  else
+    nothing
+  end
+end
+
+function affineparams3D(metadata::Metadata)
+  pixelscale = metadata.modelpixelscale
+  tiepoint = metadata.modeltiepoint
+  transformation = metadata.modeltransformation
+  if !isnothing(pixelscale) && !isnothing(tiepoint)
+    sx = pixelscale.x
+    sy = pixelscale.y
+    sz = pixelscale.z
+    (; i, j, k, x, y, z) = tiepoint
+    tx = x - i / sx
+    ty = y + j / sy
+    tz = z - k / sz
+    A = [
+      sx 0 0
+      0 -sy 0
+      0 0 sz
+    ]
+    b = [tx, ty, tz]
+    A, b
+  elseif !isnothing(transformation)
+    transformation.A, transformation.b
+  else
+    nothing
+  end
+end