Skip to content

Commit

Permalink
Allow for GeoArray to be a Matrix (#162)
Browse files Browse the repository at this point in the history 
* Allow for GeoArray to be a Matrix. Make use of Extents. Change RoundingMode.

* Commit remaining files.

* Disallow recursive GeoArrays.

* Only check first two dimension for similar.

* Update docs.
  • Loading branch information
@evetion
evetion authored Jun 11, 2024
1 parent 6c9cfba commit d53cbae
Show file tree
Hide file tree
Showing 12 changed files with 188 additions and 98 deletions.
26 changes: 26 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

## [Unreleased]

## [0.9.0] - 2024-06-11

### Changes
- [Breaking] Changed type signature of GeoArray to include the number of dimensions. This allows
for single "band" GeoArrays to be represented as matrices as well. This should make it easier to
work with single band rasters and the image ecosystem.
- [Breaking] `bbox` and friends now return an `Extents.Extent` instead of a `NamedTuple`.
- [Breaking] Reverted rename of `equals` to `Base.isequal`, now called `isgeoequal`.
- [Breaking] `getindex`, `indices` and `sample` now use the rounding mode `RoundNearestTiesUp` instead of `RoundNearest`.

### Deprecated
- [Breaking] Bounding box input for `crop`, `warp` and others is now deprecated. Use an `Extent` instead.

### Fixed
- Try to release file lock on close as soon as possible.
- Indexing a GeoArray with `[:, :]` now returns a GeoArray

### Added
- Added `enable_online_warp` function to enable online access to PROJ data for `warp`.
- Added `rounding` argument to `indices` to control rounding mode used.

## [0.8.5] - 2024-01-07
- Fix small bug in metadata
- Move GeoStatsBase into an extension
Expand Down Expand Up @@ -70,6 +91,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Fixed iterate specification so `sum` on a GeoArray is correct

[unreleased]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.1...HEAD
[0.9.0]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.5...v0.9.0
[0.8.5]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.4...v0.8.5
[0.8.4]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.3...v0.8.4
[0.8.3]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.2...v0.8.3
[0.8.2]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.1...v0.8.2
[0.8.1]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.0...v0.8.1
[0.8.0]: https://github.com/evetion/GeoArrays.jl/compare/v0.7.13...v0.8.0
[0.7.13]: https://github.com/evetion/GeoArrays.jl/compare/v0.7.12...v0.7.13
Expand Down
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "GeoArrays"
uuid = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab"
authors = ["Maarten Pronk <[email protected]>"]
version = "0.8.5"
version = "0.9.0"

[weakdeps]
GeoStatsBase = "323cb8eb-fbf6-51c0-afd0-f8fba70507b2"
Expand Down
9 changes: 5 additions & 4 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ This packages takes its inspiration from Python's [rasterio](https://github.com/
## Installation

```julia
(v1.8) pkg> add GeoArrays
(v1.10) pkg> add GeoArrays
```

## Examples
Expand All @@ -30,7 +30,7 @@ Read a GeoTIFF file and display its information, i.e. AffineMap and projection (
# Read TIF file
julia> fn = download("https://github.com/yeesian/ArchGDALDatasets/blob/master/data/utmsmall.tif?raw=true")
julia> geoarray = GeoArrays.read(fn)
100x100x1 Array{UInt8,3} with AffineMap([60.0 0.0; 0.0 -60.0], [440720.0, 3.75132e6]) and CRS PROJCS["NAD27 / UTM zone 11N"...
100x100 Matrix{UInt8} with AffineMap([60.0 0.0; 0.0 -60.0], [440720.0, 3.75132e6]) and CRS PROJCS["NAD27 / UTM zone 11N"...

# Affinemap containing offset and scaling
julia> geoarray.f
Expand Down Expand Up @@ -152,10 +152,11 @@ For example, we can vertically transform from the ellipsoid
to the EGM2008 geoid using EPSG code 3855.
Note that the underlying PROJ library needs to find the geoidgrids,
so if they're not available locally, one needs to set `ENV["PROJ_NETWORK"] = "ON"`
as early as possible, ideally before loading GeoArrays.
as early as possible, ideally before loading GeoArrays, or use `enable_online_warp`.
```julia
enable_online_warp() # If you don't have PROJ grids locally
ga = GeoArray(zeros((360, 180)))
bbox!(ga, (min_x=-180, min_y=-90, max_x=180, max_y=90))
bbox!(ga, Extent(X(-180, 180), Y(-90, 90)))
crs!(ga, GeoFormatTypes.EPSG(4979)) # WGS83 in 3D (reference to ellipsoid)
ga2 = GeoArrays.warp(ga, Dict("t_srs" => "EPSG:4326+3855"))
```
Expand Down
2 changes: 1 addition & 1 deletion src/GeoArrays.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ using StaticArrays
import GeoInterface as GI
using IterTools: partition
import DataAPI
import Extents
using Extents: Extent, intersects
using PrecompileTools # this is a small dependency

include("geoarray.jl")
Expand Down
97 changes: 56 additions & 41 deletions src/geoarray.jl
View file
Original file line number Diff line number Diff line change
@@ -1,21 +1,22 @@
const NumberOrMissing = Union{Number,Union{Missing,Number}}
const MatrixorArray = Union{<:AbstractArray{T,2},<:AbstractArray{T,3}} where {T}

"""
GeoArray{T::NumberOrMissing,A<:AbstractArray{T,3}} <: AbstractArray{T,3}
GeoArray{T::NumberOrMissing,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
A GeoArray is an AbstractArray, an AffineMap for calculating coordinates based on the axes and a CRS definition to interpret these coordinates into in the real world.
It's three dimensional and can be seen as a stack (3D) of 2D geospatial rasters (bands), the dimensions are :x, :y, and :bands.
The AffineMap and CRS (coordinates) only operate on the :x and :y dimensions.
"""
mutable struct GeoArray{T<:NumberOrMissing,A<:AbstractArray{T,3}} <: AbstractArray{T,3}
mutable struct GeoArray{T<:NumberOrMissing,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
A::A
f::CoordinateTransformations.AffineMap{StaticArrays.SMatrix{2,2,Float64,4},StaticArrays.SVector{2,Float64}}
crs::GFT.WellKnownText{GFT.CRS}
metadata::Dict{String}
end

"""
GeoArray(A::AbstractArray{T,3} where T <: NumberOrMissing)
GeoArray(A::AbstractArray{T,2|3} where T <: NumberOrMissing)
Construct a GeoArray from any Array. A default `AffineMap` and `CRS` will be generated.
Expand All @@ -25,62 +26,62 @@ julia> GeoArray(rand(10,10,1))
10x10x1 Array{Float64, 3} with AffineMap([1.0 0.0; 0.0 1.0], [0.0, 0.0]) and undefined CRS
```
"""
GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}) = GeoArray(A, geotransform_to_affine(SVector(0.0, 1.0, 0.0, 0.0, 0.0, 1.0)), "", Dict{String,Any}())
GeoArray(A::MatrixorArray where {T<:NumberOrMissing}) = GeoArray(A, geotransform_to_affine(SVector(0.0, 1.0, 0.0, 0.0, 0.0, 1.0)), "", Dict{String,Any}())
"""
GeoArray(A::AbstractArray{T,3} where T <: NumberOrMissing, f::AffineMap)
Construct a GeoArray from any Array and an `AffineMap` that specifies the coordinates. A default `CRS` will be generated.
"""
GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), ""), Dict{String,Any}())
GeoArray(A::MatrixorArray where {T<:NumberOrMissing}, f::AffineMap) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), ""), Dict{String,Any}())

