Initial implementation

Project.toml

@@ -3,5 +3,15 @@ uuid = "c470187e-1474-4b46-a068-a417f109db17"
 authors = ["Elias Carvalho <[email protected]> and contributors"]
 version = "0.1.0"
 
+[deps]
+ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
+MappedArrays = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
+TiffImages = "731e570b-9d59-4bfa-96dc-6df516fadf69"
+
 [compat]
+ColorTypes = "0.11"
+FixedPointNumbers = "0.8"
+MappedArrays = "0.4"
+TiffImages = "0.11"
 julia = "1.10"

src/GeoTIFF.jl

@@ -1,5 +1,19 @@
+# -----------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# -----------------------------------------------------------------
+
 module GeoTIFF
 
-# Write your package code here.
+using TiffImages: AbstractTIFF, DenseTaggedImage, WidePixel
+using ColorTypes: Colorant, Gray
+using MappedArrays: mappedarray
+using FixedPointNumbers: Fixed, Normed
+
+import TiffImages
+
+include("metadata.jl")
+include("image.jl")
+include("load.jl")
+include("save.jl")
 
 end

src/image.jl

@@ -0,0 +1,18 @@
+# -----------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# -----------------------------------------------------------------
+
+struct GeoTIFFImage{T,N,I<:AbstractTIFF{T,N}} <: AbstractArray{T,N}
+  tiff::I
+  metadata::Metadata
+end
+
+tiff(geotiff::GeoTIFFImage) = geotiff.tiff
+
+metadata(geotiff::GeoTIFFImage) = geotiff.metadata
+
+# AbstractArray interface
+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)

src/load.jl

@@ -0,0 +1,30 @@
+# -----------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# -----------------------------------------------------------------
+
+function load(fname; kwargs...)
+  tiff = TiffImages.load(fname; kwargs...)
+  metadata = _getmetadata(tiff)
+  GeoTIFFImage(tiff, metadata)
+end
+
+# -----------------
+# HELPER FUNCTIONS
+# -----------------
+
+function _getmetadata(tiff)
+  ifd = TiffImages.ifds(tiff)
+  geokeydirectory = _gettag(ifd, GeoKeyDirectoryTag, GeoKeyDirectory)
+  geodoubleparams = _gettag(ifd, GeoDoubleParamsTag, GeoDoubleParams)
+  geoasciiparams = _gettag(ifd, GeoAsciiParamsTag, GeoAsciiParams)
+  modelpixelscale = _gettag(ifd, ModelPixelScaleTag, ModelPixelScale)
+  modeltiepoint = _gettag(ifd, ModelTiepointTag, ModelTiepoint)
+  modeltransformation = _gettag(ifd, ModelTransformationTag, ModelTransformation)
+  Metadata(; geokeydirectory, geodoubleparams, geoasciiparams, modelpixelscale, modeltiepoint, modeltransformation)
+end
+
+function _gettag(ifd, tag, Type)
+  params = TiffImages.getdata(ifd, UInt16(tag), nothing)
+  isnothing(params) && return nothing
+  Type(params)
+end

src/metadata.jl

