-
Notifications
You must be signed in to change notification settings - Fork 33
/
mosaic.jl
240 lines (205 loc) · 7.9 KB
/
mosaic.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
"""
mosaic(f, regions...; missingval, atol)
mosaic(f, regions; missingval, atol)
Combine `regions` into a single raster.
# Arguments
- `f`: A reducing function (`mean`, `sum`, `first`, `last` etc.)
for values where `regions` overlap.
- `regions`: Iterable or splatted `Raster` or `RasterStack`.
# Keywords
- `missingval`: Fills empty areas, and defualts to the
`missingval` of the first region.
- `atol`: Absolute tolerance for comparison between index values.
This is often required due to minor differences in range values
due to floating point error. It is not applied to non-float dimensions.
A tuple of tolerances may be passed, matching the dimension order.
$FILENAME_KEYWORD
$SUFFIX_KEYWORD
If your mosaic has has apparent line errors, increase the `atol` value.
# Example
Here we cut out Australia and Africa from a stack, and join them with `mosaic`.
```jldoctest
using Rasters, RasterDataSources, ArchGDAL, Plots
st = RasterStack(WorldClim{Climate}; month=1);
africa = st[X(-20.0 .. 60.0), Y(-40.0 .. 35.0)]
a = plot(africa)
aus = st[X(100.0 .. 160.0), Y(-50.0 .. -10.0)]
b = plot(aus)
# Combine with mosaic
mos = mosaic(first, aus, africa)
c = plot(mos)
savefig(a, "docs/build/mosaic_example_africa.png")
savefig(b, "docs/build/mosaic_example_aus.png")
savefig(c, "docs/build/mosaic_example_combined.png")
nothing
# output
```
### Individual continents
![arica](../build/mosaic_example_africa.png)
![aus](../build/mosaic_example_aus.png)
### Mosaic of continents
![mosaic](../build/mosaic_example_combined.png)
$EXPERIMENTAL
"""
function mosaic(f::Function, r1::RasterStackOrArray, rs::RasterStackOrArray...; kw...)
mosaic(f, (r1, rs...); kw...)
end
mosaic(f::Function, regions; kw...) = _mosaic(f, first(regions), regions; kw...)
function _mosaic(f::Function, ::AbstractRaster, regions;
missingval=missingval(first(regions)), filename=nothing, suffix=nothing, kw...
)
missingval = missingval isa Nothing ? missing : missingval
T = Base.promote_type(typeof(missingval), Base.promote_eltype(regions...))
dims = _mosaic(Tuple(map(DD.dims, regions)))
l1 = first(regions)
A = create(filename, T, dims; name=name(l1), missingval, metadata=metadata(l1))
open(A; write=true) do a
mosaic!(f, a, regions; missingval, kw...)
end
return A
end
function _mosaic(f::Function, ::AbstractRasterStack, regions;
filename=nothing, suffix=keys(first(regions)), kw...
)
layers = map(suffix, map(values, regions)...) do s, A...
mosaic(f, A...; filename, suffix=s, kw...)
end
return DD.rebuild_from_arrays(first(regions), Tuple(layers))
end
"""
mosaic!(f, x, regions...; missingval, atol)
mosaic!(f, x, regions::Tuple; missingval, atol)
Combine `regions` in `x` using the function `f`.
# Arguments
- `f` a function (e.g. `mean`, `sum`, `first` or `last`) that is applied to
values where `regions` overlap.
- `x`: A `Raster` or `RasterStack`. May be a an opened disk-based `Raster`,
the result will be written to disk.
With the current algorithm, the read speed is slow.
- `regions`: source objects to be joined. These should be memory-backed
(use `read` first), or may experience poor performance. If all objects have
the same extent, `mosaic` is simply a merge.
# Keywords
- `missingval`: Fills empty areas, and defualts to the `missingval/
of the first layer.
- `atol`: Absolute tolerance for comparison between index values.
This is often required due to minor differences in range values
due to floating point error. It is not applied to non-float dimensions.
A tuple of tolerances may be passed, matching the dimension order.
# Example
Cut out Australia and Africa stacks, then combined them
into a single stack.
```jldoctest
using Rasters, RasterDataSources, ArchGDAL, Statistics, Plots
st = read(RasterStack(WorldClim{Climate}; month=1))
aus = st[X=100.0 .. 160.0, Y=-50.0 .. -10.0]
africa = st[X=-20.0 .. 60.0, Y=-40.0 .. 35.0]
mosaic!(first, st, aus, africa)
plot(st)
savefig("docs/build/mosaic_bang_example.png")
nothing
# output
```
![mosaic](../build/mosaic_bang_example.png)
$EXPERIMENTAL
"""
mosaic!(f::Function, x::RasterStackOrArray, regions::RasterStackOrArray...; kw...) = mosaic!(f, x, regions; kw...)
function mosaic!(f::Function, A::AbstractRaster{T}, regions;
missingval=missingval(A), atol=_default_atol(T)
) where T
_without_mapped_crs(A) do A1
broadcast!(A1, DimKeys(A1; atol)) do ds
# Get all the regions that have this point
ls = foldl(regions; init=()) do acc, l
if DD.hasselection(l, ds)
v = l[ds...]
(acc..., l)
else
acc
end
end
values = foldl(ls; init=()) do acc, l
v = l[ds...]
if isnothing(Rasters.missingval(l))
(acc..., v)
elseif ismissing(Rasters.missingval(l))
ismissing(v) ? acc : (acc..., v)
else
v === Rasters.missingval(l) ? acc : (acc..., v)
end
end
if length(values) === 0
missingval
else
f(values)
end
end
end
return A
end
function mosaic!(f::Function, st::AbstractRasterStack, regions::Tuple; kw...)
map(st, regions...) do A, r...
mosaic!(f, A, r; kw...)
end
end
_mosaic(alldims::Tuple{<:DimTuple,Vararg{<:DimTuple}}) = map(_mosaic, alldims...)
function _mosaic(dims::Dimension...)
map(dims) do d
DD.comparedims(first(dims), d; val=false, length=false, valtype=true)
end
return rebuild(first(dims), _mosaic(lookup(dims)))
end
_mosaic(lookups::LookupTuple) = _mosaic(first(lookups), lookups)
function _mosaic(lookup::Categorical, lookups::LookupTuple)
newindex = union(lookups...)
if order isa ForwardOrdered
newindex = sort(newindex; order=LA.ordering(order(lookup)))
end
return rebuild(lookup; data=newindex)
end
function _mosaic(lookup::AbstractSampled, lookups::LookupTuple)
order(lookup) isa Unordered && throw(ArgumentError("Cant mozaic an Unordered lookup"))
return _mosaic(span(lookup), lookup, lookups)
end
function _mosaic(span::Regular, lookup::AbstractSampled, lookups::LookupTuple)
newindex = if order(lookup) isa ForwardOrdered
mi = minimum(map(first, lookups))
ma = maximum(map(last, lookups))
if mi isa AbstractFloat
# Handle slight range erorrs to make sure
# we dont drop one step of the range
mi:step(span):ma + 2eps(ma)
else
mi:step(span):ma
end
else
mi = minimum(map(last, lookups))
ma = maximum(map(first, lookups))
if mi isa AbstractFloat
ma:step(span):mi - 2eps(mi)
else
ma:step(span):mi
end
end
return rebuild(lookup; data=newindex)
end
function _mosaic(::Irregular, lookup::AbstractSampled, lookups::LookupTuple)
newindex = sort(union(map(parent, lookups)...); order=LA.ordering(order(lookup)))
return rebuild(lookup; data=newindex)
end
function _mosaic(span::Explicit, lookup::AbstractSampled, lookups::LookupTuple)
# TODO make this less fragile to floating point innaccuracy
newindex = sort(union(map(parent, lookups)...); order=LA.ordering(order(lookup)))
bounds = map(val ∘ DD.span, lookups)
lower = map(b -> view(b, 1, :), bounds)
upper = map(b -> view(b, 2, :), bounds)
newlower = sort(union(lower...); order=LA.ordering(order(lookup)))
newupper = sort(union(upper...); order=LA.ordering(order(lookup)))
newbounds = vcat(permutedims(newlower), permutedims(newupper))
return rebuild(lookup; data=newindex, span=Explicit(newbounds))
end
# These are pretty random default, but seem to work
_default_atol(T::Type{<:Float32}) = 100eps(T)
_default_atol(T::Type{<:Float64}) = 1000eps(T)
_default_atol(T::Type{<:Integer}) = T(1)
_default_atol(::Type) = nothing