"""
GeoArray(A::AbstractArray{T,3} where T <: NumberOrMissing, f::AffineMap, crs::String)
GeoArray(A::AbstractArray{T,2|3} where T <: NumberOrMissing, f::AffineMap, crs::String)
Construct a GeoArray from any Array and an `AffineMap` that specifies the coordinates and `crs` string in WKT format.
"""
GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap, crs::String) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs), Dict{String,Any}())
GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap, crs::String, d::Dict{String}) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs), d)
GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap, crs::GFT.WellKnownText{GFT.CRS}) = GeoArray(A, f, crs, Dict{String,Any}())

GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap{Matrix{Float64},Vector{Float64}}, crs::GFT.WellKnownText{GFT.CRS}) = GeoArray(A, AffineMap(SMatrix{2,2}(f.linear), SVector{2}(f.translation)), crs, Dict{String,Any}())
GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap{Matrix{Float64},Vector{Float64}}, crs::GFT.WellKnownText{GFT.CRS}, d::Dict{String}) = GeoArray(A, AffineMap(SMatrix{2,2}(f.linear), SVector{2}(f.translation)), crs, d)

"""
GeoArray(A::AbstractArray{T,2} where T <: NumberOrMissing)
Construct a GeoArray from any Matrix, a third singleton dimension will be added automatically.
# Examples
```julia-repl
julia> GeoArray(rand(10,10))
10x10x1 Array{Float64, 3} with AffineMap([1.0 0.0; 0.0 1.0], [0.0, 0.0]) and undefined CRS
```
"""
GeoArray(A::AbstractArray{T,2} where {T<:NumberOrMissing}, args...) = GeoArray(reshape(A, size(A)..., 1), args...)
GeoArray(A) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs), Dict{String,Any}())
GeoArray(A, f, crs::String, d) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs), d)
GeoArray(A, f, crs) = GeoArray(A, f, crs, Dict{String,Any}())
GeoArray(A, f::AffineMap{Matrix{Float64},Vector{Float64}}, args...) = GeoArray(A, AffineMap(SMatrix{2,2}(f.linear), SVector{2}(f.translation)), args...)
GeoArray(ga::T, args...) where {T<:GeoArray} = GeoArray(parent(ga), args...)

