From 0d12b321a764e67cd6f46472c1e1590dbdfa3ef2 Mon Sep 17 00:00:00 2001 From: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> Date: Thu, 7 Nov 2024 17:19:10 -0300 Subject: [PATCH] Add utility functions to easily extract values from metadata (#3) * Add utility functions to easily extract values from metadata * Add 'affineparams2D' * Fix code --- src/GeoTIFF.jl | 1 + src/metadata.jl | 8 ++--- src/userutils.jl | 80 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 85 insertions(+), 4 deletions(-) create mode 100644 src/userutils.jl diff --git a/src/GeoTIFF.jl b/src/GeoTIFF.jl index 813b97e..756eac9 100644 --- a/src/GeoTIFF.jl +++ b/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") diff --git a/src/metadata.jl b/src/metadata.jl index 94dfb91..a52e23e 100644 --- a/src/metadata.jl +++ b/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 diff --git a/src/userutils.jl b/src/userutils.jl new file mode 100644 index 0000000..b4b6a31 --- /dev/null +++ b/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