|
| 1 | +import datetime |
| 2 | +import unittest |
| 3 | + |
1 | 4 | import numpy as np |
2 | 5 |
|
3 | 6 | import cf |
4 | 7 |
|
5 | | -f = cf.example_field(6) |
6 | | -f.del_construct("auxiliarycoordinate3") |
7 | | -f.del_construct("grid_mapping_name:latitude_longitude") |
8 | | -g = f.copy() |
9 | | -g.dump() |
10 | 8 |
|
| 9 | +class WeightsTest(unittest.TestCase): |
| 10 | + def test_weights_polygon_area_geometry(self): |
| 11 | + f = cf.example_field(6) |
| 12 | + f.del_construct("auxiliarycoordinate3") |
| 13 | + f.del_construct("grid_mapping_name:latitude_longitude") |
| 14 | + |
| 15 | + # Surface area of unit sphere |
| 16 | + sphere_area = 4 * np.pi |
| 17 | + correct_weights = np.array([sphere_area / 8, sphere_area / 4]) |
| 18 | + |
| 19 | + # Spherical polygon geometry weights with duplicated first/last |
| 20 | + # node |
| 21 | + lon = cf.AuxiliaryCoordinate() |
| 22 | + lon.standard_name = "longitude" |
| 23 | + bounds = cf.Data( |
| 24 | + [[315, 45, 45, 315, 999, 999, 999], [90, 90, 0, 45, 45, 135, 90]], |
| 25 | + "degrees_east", |
| 26 | + mask_value=999, |
| 27 | + ).reshape(2, 1, 7) |
| 28 | + lon.set_bounds(cf.Bounds(data=bounds)) |
| 29 | + lon.set_geometry("polygon") |
| 30 | + |
| 31 | + lat = cf.AuxiliaryCoordinate() |
| 32 | + lat.standard_name = "latitude" |
| 33 | + bounds = cf.Data( |
| 34 | + [[0, 0, 90, 0, 999, 999, 999], [0, 90, 0, 0, -90, 0, 0]], |
| 35 | + "degrees_north", |
| 36 | + mask_value=999, |
| 37 | + ).reshape(2, 1, 7) |
| 38 | + lat.set_bounds(cf.Bounds(data=bounds)) |
| 39 | + lat.set_geometry("polygon") |
| 40 | + |
| 41 | + f.del_construct("longitude") |
| 42 | + f.del_construct("latitude") |
| 43 | + |
| 44 | + f.set_construct(lon, axes="domainaxis0", copy=False) |
| 45 | + f.set_construct(lat, axes="domainaxis0", copy=False) |
| 46 | + |
| 47 | + w = f.weights("X", great_circle=True) |
| 48 | + self.assertTrue((w.array == correct_weights).all()) |
| 49 | + |
| 50 | + r = 2 |
| 51 | + w = f.weights("area", great_circle=True, measure=True, radius=r) |
| 52 | + self.assertTrue((w.array == (r**2) * correct_weights).all()) |
| 53 | + |
| 54 | + # Spherical polygon geometry weights without duplicated |
| 55 | + # first/last node |
| 56 | + bounds = cf.Data( |
| 57 | + [[315, 45, 45, 999, 999, 999], [90, 90, 0, 45, 45, 135]], |
| 58 | + "degrees_east", |
| 59 | + mask_value=999, |
| 60 | + ).reshape(2, 1, 6) |
| 61 | + lon.set_bounds(cf.Bounds(data=bounds)) |
| 62 | + |
| 63 | + bounds = cf.Data( |
| 64 | + [[0, 0, 90, 999, 999, 999], [0, 90, 0, 0, -90, 0]], |
| 65 | + "degrees_north", |
| 66 | + mask_value=999, |
| 67 | + ).reshape(2, 1, 6) |
| 68 | + lat.set_bounds(cf.Bounds(data=bounds)) |
| 69 | + |
| 70 | + w = f.weights("X", great_circle=True) |
| 71 | + self.assertTrue((w.array == correct_weights).all()) |
| 72 | + |
| 73 | + r = 2 |
| 74 | + w = f.weights("area", great_circle=True, measure=True, radius=r) |
| 75 | + self.assertTrue((w.array == (r**2) * correct_weights).all()) |
| 76 | + |
| 77 | + # Plane polygon geometry weights with no duplicated first/last |
| 78 | + # nodes, and an interior ring |
| 79 | + lon.override_units("m", inplace=True) |
| 80 | + lon.standard_name = "projection_x_coordinate" |
| 81 | + bounds = cf.Data( |
| 82 | + [ |
| 83 | + [ |
| 84 | + [2, 2, 0, 0, 999, 999, 999, 999], # anticlockwise |
| 85 | + [0.5, 1.5, 1.5, 0.5, 999, 999, 999, 999], |
| 86 | + ], # clockwise |
| 87 | + [ |
| 88 | + [2, 2, 0, 0, 1, 1, 3, 3], |
| 89 | + [999, 999, 999, 999, 999, 999, 999, 999], |
| 90 | + ], |
| 91 | + ], |
| 92 | + "m", |
| 93 | + mask_value=999, |
| 94 | + ).reshape(2, 2, 8) |
| 95 | + lon.set_bounds(cf.Bounds(data=bounds)) |
| 96 | + lon.set_interior_ring(cf.InteriorRing(data=[[0, 1], [0, 0]])) |
| 97 | + |
| 98 | + lat.override_units("m", inplace=True) |
| 99 | + lat.standard_name = "projection_y_coordinate" |
| 100 | + bounds = cf.Data( |
| 101 | + [ |
| 102 | + [ |
| 103 | + [-1, 1, 1, -1, 999, 999, 999, 999], # anticlockwise |
| 104 | + [0.5, 0.5, -0.5, -0.5, 999, 999, 999, 999], |
| 105 | + ], # clockwise |
| 106 | + [ |
| 107 | + [-1, 1, 1, -1, -1, -3, -3, -1], |
| 108 | + [999, 999, 999, 999, 999, 999, 999, 999], |
| 109 | + ], |
| 110 | + ], |
| 111 | + "m", |
| 112 | + mask_value=999, |
| 113 | + ).reshape(2, 2, 8) |
| 114 | + lat.set_bounds(cf.Bounds(data=bounds)) |
| 115 | + lat.set_interior_ring(cf.InteriorRing(data=[[0, 1], [0, 0]])) |
| 116 | + |
| 117 | + correct_weights = np.array([3, 8]) |
| 118 | + w = f.weights("area") |
| 119 | + self.assertTrue((w.array == correct_weights).all()) |
| 120 | + |
| 121 | + def test_weights_polygon_area_ugrid(self): |
| 122 | + f = cf.example_field(8) |
| 123 | + f = f[..., [0, 2]] |
11 | 124 |
|
12 | | -lon = cf.AuxiliaryCoordinate() |
13 | | -lon.standard_name = "longitude" |
14 | | -bounds = cf.Data( |
15 | | - [[315, 45, 45, 315, 999, 999, 999], [90, 90, 0, 45, 45, 135, 90]], |
16 | | - "degrees_east", |
17 | | - masked_values=999, |
18 | | -).reshape(2, 1, 7) |
19 | | -lon.set_bounds(cf.Bounds(data=bounds)) |
20 | | -lon.set_geometry("polygon") |
| 125 | + # Surface area of unit sphere |
| 126 | + sphere_area = 4 * np.pi |
| 127 | + correct_weights = np.array([sphere_area / 8, sphere_area / 4]) |
21 | 128 |
|
22 | | -lat = cf.AuxiliaryCoordinate() |
23 | | -lat.standard_name = "latitude" |
24 | | -bounds = cf.Data( |
25 | | - [[0, 0, 90, 0, 999, 999, 999], [0, 90, 0, 0, -90, 0, 0]], |
26 | | - "degrees_north", |
27 | | - masked_values=999, |
28 | | -).reshape(2, 1, 7) |
29 | | -lat.set_bounds(cf.Bounds(data=bounds)) |
30 | | -lat.set_geometry("polygon") |
31 | | -lon.dump() |
32 | | -f.del_construct("longitude") |
33 | | -f.del_construct("latitude") |
| 129 | + # Spherical polygon weights |
| 130 | + lon = f.auxiliary_coordinate("X") |
| 131 | + lon.del_data() |
| 132 | + bounds = cf.Data( |
| 133 | + [[315, 45, 45, 999, 999, 999], [90, 90, 0, 45, 45, 135]], |
| 134 | + "degrees_east", |
| 135 | + mask_value=999, |
| 136 | + ).reshape(2, 6) |
| 137 | + lon.set_bounds(cf.Bounds(data=bounds)) |
34 | 138 |
|
35 | | -# f.dump() |
36 | | -f.set_construct(lon, axes="domainaxis0") |
37 | | -f.set_construct(lat, axes="domainaxis0") |
38 | | -print(f.constructs) |
| 139 | + lat = f.auxiliary_coordinate("Y") |
| 140 | + bounds = cf.Data( |
| 141 | + [[0, 0, 90, 999, 999, 999], [0, 90, 0, 0, -90, 0]], |
| 142 | + "degrees_north", |
| 143 | + mask_value=999, |
| 144 | + ).reshape(2, 6) |
| 145 | + lat.set_bounds(cf.Bounds(data=bounds)) |
39 | 146 |
|
40 | | -f.dump() |
| 147 | + w = f.weights("X", great_circle=True) |
| 148 | + self.assertTrue((w.array == correct_weights).all()) |
41 | 149 |
|
42 | | -sphere_area = 4 * np.pi |
| 150 | + r = 2 |
| 151 | + w = f.weights("area", great_circle=True, measure=True, radius=r) |
| 152 | + self.assertTrue((w.array == (r**2) * correct_weights).all()) |
43 | 153 |
|
| 154 | + # Plane polygon weights |
| 155 | + lon.override_units("m", inplace=True) |
| 156 | + lon.standard_name = "projection_x_coordinate" |
| 157 | + bounds = cf.Data( |
| 158 | + [[2, 2, 0, 0, 999, 999, 999, 999], [2, 2, 0, 0, 1, 1, 3, 3]], |
| 159 | + "m", |
| 160 | + mask_value=999, |
| 161 | + ).reshape(2, 8) |
| 162 | + lon.set_bounds(cf.Bounds(data=bounds)) |
44 | 163 |
|
45 | | -w = f.weights("X", great_circle=True) |
46 | | -print(repr(w), w.array, np.array([sphere_area / 8, sphere_area / 4])) |
| 164 | + lat.override_units("m", inplace=True) |
| 165 | + lat.standard_name = "projection_y_coordinate" |
| 166 | + bounds = cf.Data( |
| 167 | + [ |
| 168 | + [-1, 1, 1, -1, 999, 999, 999, 999], |
| 169 | + [-1, 1, 1, -1, -1, -3, -3, -1], |
| 170 | + ], |
| 171 | + "m", |
| 172 | + mask_value=999, |
| 173 | + ).reshape(2, 8) |
| 174 | + lat.set_bounds(cf.Bounds(data=bounds)) |
47 | 175 |
|
48 | | -print(w.array == np.array([sphere_area / 8, sphere_area / 4])) |
| 176 | + correct_weights = np.array([4, 8]) |
| 177 | + w = f.weights("area") |
| 178 | + self.assertTrue((w.array == correct_weights).all()) |
49 | 179 |
|
50 | | -w = f.weights("X", great_circle=True, measure=True, radius=2) |
51 | | -print(repr(w), w.array, 4 * np.array([sphere_area / 8, sphere_area / 4])) |
52 | 180 |
|
53 | | -print(w.array - 4 * np.array([sphere_area / 8, sphere_area / 4])) |
| 181 | +if __name__ == "__main__": |
| 182 | + print("Run date:", datetime.datetime.now()) |
| 183 | + cf.environment() |
| 184 | + print("") |
| 185 | + unittest.main(verbosity=2) |
0 commit comments