@@ -0,0 +1,214 @@
+# -----------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# -----------------------------------------------------------------
+
+@enum GeoTag::UInt16 begin
+  GeoKeyDirectoryTag = 34735
+  GeoDoubleParamsTag = 34736
+  GeoAsciiParamsTag = 34737
+  ModelPixelScaleTag = 33550
+  ModelTiepointTag = 33922
+  ModelTransformationTag = 34264
+end
+
+@enum GeoKeyID::UInt16 begin
+  GTRasterTypeGeoKey = 1025
+  GTModelTypeGeoKey = 1024
+  ProjectedCRSGeoKey = 3072
+  GeodeticCRSGeoKey = 2048
+  VerticalGeoKey = 4096
+  GTCitationGeoKey = 1026
+  GeodeticCitationGeoKey = 2049
+  ProjectedCitationGeoKey = 3073
+  VerticalCitationGeoKey = 4097
+  GeogAngularUnitsGeoKey = 2054
+  GeogAzimuthUnitsGeoKey = 2060
+  GeogLinearUnitsGeoKey = 2052
+  ProjLinearUnitsGeoKey = 3076
+  VerticalUnitsGeoKey = 4099
+  GeogAngularUnitSizeGeoKey = 2055
+  GeogLinearUnitSizeGeoKey = 2053
+  ProjLinearUnitSizeGeoKey = 3077
+  GeodeticDatumGeoKey = 2050
+  PrimeMeridianGeoKey = 2051
+  PrimeMeridianLongitudeGeoKey = 2061
+  EllipsoidGeoKey = 2056
+  EllipsoidSemiMajorAxisGeoKey = 2057
+  EllipsoidSemiMinorAxisGeoKey = 2058
+  EllipsoidInvFlatteningGeoKey = 2059
+  VerticalDatumGeoKey = 4098
+  ProjectionGeoKey = 3074
+  ProjMethodGeoKey = 3075
+  ProjStdParallel1GeoKey = 3078
+  ProjStdParallel2GeoKey = 3079
+  ProjNatOriginLongGeoKey = 3080
+  ProjNatOriginLatGeoKey = 3081
+  ProjFalseOriginLongGeoKey = 3084
+  ProjFalseOriginLatGeoKey = 3085
+  ProjCenterLongGeoKey = 3088
+  ProjCenterLatGeoKey = 3089
+  ProjStraightVertPoleLongGeoKey = 3095
+  ProjAzimuthAngleGeoKey = 3094
+  ProjFalseEastingGeoKey = 3082
+  ProjFalseNorthingGeoKey = 3083
+  ProjFalseOriginEastingGeoKey = 3086
+  ProjFalseOriginNorthingGeoKey = 3087
+  ProjCenterEastingGeoKey = 3090
+  ProjCenterNorthingGeoKey = 3091
+  ProjScaleAtNatOriginGeoKey = 3092
+  ProjScaleAtCenterGeoKey = 3093
+end
+
+# Corresponding names in the GeoTIFF specification:
+# id - KeyID
+# tag - TIFFTagLocation
+# count - Count
+# value - ValueOffset
+struct GeoKey
+  id::GeoKeyID
+  tag::UInt16
+  count::UInt16
+  value::UInt16
+end
+
+# Corresponding names in the GeoTIFF specification:
+# version - KeyDirectoryVersion
+# revision - KeyRevision
+# minor - MinorRevision
+# nkeys - NumberOfKeys
+# geokeys - Key Entry Set
+struct GeoKeyDirectory
+  version::UInt16
+  revision::UInt16
+  minor::UInt16
+  nkeys::UInt16
+  geokeys::Vector{GeoKey}
+end
+
+GeoKeyDirectory(; version=1, revision=1, minor=1, geokeys=GeoKey[]) =
+  GeoKeyDirectory(version, revision, minor, length(geokeys), geokeys)
+
+function GeoKeyDirectory(params::Vector{UInt16})
+  geokeyvalues = @view params[5:end]
+  geokeys = map(Iterators.partition(geokeyvalues, 4)) do geokey
+    id = GeoKeyID(geokey[1])
+    tag = geokey[2]
+    count = geokey[3]
+    value = geokey[4]
+    GeoKey(id, tag, count, value)
+  end
+
+  version = params[1]
+  revision = params[2]
+  minor = params[3]
+  nkeys = params[4]
+  GeoKeyDirectory(version, revision, minor, nkeys, geokeys)
+end
+
+function params(geokeydir::GeoKeyDirectory)
+  params = [geokeydir.version, geokeydir.revision, geokeydir.minor, geokeydir.nkeys]
+  for geokey in geokeydir.geokeys
+    append!(params, [UInt16(geokey.id), geokey.tag, geokey.count, geokey.value])
+  end
+  params
+end
+
+struct GeoDoubleParams
+  params::Vector{Float64}
+end
+
+params(geodouble::GeoDoubleParams) = geodouble.params
+
+struct GeoAsciiParams
+  params::String
+end
+
+params(geoascii::GeoAsciiParams) = geoascii.params
+
+struct ModelPixelScale
+  x::Float64
+  y::Float64
+  z::Float64
+end
+
+ModelPixelScale(; x=1.0, y=1.0, z=1.0) = ModelPixelScale(x, y, z)
+
+ModelPixelScale(params::Vector{Float64}) = ModelPixelScale(params[1], params[2], params[3])
+
+params(scale::ModelPixelScale) = [scale.x, scale.y, scale.z]
+
+struct ModelTiepoint
+  i::Float64
+  j::Float64
+  k::Float64
+  x::Float64
+  y::Float64
+  z::Float64
+end
+
+ModelTiepoint(; i=0.0, j=0.0, k=0.0, x=0.0, y=0.0, z=0.0) = ModelTiepoint(i, j, k, x, y, z)
+
+ModelTiepoint(params::Vector{Float64}) = ModelTiepoint(params[1], params[2], params[3], params[4], params[5], params[6])
+
+params(tiepoint::ModelTiepoint) = [tiepoint.i, tiepoint.j, tiepoint.k, tiepoint.x, tiepoint.y, tiepoint.z]
+
+struct ModelTransformation
+  A::Matrix{Float64}
+  b::Vector{Float64}
+end
+
+function ModelTransformation(; A=[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], b=[0.0, 0.0, 0.0])
+  sz = size(A)
+  if !allequal(sz)
+    throw(ArgumentError("`A` must be a square matrix"))
+  end
+  dim = first(sz)
+  if dim ∉ (2, 3)
+    throw(ArgumentError("only 2D and 3D transformations are supported"))
+  end
+  if dim ≠ length(b)
+    throw(ArgumentError("`A` and `b` must have the same dimension"))
+  end
+  A′, b′ = if dim == 2
+    [A[1, 1] A[1, 2] 0; A[2, 1] A[2, 2] 0; 0 0 0], [b[1], b[2], 0]
+  elseif dim == 3
+    A, b
+  else
+    throw(ArgumentError("only 2D and 3D transformations are supported"))
+  end
+  ModelTransformation(A′, b′)
+end
+
+function ModelTransformation(params::Vector{Float64})
+  A = [
+    params[1] params[2] params[3]
+    params[5] params[6] params[7]
+    params[9] params[10] params[11]
+  ]
+  b = [params[4], params[8], params[12]]
+  ModelTransformation(A, b)
+end
+
+function params(transform::ModelTransformation)
+  A = transform.A
+  b = transform.b
+  [A[1, 1], A[1, 2], A[1, 3], b[1], A[2, 1], A[2, 2], A[2, 3], b[2], A[3, 1], A[3, 2], A[3, 3], b[3], 0, 0, 0, 1]
+end
+
+struct Metadata
+  geokeydirectory::GeoKeyDirectory
+  geodoubleparams::Union{GeoDoubleParams,Nothing}
+  geoasciiparams::Union{GeoAsciiParams,Nothing}
+  modelpixelscale::Union{ModelPixelScale,Nothing}
+  modeltiepoint::Union{ModelTiepoint,Nothing}
+  modeltransformation::Union{ModelTransformation,Nothing}
+end
+
+Metadata(;
+  geokeydirectory=GeoKeyDirectory(),
+  geodoubleparams=nothing,
+  geoasciiparams=nothing,
+  modelpixelscale=nothing,
+  modeltiepoint=nothing,
+  modeltransformation=ModelTransformation()
+) = Metadata(geokeydirectory, geodoubleparams, geoasciiparams, modelpixelscale, modeltiepoint, modeltransformation)

