Skip to content

Commit a10a22c

Browse files
juliohmjuliohm
committed
Fix cartesianrange/slice
1 parent a526d38 commit a10a22c

File tree

4 files changed

+152
-345
lines changed

4 files changed

+152
-345
lines changed

‎src/indices.jl‎

Lines changed: 114 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,10 @@
22
# Licensed under the MIT License. See LICENSE in the project root.
33
# ------------------------------------------------------------------
44

5+
# ---------------
6+
# LINEAR INDICES
7+
# ---------------
8+
59
"""
610
indices(domain, geometry)
711
@@ -54,7 +58,7 @@ end
5458

5559
function indices(grid::CartesianGrid, box::Box)
5660
# cartesian range
57-
range = cartesianrange(grid, _boxlimits(box))
61+
range = cartesianrange(grid, box)
5862

5963
# convert to linear indices
6064
LinearIndices(size(grid))[range] |> vec
@@ -64,26 +68,125 @@ indices(grid::CartesianGrid, multi::Multi) = mapreduce(geom -> indices(grid, geo
6468

6569
function indices(grid::RectilinearGrid, box::Box)
6670
# cartesian range
67-
range = cartesianrange(grid, _boxlimits(box))
71+
range = cartesianrange(grid, box)
6872

6973
# convert to linear indices
7074
LinearIndices(size(grid))[range] |> vec
7175
end
7276

73-
# -----------------
74-
# HELPER FUNCTIONS
75-
# -----------------
77+
# ----------------
78+
# CARTESIAN RANGE
79+
# ----------------
80+
81+
"""
82+
cartesianrange(grid, box)
83+
84+
Return the Cartesian range of the elements of the `grid` that intersect with the `box`.
85+
"""
86+
cartesianrange(grid::Grid{M}, box::Box{M}) where {M} = _manifoldrange(M, grid, box)
87+
88+
_manifoldrange(::Type{<:𝔼}, grid::Grid, box::Box) = _euclideanrange(grid, box)
89+
90+
_manifoldrange(::Type{<:🌐}, grid::Grid, box::Box) = _geodesicrange(grid, box)
91+
92+
function _euclideanrange(grid::CartesianGrid, box::Box)
93+
# grid properties
94+
or = minimum(grid)
95+
sp = spacing(grid)
96+
sz = size(grid)
97+
98+
# intersection of boxes
99+
lo, up = extrema(boundingbox(grid) ∩ box)
100+
101+
# Cartesian indices of new corners
102+
ijkₛ = max.(ceil.(Int, (lo - or) ./ sp), 1)
103+
ijkₑ = min.(floor.(Int, (up - or) ./ sp) .+ 1, sz)
76104

77-
function _boxlimits(box::Box{𝔼{2}})
78-
min, max = convert.(Cartesian, coords.(extrema(box)))
79-
(min.x, max.x), (min.y, max.y)
105+
# Cartesian range from corner to corner
106+
CartesianIndex(Tuple(ijkₛ)):CartesianIndex(Tuple(ijkₑ))
80107
end
81108

82-
function _boxlimits(box::Box{𝔼{3}})
83-
min, max = convert.(Cartesian, coords.(extrema(box)))
84-
(min.x, max.x), (min.y, max.y), (min.z, max.z)
109+
function _euclideanrange(grid::RectilinearGrid, box::Box)
110+
# grid properties
111+
nd = paramdim(grid)
112+
113+
# intersection of boxes
114+
lo, up = to.(extrema(boundingbox(grid) ∩ box))
115+
116+
# integer coordinates of lower point
117+
ijkₛ = ntuple(nd) do i
118+
findlast(x -> x ≤ lo[i], xyz(grid)[i])
119+
end
120+
121+
# integer coordinates of upper point
122+
ijkₑ = ntuple(nd) do i
123+
findfirst(x -> x ≥ up[i], xyz(grid)[i])
124+
end
125+
126+
# integer coordinates of elements
127+
CartesianIndex(ijkₛ):CartesianIndex(ijkₑ .- 1)
85128
end
86129

130+
function _geodesicrange(grid::Grid, box::Box)
131+
nlon, nlat = vsize(grid)
132+
133+
boxmin = convert(LatLon, coords(minimum(box)))
134+
boxmax = convert(LatLon, coords(maximum(box)))
135+
136+
a = convert(LatLon, coords(vertex(grid, (1, 1))))
137+
b = convert(LatLon, coords(vertex(grid, (nlon, 1))))
138+
c = convert(LatLon, coords(vertex(grid, (1, nlat))))
139+
140+
swaplon = a.lon > b.lon
141+
swaplat = a.lat > c.lat
142+
143+
loninds = swaplon ? (nlon:-1:1) : (1:1:nlon)
144+
latinds = swaplat ? (nlat:-1:1) : (1:1:nlat)
145+
146+
gridlonₛ, gridlonₑ = swaplon ? (b.lon, a.lon) : (a.lon, b.lon)
147+
gridlatₛ, gridlatₑ = swaplat ? (c.lat, a.lat) : (a.lat, c.lat)
148+
149+
lonmin = max(boxmin.lon, gridlonₛ)
150+
latmin = max(boxmin.lat, gridlatₛ)
151+
lonmax = min(boxmax.lon, gridlonₑ)
152+
latmax = min(boxmax.lat, gridlatₑ)
153+
154+
iₛ = findlast(loninds) do i
155+
p = vertex(grid, (i, 1))
156+
c = convert(LatLon, coords(p))
157+
c.lon ≤ lonmin
158+
end
159+
iₑ = findfirst(loninds) do i
160+
p = vertex(grid, (i, 1))
161+
c = convert(LatLon, coords(p))
162+
c.lon ≥ lonmax
163+
end
164+
165+
jₛ = findlast(latinds) do i
166+
p = vertex(grid, (1, i))
167+
c = convert(LatLon, coords(p))
168+
c.lat ≤ latmin
169+
end
170+
jₑ = findfirst(latinds) do i
171+
p = vertex(grid, (1, i))
172+
c = convert(LatLon, coords(p))
173+
c.lat ≥ latmax
174+
end
175+
176+
if iₛ == iₑ || jₛ == jₑ
177+
throw(ArgumentError("the passed limits are not valid for the grid"))
178+
end
179+
180+
iₛ, iₑ = swaplon ? (iₑ, iₛ) : (iₛ, iₑ)
181+
jₛ, jₑ = swaplat ? (jₑ, jₛ) : (jₛ, jₑ)
182+
183+
CartesianIndex(loninds[iₛ], latinds[jₛ]):CartesianIndex(loninds[iₑ] - 1, latinds[jₑ] - 1)
184+
end
185+
186+
# -----------------
187+
# HELPER FUNCTIONS
188+
# -----------------
189+
87190
function _fill!(mask, grid, val, triangle)
88191
v = vertices(triangle)
89192

‎src/transforms/slice.jl‎

Lines changed: 38 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66
Slice(x=(xmin, xmax), y=(ymin, ymax), z=(zmin, zmax))
77
Slice(lat=(latmin, latmax), lon=(lonmin, lonmax))
88
9-
Retain the grid elements within `x` limits [`xmax`,`xmax`],
9+
Retain the domain elements within `x` limits [`xmax`,`xmax`],
1010
`y` limits [`ymax`,`ymax`] and `z` limits [`zmin`,`zmax`]
1111
in length units (default to meters), or within `lat` limits
1212
[`latmin`,`latmax`] and `lon` limits [`lonmin`,`lonmax`]
@@ -15,11 +15,10 @@ in degree units.
1515
## Examples
1616
1717
```julia
18-
Slice(x=(2, 4))
19-
Slice(x=(1u"km", 3u"km"))
20-
Slice(y=(1.2, 1.8), z=(2.4, 3.0))
21-
Slice(lat=(30, 60))
22-
Slice(lon=(45u"°", 90u"°"))
18+
Slice(x=(1000km, 3000km))
19+
Slice(x=(1000km, 2000km), y=(2000km, 5000km))
20+
Slice(lon=(0°, 90°))
21+
Slice(lon=(0°, 45°), lat=(0°, 45°))
2322
```
2423
"""
2524
struct Slice{T} <: GeometricTransform
@@ -30,54 +29,51 @@ Slice(; kwargs...) = Slice(values(kwargs))
3029

