-
Notifications
You must be signed in to change notification settings - Fork 33
/
reproject.jl
81 lines (69 loc) · 3.24 KB
/
reproject.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
"""
reproject(obj; crs)
Reproject the lookups of `obj` to a different crs.
This is a lossless operation for the raster data, as only the
lookup values change. This is only possible when the axes of source
and destination projections are alligned: the change is usually from
a [`Regular`](@ref) and an [`Irregular`](@ref) lookup spans.
For converting between projections that are rotated,
skewed or warped in any way, use [`resample`](@ref).
Dimensions without an `AbstractProjected` lookup (such as a `Ti` dimension)
are silently returned without modification.
# Arguments
- `obj`: a `Lookup`, `Dimension`, `Tuple` of `Dimension`, `Raster` or `RasterStack`.
$CRS_KEYWORD
"""
reproject(x; crs::GeoFormat) = reproject(crs, x)
reproject(target::GeoFormat, x) = rebuild(x; dims=reproject(target, dims(x)))
reproject(target::GeoFormat, dims::Tuple) = map(d -> reproject(target, d), dims)
reproject(target::GeoFormat, l::Lookup) = l
reproject(target::GeoFormat, dim::Dimension) = rebuild(dim, reproject(target, lookup(dim)))
function reproject(target::GeoFormat, l::AbstractProjected)
source = crs(l)
newdata = reproject(source, target, l.dim, parent(l))
newlookup = rebuild(l; data=newdata, crs=target)
if _checkregular(newdata)
return set(newlookup, Regular(stepof(newdata)))
else
newbounds = reproject(crs(l), target, l.dim, bounds(l))
return set(newlookup, Irregular(newbounds))
end
end
"""
reproject(source::GeoFormat, target::GeoFormat, dim::Dimension, val)
`reproject` uses ArchGDAL.reproject, but implemented for a reprojecting
a value array of values, a single dimension at a time.
"""
function reproject(source, target, dim, val)
if source == target
return val
else
return _reproject(source, target, dim, val)
end
end
_reproject(source::GeoFormat, target::GeoFormat, dim, val) = val
_reproject(source::Nothing, target::GeoFormat, dim, val) = val
_reproject(source::GeoFormat, target::Nothing, dim, val) = val
_reproject(source::Nothing, target::Nothing, dim, val) = val
function _reproject(source::GeoFormat, target::GeoFormat, dim::Union{XDim,YDim}, val::Number)
# This is a dumb way to do this. But it save having to inspect crs,
# and prevents reprojections that dont make sense from working.
# A better method for this should be implemented in future.
return first(reproject(source, target, dim, [val]))
end
function _reproject(source::GeoFormat, target::GeoFormat, dim, vals::NTuple{N}) where N
reps = reproject(source, target, dim, [vals...])
return ntuple(x -> reps[x], N)
end
function _reproject(source::GeoFormat, target::GeoFormat, dim::Union{XDim,YDim}, vals::AbstractArray)
reshape(reproject(source, target, dim, vec(vals)), size(vals))
end
function _reproject(source::Nothing, target::GeoFormat, dim::Union{XDim,YDim}, vals::AbstractArray)
reshape(reproject(source, target, dim, vec(vals)), size(vals))
end
function _reproject(source::GeoFormat, target::GeoFormat, dim::Union{XDim,YDim}, vals::AbstractVector)
error("Rasters.jl requires backends to be loaded externally as of version 0.8. Run `using ArchGDAL` to use `reproject`")
end
# Guess the step for arrays
stepof(A::AbstractArray) = (last(A) - first(A)) / (length(A) - 1)
stepof(A::AbstractRange) = step(A)