src/save.jl

@@ -0,0 +1,52 @@
+# -----------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# -----------------------------------------------------------------
+
+const GeoTIFFType = Union{AbstractFloat,Signed,Unsigned}
+
+const WidePixelOrColorant = Union{WidePixel,Colorant}
+
+save(fname, geotiff::GeoTIFFImage) = TiffImages.save(fname, tiff(geotiff))
+
+function save(fname, tiff::AbstractTIFF; metadata=Metadata())
+  _setmetadata!(tiff, metadata)
+  TiffImages.save(fname, tiff)
+end
+
+save(fname, img::AbstractArray{<:WidePixelOrColorant}; kwargs...) = save(fname, DenseTaggedImage(img); kwargs...)
+
+function save(fname, channel₁::AbstractArray{T}, channelₙ::AbstractArray{T}...; kwargs...) where {T<:GeoTIFFType}
+  CT = _colordatatype(T)
+  colors = reinterpret(Gray{CT}, channel₁)
+  img = if length(channelₙ) > 0
+    extras = reinterpret.(CT, channelₙ)
+    mappedarray((color, extra...) -> WidePixel(color, extra), colors, extras...)
+  else
+    colors
+  end
+  save(fname, img; kwargs...)
+end
+
+# -----------------
+# HELPER FUNCTIONS
+# -----------------
+
+_colordatatype(::Type{T}) where {T<:AbstractFloat} = T
+_colordatatype(::Type{T}) where {T<:Unsigned} = Normed{T,sizeof(T) * 8}
+_colordatatype(::Type{T}) where {T<:Signed} = Fixed{T,sizeof(T) * 8 - 1}
+
+function _setmetadata!(tiff, metadata)
+  ifd = TiffImages.ifds(tiff)
+  _settag!(ifd, GeoKeyDirectoryTag, metadata.geokeydirectory)
+  _settag!(ifd, GeoDoubleParamsTag, metadata.geodoubleparams)
+  _settag!(ifd, GeoAsciiParamsTag, metadata.geoasciiparams)
+  _settag!(ifd, ModelPixelScaleTag, metadata.modelpixelscale)
+  _settag!(ifd, ModelTiepointTag, metadata.modeltiepoint)
+  _settag!(ifd, ModelTransformationTag, metadata.modeltransformation)
+end
+
+function _settag!(ifd, tag, geotag)
+  if !isnothing(geotag)
+    ifd[UInt16(tag)] = params(geotag)
+  end
+end