"""
GeoArray(A::AbstractArray{T,3} where T <: NumberOrMissing, x::AbstractRange, y::AbstractRange, args...)
GeoArray(A::AbstractArray{T,2|3} where T <: NumberOrMissing, x::AbstractRange, y::AbstractRange, args...)
Construct a GeoArray any Array and it's coordinates from `AbstractRange`s for each dimension.
"""
function GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, x::AbstractRange, y::AbstractRange, args...)
function GeoArray(A::MatrixorArray where {T<:NumberOrMissing}, x::AbstractRange, y::AbstractRange, args...)
size(A)[1:2] != (length(x), length(y)) && error("Size of `GeoArray` $(size(A)) does not match size of (x,y): $((length(x), length(y))). Note that this function takes *center coordinates*.")
f = unitrange_to_affine(x, y)
GeoArray(A, f, args...)
end

# Behave like an Array
Base.size(ga::GeoArray) = size(ga.A)
_size(ga::GeoArray{T,2}) where {T} = (size(ga.A)..., 1)
_size(ga::GeoArray{T,3}) where {T} = size(ga.A)
Base.IndexStyle(::Type{<:GeoArray}) = IndexCartesian()
Base.similar(ga::GeoArray, t::Type) = GeoArray(similar(ga.A, t), ga.f, ga.crs, ga.metadata)
Base.similar(ga::GeoArray, t::Type) = GeoArray(similar(parent(ga), t), ga.f, ga.crs, ga.metadata)
function Base.similar(ga::GeoArray, A::MatrixorArray)
size(ga.A)[1:2] == size(A)[1:2] || error(lazy"Size of `GeoArray` $(size(ga.A)) does not match size of `A`: $(size(A)).")
GeoArray(A, ga.f, ga.crs, ga.metadata)
end
Base.similar(ga::GeoArray, A::GeoArray) = similar(ga, parent(A))
Base.iterate(ga::GeoArray) = iterate(ga.A)
Base.iterate(ga::GeoArray, state) = iterate(ga.A, state)
Base.length(ga::GeoArray) = length(ga.A)
Base.parent(ga::GeoArray) = ga.A
Base.eltype(::Type{GeoArray{T}}) where {T} = T
Base.eltype(::Type{GeoArray{T,N}}) where {T,N} = T
Base.show(io::IO, ::MIME"text/plain", ga::GeoArray) = show(io, ga)

