-
Notifications
You must be signed in to change notification settings - Fork 95
/
data_reduce.py
307 lines (255 loc) · 10.5 KB
/
data_reduce.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
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Copyright (C) 2010-2021 Pyresample developers
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program. If not, see <http://www.gnu.org/licenses/>.
"""Reduce data sets based on geographical information."""
from __future__ import absolute_import
import numpy as np
# Earth radius
R = 6370997.0
def swath_from_cartesian_grid(cart_grid, lons, lats, data,
radius_of_influence):
"""Make coarse data reduction of swath data by comparison with cartesian grid.
Parameters
----------
chart_grid : numpy array
Grid of area cartesian coordinates
lons : numpy array
Swath lons
lats : numpy array
Swath lats
data : numpy array
Swath data
radius_of_influence : float
Cut off distance in meters
Returns
-------
(lons, lats, data) : list of numpy arrays
Reduced swath data and coordinate set
"""
valid_index = get_valid_index_from_cartesian_grid(cart_grid, lons, lats,
radius_of_influence)
lons = lons[valid_index]
lats = lats[valid_index]
data = data[valid_index]
return lons, lats, data
def get_valid_index_from_cartesian_grid(cart_grid, lons, lats,
radius_of_influence):
"""Calculate relevant data indices using coarse data reduction of swath data by comparison with cartesian grid.
Parameters
----------
chart_grid : numpy array
Grid of area cartesian coordinates
lons : numpy array
Swath lons
lats : numpy array
Swath lats
data : numpy array
Swath data
radius_of_influence : float
Cut off distance in meters
Returns
-------
valid_index : numpy array
Boolean array of same size as lons and lats indicating relevant indices
"""
def _get_lons(x, y):
return np.rad2deg(np.arccos(x / np.sqrt(x ** 2 + y ** 2))) * np.sign(y)
def _get_lats(z):
return 90 - np.rad2deg(np.arccos(z / R))
# Get sides of target grid and transform to lon lats
lons_side1 = _get_lons(cart_grid[0, :, 0], cart_grid[0, :, 1])
lons_side2 = _get_lons(cart_grid[:, -1, 0], cart_grid[:, -1, 1])
lons_side3 = _get_lons(cart_grid[-1, ::-1, 0], cart_grid[-1, ::-1, 1])
lons_side4 = _get_lons(cart_grid[::-1, 0, 0], cart_grid[::-1, 0, 1])
lats_side1 = _get_lats(cart_grid[0, :, 2])
lats_side2 = _get_lats(cart_grid[:, -1, 2])
lats_side3 = _get_lats(cart_grid[-1, ::-1, 2])
lats_side4 = _get_lats(cart_grid[::-1, 0, 2])
valid_index = _get_valid_index(lons_side1, lons_side2, lons_side3, lons_side4,
lats_side1, lats_side2, lats_side3, lats_side4,
lons, lats, radius_of_influence)
return valid_index
def swath_from_lonlat_grid(grid_lons, grid_lats, lons, lats, data,
radius_of_influence):
"""Make coarse data reduction of swath data by comparison with lon lat grid.
Parameters
----------
grid_lons : numpy array
Grid of area lons
grid_lats : numpy array
Grid of area lats
lons : numpy array
Swath lons
lats : numpy array
Swath lats
data : numpy array
Swath data
radius_of_influence : float
Cut off distance in meters
Returns
-------
(lons, lats, data) : list of numpy arrays
Reduced swath data and coordinate set
"""
valid_index = get_valid_index_from_lonlat_grid(
grid_lons, grid_lats, lons, lats, radius_of_influence)
lons = lons[valid_index]
lats = lats[valid_index]
data = data[valid_index]
return lons, lats, data
def swath_from_lonlat_boundaries(boundary_lons, boundary_lats, lons, lats, data,
radius_of_influence):
"""Make coarse data reduction of swath data by comparison with lon lat boundary.
Parameters
----------
boundary_lons : numpy array
Grid of area lons
boundary_lats : numpy array
Grid of area lats
lons : numpy array
Swath lons
lats : numpy array
Swath lats
data : numpy array
Swath data
radius_of_influence : float
Cut off distance in meters
Returns
-------
(lons, lats, data) : list of numpy arrays
Reduced swath data and coordinate set
"""
valid_index = get_valid_index_from_lonlat_boundaries(boundary_lons,
boundary_lats, lons, lats, radius_of_influence)
lons = lons[valid_index]
lats = lats[valid_index]
data = data[valid_index]
return lons, lats, data
def get_valid_index_from_lonlat_grid(grid_lons, grid_lats, lons, lats, radius_of_influence):
"""Calculate relevant data indices using coarse data reduction of swath data by comparison with lon lat grid.
Parameters
----------
chart_grid : numpy array
Grid of area cartesian coordinates
lons : numpy array
Swath lons
lats : numpy array
Swath lats
data : numpy array
Swath data
radius_of_influence : float
Cut off distance in meters
Returns
-------
valid_index : numpy array
Boolean array of same size as lon and lat indicating relevant indices
"""
# Get sides of target grid
lons_side1 = grid_lons[0, :]
lons_side2 = grid_lons[:, -1]
lons_side3 = grid_lons[-1, ::-1]
lons_side4 = grid_lons[::-1, 0]
lats_side1 = grid_lats[0, :]
lats_side2 = grid_lats[:, -1]
lats_side3 = grid_lats[-1, :]
lats_side4 = grid_lats[:, 0]
valid_index = _get_valid_index(lons_side1, lons_side2, lons_side3, lons_side4,
lats_side1, lats_side2, lats_side3, lats_side4,
lons, lats, radius_of_influence)
return valid_index
def get_valid_index_from_lonlat_boundaries(boundary_lons, boundary_lats, lons, lats, radius_of_influence):
"""Find relevant indices from grid boundaries using the winding number theorem."""
valid_index = _get_valid_index(boundary_lons.side1, boundary_lons.side2,
boundary_lons.side3, boundary_lons.side4,
boundary_lats.side1, boundary_lats.side2,
boundary_lats.side3, boundary_lats.side4,
lons, lats, radius_of_influence)
return valid_index
def _get_valid_index(lons_side1, lons_side2, lons_side3, lons_side4,
lats_side1, lats_side2, lats_side3, lats_side4,
lons, lats, radius_of_influence):
"""Find relevant indices from grid boundaries using the winding number theorem."""
# Coarse reduction of data based on extrema analysis of the boundary
# lon lat values of the target grid
illegal_lons = (((lons_side1 < -180) | (lons_side1 > 180)).any() or
((lons_side2 < -180) | (lons_side2 > 180)).any() or
((lons_side3 < -180) | (lons_side3 > 180)).any() or
((lons_side4 < -180) | (lons_side4 > 180)).any())
illegal_lats = (((lats_side1 < -90) | (lats_side1 > 90)).any() or
((lats_side2 < -90) | (lats_side2 > 90)).any() or
((lats_side3 < -90) | (lats_side3 > 90)).any() or
((lats_side4 < -90) | (lats_side4 > 90)).any())
if illegal_lons or illegal_lats:
# Grid boundaries are not safe to operate on
return np.ones(lons.size, dtype=bool)
# Find sum angle sum of grid boundary
angle_sum = 0
for side in (lons_side1, lons_side2, lons_side3, lons_side4):
prev = None
side_sum = 0
for lon in side:
if prev:
delta = lon - prev
if abs(delta) > 180:
delta = (abs(delta) - 360) * (delta // abs(delta))
angle_sum += delta
side_sum += delta
prev = lon
# Buffer min and max lon and lat of interest with radius of interest
lat_min = min(lats_side1.min(), lats_side2.min(), lats_side3.min(),
lats_side4.min())
lat_min_buffered = lat_min - np.degrees(float(radius_of_influence) / R)
lat_max = max(lats_side1.max(), lats_side2.max(), lats_side3.max(),
lats_side4.max())
lat_max_buffered = lat_max + np.degrees(float(radius_of_influence) / R)
max_angle_s2 = max(abs(lats_side2.max()), abs(lats_side2.min()))
max_angle_s4 = max(abs(lats_side4.max()), abs(lats_side4.min()))
lon_min_buffered = (lons_side4.min() -
np.degrees(float(radius_of_influence) /
(np.sin(np.radians(max_angle_s4)) * R)))
lon_max_buffered = (lons_side2.max() +
np.degrees(float(radius_of_influence) /
(np.sin(np.radians(max_angle_s2)) * R)))
# From the winding number theorem follows:
# angle_sum possiblilities:
# -360: area covers north pole
# 360: area covers south pole
# 0: area covers no poles
# else: area covers both poles
if round(angle_sum) == -360:
# Covers NP
valid_index = (lats >= lat_min_buffered)
elif round(angle_sum) == 360:
# Covers SP
valid_index = (lats <= lat_max_buffered)
elif round(angle_sum) == 0:
# Covers no poles
valid_lats = (lats >= lat_min_buffered) * (lats <= lat_max_buffered)
if lons_side2.min() > lons_side4.max():
# No date line crossing
valid_lons = (lons >= lon_min_buffered) * \
(lons <= lon_max_buffered)
else:
# Date line crossing
seg1 = (lons >= lon_min_buffered) * (lons <= 180)
seg2 = (lons <= lon_max_buffered) * (lons >= -180)
valid_lons = seg1 + seg2
valid_index = valid_lats * valid_lons
else:
# Covers both poles don't reduce
valid_index = np.ones(lons.size, dtype=bool)
return valid_index