|
12 | 12 | )
|
13 | 13 |
|
14 | 14 | ds_no_attrs = airds.copy(deep=True)
|
15 |
| -for variable in ds_no_attrs.variables: |
16 |
| - ds_no_attrs[variable].attrs = {} |
| 15 | +for _variable in ds_no_attrs.variables: |
| 16 | + ds_no_attrs[_variable].attrs = {} |
17 | 17 |
|
18 | 18 |
|
19 | 19 | # POM dataset
|
|
154 | 154 | )
|
155 | 155 |
|
156 | 156 |
|
157 |
| -# Dataset with random data on a grid that is some sort of Mollweide projection |
158 |
| -XX, YY = np.mgrid[:11, :11] * 5 - 25 |
159 |
| -XX_bnds, YY_bnds = np.mgrid[:12, :12] * 5 - 27.5 |
| 157 | +def _create_mollw_ds(): |
| 158 | + # Dataset with random data on a grid that is some sort of Mollweide projection |
| 159 | + XX, YY = np.mgrid[:11, :11] * 5 - 25 |
| 160 | + XX_bnds, YY_bnds = np.mgrid[:12, :12] * 5 - 27.5 |
160 | 161 |
|
161 |
| -R = 50 |
162 |
| -theta = np.arcsin(YY / (R * np.sqrt(2))) |
163 |
| -lat = np.rad2deg(np.arcsin((2 * theta + np.sin(2 * theta)) / np.pi)) |
164 |
| -lon = np.rad2deg(XX * np.pi / (R * 2 * np.sqrt(2) * np.cos(theta))) |
| 162 | + R = 50 |
| 163 | + theta = np.arcsin(YY / (R * np.sqrt(2))) |
| 164 | + lat = np.rad2deg(np.arcsin((2 * theta + np.sin(2 * theta)) / np.pi)) |
| 165 | + lon = np.rad2deg(XX * np.pi / (R * 2 * np.sqrt(2) * np.cos(theta))) |
165 | 166 |
|
166 |
| -theta_bnds = np.arcsin(YY_bnds / (R * np.sqrt(2))) |
167 |
| -lat_vertices = np.rad2deg(np.arcsin((2 * theta_bnds + np.sin(2 * theta_bnds)) / np.pi)) |
168 |
| -lon_vertices = np.rad2deg(XX_bnds * np.pi / (R * 2 * np.sqrt(2) * np.cos(theta_bnds))) |
169 |
| - |
170 |
| -lon_bounds = np.stack( |
171 |
| - ( |
172 |
| - lon_vertices[:-1, :-1], |
173 |
| - lon_vertices[:-1, 1:], |
174 |
| - lon_vertices[1:, 1:], |
175 |
| - lon_vertices[1:, :-1], |
176 |
| - ), |
177 |
| - axis=0, |
178 |
| -) |
179 |
| -lat_bounds = np.stack( |
180 |
| - ( |
181 |
| - lat_vertices[:-1, :-1], |
182 |
| - lat_vertices[:-1, 1:], |
183 |
| - lat_vertices[1:, 1:], |
184 |
| - lat_vertices[1:, :-1], |
185 |
| - ), |
186 |
| - axis=0, |
187 |
| -) |
| 167 | + theta_bnds = np.arcsin(YY_bnds / (R * np.sqrt(2))) |
| 168 | + lat_vertices = np.rad2deg( |
| 169 | + np.arcsin((2 * theta_bnds + np.sin(2 * theta_bnds)) / np.pi) |
| 170 | + ) |
| 171 | + lon_vertices = np.rad2deg( |
| 172 | + XX_bnds * np.pi / (R * 2 * np.sqrt(2) * np.cos(theta_bnds)) |
| 173 | + ) |
188 | 174 |
|
189 |
| -mollwds = xr.Dataset( |
190 |
| - coords=dict( |
191 |
| - lon=xr.DataArray( |
192 |
| - lon, |
193 |
| - dims=("x", "y"), |
194 |
| - attrs={"units": "degrees_east", "bounds": "lon_bounds"}, |
| 175 | + lon_bounds = np.stack( |
| 176 | + ( |
| 177 | + lon_vertices[:-1, :-1], |
| 178 | + lon_vertices[:-1, 1:], |
| 179 | + lon_vertices[1:, 1:], |
| 180 | + lon_vertices[1:, :-1], |
195 | 181 | ),
|
196 |
| - lat=xr.DataArray( |
197 |
| - lat, |
198 |
| - dims=("x", "y"), |
199 |
| - attrs={"units": "degrees_north", "bounds": "lat_bounds"}, |
| 182 | + axis=0, |
| 183 | + ) |
| 184 | + lat_bounds = np.stack( |
| 185 | + ( |
| 186 | + lat_vertices[:-1, :-1], |
| 187 | + lat_vertices[:-1, 1:], |
| 188 | + lat_vertices[1:, 1:], |
| 189 | + lat_vertices[1:, :-1], |
200 | 190 | ),
|
201 |
| - ), |
202 |
| - data_vars=dict( |
203 |
| - lon_bounds=xr.DataArray( |
204 |
| - lon_bounds, dims=("bounds", "x", "y"), attrs={"units": "degrees_east"} |
| 191 | + axis=0, |
| 192 | + ) |
| 193 | + |
| 194 | + mollwds = xr.Dataset( |
| 195 | + coords=dict( |
| 196 | + lon=xr.DataArray( |
| 197 | + lon, |
| 198 | + dims=("x", "y"), |
| 199 | + attrs={"units": "degrees_east", "bounds": "lon_bounds"}, |
| 200 | + ), |
| 201 | + lat=xr.DataArray( |
| 202 | + lat, |
| 203 | + dims=("x", "y"), |
| 204 | + attrs={"units": "degrees_north", "bounds": "lat_bounds"}, |
| 205 | + ), |
205 | 206 | ),
|
206 |
| - lat_bounds=xr.DataArray( |
207 |
| - lat_bounds, dims=("bounds", "x", "y"), attrs={"units": "degrees_north"} |
| 207 | + data_vars=dict( |
| 208 | + lon_bounds=xr.DataArray( |
| 209 | + lon_bounds, dims=("bounds", "x", "y"), attrs={"units": "degrees_east"} |
| 210 | + ), |
| 211 | + lat_bounds=xr.DataArray( |
| 212 | + lat_bounds, dims=("bounds", "x", "y"), attrs={"units": "degrees_north"} |
| 213 | + ), |
| 214 | + lon_vertices=xr.DataArray(lon_vertices, dims=("x_vertices", "y_vertices")), |
| 215 | + lat_vertices=xr.DataArray(lat_vertices, dims=("x_vertices", "y_vertices")), |
208 | 216 | ),
|
209 |
| - lon_vertices=xr.DataArray(lon_vertices, dims=("x_vertices", "y_vertices")), |
210 |
| - lat_vertices=xr.DataArray(lat_vertices, dims=("x_vertices", "y_vertices")), |
211 |
| - ), |
212 |
| -) |
| 217 | + ) |
| 218 | + |
| 219 | + return mollwds |
| 220 | + |
| 221 | + |
| 222 | +mollwds = _create_mollw_ds() |
213 | 223 |
|
214 | 224 | forecast = xr.decode_cf(
|
215 | 225 | xr.Dataset.from_dict(
|
|
0 commit comments