Base.convert(::Type{GeoArray{T}}, ga::GeoArray) where {T} = GeoArray(convert(Array{T}, ga.A), ga.f, ga.crs, ga.metadata)
function Base.convert(::Type{GeoArray{T}}, ga::GeoArray{X,N}) where {T,N,X}
A = convert(Array{T,N}, ga.A)
GeoArray(A, ga.f, ga.crs, ga.metadata)
end
function Base.convert(::Type{GeoArray{T,N}}, ga::GeoArray) where {T,N}
A = convert(Array{T,N}, ga.A)
GeoArray(A, ga.f, ga.crs, ga.metadata)
end

Base.BroadcastStyle(::Type{<:GeoArray}) = Broadcast.ArrayStyle{GeoArray}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GeoArray}}, ::Type{ElType}) where {ElType}
Expand Down Expand Up @@ -114,17 +115,26 @@ julia> ga[2:3,2:3,1]
2x2x1 Array{Float64, 3} with AffineMap([1.0 0.0; 0.0 1.0], [1.0, 1.0]) and undefined CRS
```
"""
function Base.getindex(ga::GeoArray, i::AbstractRange, j::AbstractRange, k::Union{Colon,AbstractRange,Integer})
A = getindex(ga.A, i, j, k)
function Base.getindex(ga::GeoArray{T,N}, i::AbstractRange, j::AbstractRange, k::Union{Colon,AbstractRange,Integer}=Colon()) where {T,N}
A = N == 2 ? getindex(ga.A, i, j) : getindex(ga.A, i, j, k)
x, y = first(i) - 1, first(j) - 1
t = ga.f(SVector(x, y))
l = ga.f.linear * SMatrix{2,2}([step(i) 0; 0 step(j)])
GeoArray(A, AffineMap(l, t), crs(ga), ga.metadata)
end
Base.getindex(ga::GeoArray, i::AbstractRange, j::AbstractRange) = Base.getindex(ga, i, j, :)
function Base.getindex(ga::GeoArray{T,N}, i::Colon, j::Colon, k::Union{Colon,AbstractRange,Integer}=Colon()) where {T,N}
A = N == 2 ? getindex(ga.A, i, j) : getindex(ga.A, i, j, k)
GeoArray(A, ga.f, ga.crs, ga.metadata)
end

Base.getindex(ga::GeoArray{T,3}, i::AbstractRange, j::AbstractRange) where {T} = Base.getindex(ga, i, j, :)

Base.getindex(ga::GeoArray, I::Vararg{Integer,2}) = getindex(ga.A, I...)
Base.getindex(ga::GeoArray, I::Vararg{Integer,3}) = getindex(ga.A, I...)
function Base.getindex(ga::GeoArray{T,2}, I::Vararg{Integer,3}) where {T}
I[3] != 1 && throw(BoundsError(ga, I))
getindex(ga.A, I[1], I[2])
end

# Getindex and setindex! with floats
"""
Expand All @@ -146,10 +156,14 @@ function Base.getindex(ga::GeoArray, I::SVector{2,<:AbstractFloat})
end
Base.getindex(ga::GeoArray, I::Vararg{AbstractFloat,2}) = getindex(ga, SVector{2}(I))

function Base.setindex!(ga::GeoArray, v, I::SVector{2,AbstractFloat})
i = indices(ga, I, Center())
function Base.setindex!(ga::GeoArray{T,3}, v, I::SVector{2,AbstractFloat}) where {T}
i = indices(ga, I, Center(), RoundNearestTiesUp)
ga.A[i, :] .= v
end
function Base.setindex!(ga::GeoArray{T,2}, v, I::SVector{2,AbstractFloat}) where {T}
i = indices(ga, I, Center(), RoundNearestTiesUp)
ga.A[i] = v
end
Base.setindex!(ga::GeoArray, v, I::Vararg{AbstractFloat,2}) = setindex!(ga, v, SVector{2}(I))
Base.setindex!(ga::GeoArray, v, I::Vararg{Union{<:Integer,<:AbstractRange{<:Integer}},2}) = setindex!(ga.A, v, I...)
Base.setindex!(ga::GeoArray, v, I::Vararg{Union{<:Integer,<:AbstractRange{<:Integer}},3}) = setindex!(ga.A, v, I...)
Expand Down Expand Up @@ -193,22 +207,23 @@ coords(ga::GeoArray, p::Tuple{<:Integer,<:Integer}, strategy::AbstractStrategy=C
coords(ga::GeoArray, p::CartesianIndex{2}, strategy::AbstractStrategy=Center()) = coords(ga, SVector{2}(p.I), strategy)

"""
indices(ga::GeoArray, p::SVector{2,<:Real}, strategy::AbstractStrategy)
indices(ga::GeoArray, p::SVector{2,<:Real}, strategy::AbstractStrategy, rounding::RoundingMode)
Retrieve logical indices of the cell represented by coordinates `p`.
`strategy` can be used to define whether the coordinates represent the center (`Center`) or the top left corner (`Vertex`) of the cell.
`rounding` can be used to define how the coordinates are rounded to the nearest integer index.
See `coords` for the inverse function.
"""
function indices(ga::GeoArray, p::SVector{2,<:Real}, strategy::AbstractStrategy)
CartesianIndex(Tuple(round.(Int, inv(ga.f)(p) .+ strategy.offset)))
function indices(ga::GeoArray, p::SVector{2,<:Real}, strategy::AbstractStrategy, rounding::RoundingMode)
CartesianIndex(Tuple(round.(Int, inv(ga.f)(p) .+ strategy.offset, rounding)))
end
indices(ga::GeoArray, p::AbstractVector{<:Real}, strategy::AbstractStrategy=Center()) = indices(ga, SVector{2}(p), strategy)
indices(ga::GeoArray, p::Tuple{<:Real,<:Real}, strategy::AbstractStrategy=Center()) = indices(ga, SVector{2}(p), strategy)
indices(ga::GeoArray, p::AbstractVector{<:Real}, strategy::AbstractStrategy=Center(), rounding::RoundingMode=RoundNearestTiesUp) = indices(ga, SVector{2}(p), strategy, rounding)
indices(ga::GeoArray, p::Tuple{<:Real,<:Real}, strategy::AbstractStrategy=Center(), rounding::RoundingMode=RoundNearestTiesUp) = indices(ga, SVector{2}(p), strategy, rounding)


# Generate coordinates for complete GeoArray
function coords(ga::GeoArray, strategy::AbstractStrategy=Center())
(ui, uj) = size(ga)[1:2]
(ui, uj) = size(ga)
extra = typeof(strategy) == Center ? 0 : 1
(coords(ga, SVector{2}(i, j), strategy) for i in 1:ui+extra, j in 1:uj+extra)
end
Expand Down
6 changes: 1 addition & 5 deletions src/geointerface.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,2 @@
function GI.extent(ga::GeoArray)
bb = bbox(ga)
Extents.Extent(X=(bb.min_x, bb.max_x), Y=(bb.min_y, bb.max_y))
end

GI.extent(ga::GeoArray) = bbox(ga)
GI.crs(ga::GeoArray) = isempty(GFT.val(crs(ga))) ? nothing : crs(ga)
Loading

2 comments on commit d53cbae

@evetion
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/108709

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.9.0 -m "<description of version>" d53cbae653f168f9bca1ff7723ac3198b2878f1f
git push origin v0.9.0

Please sign in to comment.