-
Notifications
You must be signed in to change notification settings - Fork 282
/
geometry.py
229 lines (186 loc) · 8.24 KB
/
geometry.py
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
# Copyright Iris contributors
#
# This file is part of Iris and is released under the BSD license.
# See LICENSE in the root of the repository for full licensing details.
"""Various utilities related to geometric operations.
.. note::
This module requires :mod:`shapely`.
"""
import warnings
import numpy as np
from shapely.geometry import Polygon
import iris.exceptions
import iris.warnings
def _extract_relevant_cube_slice(cube, geometry):
"""Calculate geometry intersection with spatial region defined by cube.
This helper method returns the tuple
(subcube, x_coord_of_subcube, y_coord_of_subcube,
(min_x_index, min_y_index, max_x_index, max_y_index)).
If cube and geometry don't overlap, returns None.
"""
# Validate the input parameters
if not cube.coords(axis="x") or not cube.coords(axis="y"):
raise ValueError("The cube must contain x and y axes.")
x_coords = cube.coords(axis="x")
y_coords = cube.coords(axis="y")
if len(x_coords) != 1 or len(y_coords) != 1:
raise ValueError(
"The cube must contain one, and only one, coordinate "
"for each of the x and y axes."
)
x_coord = x_coords[0]
y_coord = y_coords[0]
if not (x_coord.has_bounds() and y_coord.has_bounds()):
raise ValueError("Both horizontal coordinates must have bounds.")
if x_coord.ndim != 1:
raise iris.exceptions.CoordinateMultiDimError(x_coord)
if y_coord.ndim != 1:
raise iris.exceptions.CoordinateMultiDimError(y_coord)
# bounds of cube dimensions
x_bounds = x_coord.bounds
y_bounds = y_coord.bounds
# identify ascending/descending coordinate dimensions
x_ascending = x_coord.points[1] - x_coord.points[0] > 0.0
y_ascending = y_coord.points[1] - y_coord.points[0] > 0.0
# identify upper/lower bounds of coordinate dimensions
x_bounds_lower = x_bounds[:, 0] if x_ascending else x_bounds[:, 1]
y_bounds_lower = y_bounds[:, 0] if y_ascending else y_bounds[:, 1]
x_bounds_upper = x_bounds[:, 1] if x_ascending else x_bounds[:, 0]
y_bounds_upper = y_bounds[:, 1] if y_ascending else y_bounds[:, 0]
# find indices of coordinate bounds to fully cover geometry
x_min_geom, y_min_geom, x_max_geom, y_max_geom = geometry.bounds
try:
x_min_ix = np.where(x_bounds_lower <= x_min_geom)[0]
x_min_ix = x_min_ix[np.argmax(x_bounds_lower[x_min_ix])]
except ValueError:
warnings.warn(
"The geometry exceeds the cube's x dimension at the lower end.",
category=iris.warnings.IrisGeometryExceedWarning,
)
x_min_ix = 0 if x_ascending else x_coord.points.size - 1
try:
x_max_ix = np.where(x_bounds_upper >= x_max_geom)[0]
x_max_ix = x_max_ix[np.argmin(x_bounds_upper[x_max_ix])]
except ValueError:
warnings.warn(
"The geometry exceeds the cube's x dimension at the upper end.",
category=iris.warnings.IrisGeometryExceedWarning,
)
x_max_ix = x_coord.points.size - 1 if x_ascending else 0
try:
y_min_ix = np.where(y_bounds_lower <= y_min_geom)[0]
y_min_ix = y_min_ix[np.argmax(y_bounds_lower[y_min_ix])]
except ValueError:
warnings.warn(
"The geometry exceeds the cube's y dimension at the lower end.",
category=iris.warnings.IrisGeometryExceedWarning,
)
y_min_ix = 0 if y_ascending else y_coord.points.size - 1
try:
y_max_ix = np.where(y_bounds_upper >= y_max_geom)[0]
y_max_ix = y_max_ix[np.argmin(y_bounds_upper[y_max_ix])]
except ValueError:
warnings.warn(
"The geometry exceeds the cube's y dimension at the upper end.",
category=iris.warnings.IrisGeometryExceedWarning,
)
y_max_ix = y_coord.points.size - 1 if y_ascending else 0
# extract coordinate values at these indices
x_min = x_bounds_lower[x_min_ix]
x_max = x_bounds_upper[x_max_ix]
y_min = y_bounds_lower[y_min_ix]
y_max = y_bounds_upper[y_max_ix]
# switch min and max if necessary, to create slice objects later on
if x_min_ix > x_max_ix:
x_min_ix, x_max_ix = x_max_ix, x_min_ix
if y_min_ix > y_max_ix:
y_min_ix, y_max_ix = y_max_ix, y_min_ix
bnds_ix = x_min_ix, y_min_ix, x_max_ix, y_max_ix
# cut the relevant part from the original cube
coord_constr = {
x_coord.name(): lambda x: x_min <= x.point <= x_max,
y_coord.name(): lambda y: y_min <= y.point <= y_max,
}
constraint = iris.Constraint(coord_values=coord_constr)
subcube = cube.extract(constraint)
if subcube is None:
return None
else:
x_coord = subcube.coord(axis="x")
y_coord = subcube.coord(axis="y")
return subcube, x_coord, y_coord, bnds_ix
def geometry_area_weights(cube, geometry, normalize=False):
"""Return the array of weights corresponding to the area of overlap.
Return the array of weights corresponding to the area of overlap between
the cells of cube's horizontal grid, and the given shapely geometry.
The returned array is suitable for use with :const:`iris.analysis.MEAN`.
The cube must have bounded horizontal coordinates.
.. note::
This routine works in Euclidean space. Area calculations do not
account for the curvature of the Earth. And care must be taken to
ensure any longitude values are expressed over a suitable interval.
.. note::
This routine currently does not handle all out-of-bounds cases
correctly. In cases where both the coordinate bounds and the
geometry's bounds lie outside the physically realistic range
(i.e., abs(latitude) > 90., as it is commonly the case when
bounds are constructed via guess_bounds()), the weights
calculation might be wrong. In this case, a UserWarning will
be issued.
.. note::
This function does not maintain laziness when called; it realises data.
See more at :doc:`/userguide/real_and_lazy_data`.
Parameters
----------
cube : :class:`iris.cube.Cube`
A Cube containing a bounded, horizontal grid definition.
geometry : shapely geometry instance
The geometry of interest. To produce meaningful results this geometry
must have a non-zero area. Typically a Polygon or MultiPolygon.
normalize : bool, default=False
Calculate each individual cell weight as the cell area overlap between
the cell and the given shapely geometry divided by the total cell area.
Default is False.
"""
# extract smallest subcube containing geometry
shape = cube.shape
extraction_results = _extract_relevant_cube_slice(cube, geometry)
# test if there is overlap between cube and geometry
if extraction_results is None:
return np.zeros(shape)
subcube, subx_coord, suby_coord, bnds_ix = extraction_results
x_min_ix, y_min_ix, x_max_ix, y_max_ix = bnds_ix
# prepare the weights array
subshape = list(cube.shape)
x_dim = cube.coord_dims(subx_coord)[0]
y_dim = cube.coord_dims(suby_coord)[0]
subshape[x_dim] = subx_coord.shape[0]
subshape[y_dim] = suby_coord.shape[0]
subx_bounds = subx_coord.bounds
suby_bounds = suby_coord.bounds
subweights = np.empty(subshape, np.float32)
# calculate the area weights
for nd_index in np.ndindex(subweights.shape):
xi = nd_index[x_dim]
yi = nd_index[y_dim]
x0, x1 = subx_bounds[xi]
y0, y1 = suby_bounds[yi]
polygon = Polygon([(x0, y0), (x0, y1), (x1, y1), (x1, y0)])
subweights[nd_index] = polygon.intersection(geometry).area
if normalize:
subweights[nd_index] /= polygon.area
# pad the calculated weights with zeros to match original cube shape
weights = np.zeros(shape, np.float32)
slices = []
for i in range(weights.ndim):
if i == x_dim:
slices.append(slice(x_min_ix, x_max_ix + 1))
elif i == y_dim:
slices.append(slice(y_min_ix, y_max_ix + 1))
else:
slices.append(slice(None))
weights[tuple(slices)] = subweights
# Fix for the limitation of iris.analysis.MEAN weights handling.
# Broadcast the array to the full shape of the cube
weights = np.broadcast_arrays(weights, cube.data)[0]
return weights