3130
parameters(t::Slice) = (; limits=t.limits)
3231

33-
preprocess(t::Slice, g::Grid) = cartesianrange(g, _fixlimits(boundingbox(g), t.limits))
32+
preprocess(t::Slice, d::Domain) = _sliceinds(d, _slicebox(boundingbox(d), t.limits))
3433

35-
apply(t::Slice, g::Grid) = g[preprocess(t, g)], nothing
34+
apply(t::Slice, d::Domain) = _slice(d, preprocess(t, d)), nothing
3635

3736
# -----------------
3837
# HELPER FUNCTIONS
3938
# -----------------
4039

41-
function _fixlimits(box::Box{<:𝔼}, limits)
42-
lims = _xyzlimits(limits)
43-
min = convert(Cartesian, coords(minimum(box)))
44-
max = convert(Cartesian, coords(maximum(box)))
45-
_minmax(min, max, lims)
46-
end
47-
48-
function _fixlimits(box::Box{🌐}, limits)
49-
lims = _latlonlimits(limits)
50-
min = convert(LatLon, coords(minimum(box)))
51-
max = convert(LatLon, coords(maximum(box)))
52-
_minmax(min, max, lims)
53-
end
40+
_slice(d::Domain, inds) = view(d, inds)
41+
_slice(g::Grid, inds::CartesianIndices) = getindex(g, inds)
5442

55-
_xyzlimits(limits) = (
56-
x=haskey(limits, :x) ? _aslen.(limits.x) : nothing,
57-
y=haskey(limits, :y) ? _aslen.(limits.y) : nothing,
58-
z=haskey(limits, :z) ? _aslen.(limits.z) : nothing
59-
)
43+
_sliceinds(d::Domain, b) = indices(d, b)
44+
_sliceinds(g::CartesianGrid, b) = cartesianrange(g, b)
45+
_sliceinds(g::RectilinearGrid, b) = cartesianrange(g, b)
46+
_sliceinds(g::Grid{🌐}, b::Box{🌐}) = cartesianrange(g, b)
6047

61-
_latlonlimits(limits) =
62-
(lat=haskey(limits, :lat) ? _asdeg.(limits.lat) : nothing, lon=haskey(limits, :lon) ? _asdeg.(limits.lon) : nothing)
63-
64-
function _minmax(min::Cartesian2D, max::Cartesian2D, lims)
65-
xmin, xmax = isnothing(lims.x) ? (min.x, max.x) : lims.x
66-
ymin, ymax = isnothing(lims.y) ? (min.y, max.y) : lims.y
67-
(xmin, xmax), (ymin, ymax)
48+
function _slicebox(box::Box{𝔼{2}}, limits)
49+
min = convert(Cartesian, coords(minimum(box)))
50+
max = convert(Cartesian, coords(maximum(box)))
51+
xmin, xmax = get(limits, :x, (min.x, max.x))
52+
ymin, ymax = get(limits, :y, (min.y, max.y))
53+
bmin = _aslen.((xmin, ymin))
54+
bmax = _aslen.((xmax, ymax))
55+
Box(withcrs(box, bmin), withcrs(box, bmax))
6856
end
6957

70-
function _minmax(min::Cartesian3D, max::Cartesian3D, lims)
71-
xmin, xmax = isnothing(lims.x) ? (min.x, max.x) : lims.x
72-
ymin, ymax = isnothing(lims.y) ? (min.y, max.y) : lims.y
73-
zmin, zmax = isnothing(lims.z) ? (min.z, max.z) : lims.z
74-
(xmin, xmax), (ymin, ymax), (zmin, zmax)
58+
function _slicebox(box::Box{𝔼{3}}, limits)
59+
min = convert(Cartesian, coords(minimum(box)))
60+
max = convert(Cartesian, coords(maximum(box)))
61+
xmin, xmax = get(limits, :x, (min.x, max.x))
62+
ymin, ymax = get(limits, :y, (min.y, max.y))
63+
zmin, zmax = get(limits, :z, (min.z, max.z))
64+
bmin = _aslen.((xmin, ymin, zmin))
65+
bmax = _aslen.((xmax, ymax, zmax))
66+
Box(withcrs(box, bmin), withcrs(box, bmax))
7567
end
7668

77-
function _minmax(min::LatLon, max::LatLon, lims)
78-
lonmin, lonmax = isnothing(lims.lon) ? (min.lon, max.lon) : lims.lon
79-
latmin, latmax = isnothing(lims.lat) ? (min.lat, max.lat) : lims.lat
80-
(lonmin, lonmax), (latmin, latmax)
69+
function _slicebox(box::Box{🌐}, limits)
70+
min = convert(LatLon, coords(minimum(box)))
71+
max = convert(LatLon, coords(maximum(box)))
72+
latmin, latmax = get(limits, :lat, (min.lat, max.lat))
73+
lonmin, lonmax = get(limits, :lon, (min.lon, max.lon))
74+
bmin = _asdeg.((latmin, lonmin))
75+
bmax = _asdeg.((latmax, lonmax))
76+
Box(withcrs(box, bmin, LatLon), withcrs(box, bmax, LatLon))
8177
end
8278

8379
_aslen(x::Len) = float(x)

0 commit comments

Comments
 (0)