1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7
8 #include <cstring>
9
10 #include <cdi.h>
11
12 #include "process_int.h"
13 #include "cdo_cdi_wrapper.h"
14 #include <mpim_grid.h>
15
16 /*
17 @Function cdo_define_destagered_grid
18 @Title Define a de-staggered grid for U and V
19
20 @Prototype int cdo_define_destagered_grid(int gridID_u_stag, int gridID_v_stag,
21 double *destagGridOffsets)
22 @Parameter
23 @Item grid_u_stag Staggered grid of u-wind component
24 @Item grid_v_stag Staggered grid of v-wind component
25 @Item grid_uv_destag Destaggered grid of uv-wind
26
27 @Description
28 The function @func{cdo_define_destagered_grid} defines a de-staggered grid for U
29 and V
30
31 @EndFunction
32 */
33 int
cdo_define_destagered_grid(int gridID_u_stag,int gridID_v_stag,double * destagGridOffsets)34 cdo_define_destagered_grid(int gridID_u_stag, int gridID_v_stag, double *destagGridOffsets)
35 {
36 /* Example of horizontal grids (Hirlam LAMH_D11):
37
38 U : lonlat > size : dim = 399300 nlon = 726 nlat = 550
39 rlon : first = -30.15 last = 42.35 inc = 0.1 degrees
40 rlat : first = -30.8 last = 24.1 inc = 0.1 degrees
41 northpole : lon = -195 lat = 30
42 V : lonlat > size : dim = 399300 nlon = 726 nlat = 550
43 rlon : first = -30.2 last = 42.3 inc = 0.1 degrees
44 rlat : first = -30.75 last = 24.15 inc = 0.1 degrees
45 northpole : lon = -195 lat = 30
46 => RESULT:
47 R : lonlat > size : dim = 399300 nlon = 726 nlat = 550
48 rlon : first = -30.2 last = 42.3 inc = 0.1 degrees
49 rlat : first = -30.8 last = 24.1 inc = 0.1 degrees
50 northpole : lon = -195 lat = 30
51 */
52 if (cdoDebugExt)
53 cdo_print("%s(gridID_u=%d,gridID_v=%d,destagGridOffsets(%02.1f,%02.1f)) ...\n", __func__, gridID_u_stag, gridID_v_stag,
54 destagGridOffsets[0], destagGridOffsets[1]);
55
56 if (cdoDebugExt > 1)
57 {
58 cdoPrintGrid(gridID_u_stag, 1);
59 cdoPrintGrid(gridID_v_stag, 1);
60 }
61
62 int gridtype = gridInqType(gridID_u_stag);
63 int xsize = gridInqXsize(gridID_u_stag);
64 int ysize = gridInqYsize(gridID_u_stag);
65
66 double xfirst_U = gridInqXval(gridID_u_stag, 0); // staggered grid of u-wind
67 double yfirst_U = gridInqYval(gridID_u_stag, 0);
68 double xlast_U = gridInqXval(gridID_u_stag, xsize - 1);
69 double ylast_U = gridInqYval(gridID_u_stag, ysize - 1);
70 double xfirst_V = gridInqXval(gridID_v_stag, 0); // staggered grid of v-wind
71 double yfirst_V = gridInqYval(gridID_v_stag, 0);
72 double xlast_V = gridInqXval(gridID_v_stag, xsize - 1);
73 double ylast_V = gridInqYval(gridID_v_stag, ysize - 1);
74 double xinc = gridInqXinc(gridID_u_stag);
75 double yinc = gridInqYinc(gridID_u_stag);
76
77 int gridID_uv_destag = gridDuplicate(gridID_u_stag);
78
79 if (cdoDebugExt)
80 {
81 cdoPrintGrid(gridID_uv_destag, 1);
82
83 cdo_print("%s(): (gridXsize=%d, gridYsize=%d)", __func__, xsize, ysize);
84 cdo_print("%s(): (xfirst_U = %3.2f; yfirst_U = %3.2f); (xfirst_V = %3.2f; yfirst_V = %3.2f)", __func__, xfirst_U, yfirst_U,
85 xfirst_V, yfirst_V);
86 cdo_print("%s(): (xlast_U = %3.2f; ylast_U = %3.2f); (xlast_V = %3.2f; ylast_V = %3.2f)", __func__, xlast_U, ylast_U,
87 xlast_V, ylast_V);
88 }
89
90 double xfirst = 0, xlast = 0, yfirst = 0, ylast = 0;
91 if (IS_EQUAL(destagGridOffsets[0], -0.5) && IS_EQUAL(destagGridOffsets[1], -0.5))
92 {
93 xfirst = xfirst_V;
94 xlast = xlast_V;
95
96 yfirst = yfirst_U;
97 ylast = ylast_U;
98 }
99 else if (IS_EQUAL(destagGridOffsets[0], 0.5) && IS_EQUAL(destagGridOffsets[1], 0.5))
100 {
101 xfirst = xfirst_V + xinc * destagGridOffsets[0];
102 xlast = xlast_V + xinc * destagGridOffsets[0];
103
104 yfirst = yfirst_U + yinc * destagGridOffsets[1];
105 ylast = ylast_U + yinc * destagGridOffsets[1];
106 }
107 else
108 cdo_abort("%s() Unsupported destaggered grid offsets! We support only: (-0.5,-0.5) or (0.5,0.5)", __func__);
109
110 std::vector<double> xvals(xsize);
111 grid_gen_xvals(xsize, xfirst, xlast, xinc, xvals.data());
112 gridDefXvals(gridID_uv_destag, xvals.data());
113
114 std::vector<double> yvals(ysize);
115 grid_gen_yvals(gridtype, ysize, yfirst, ylast, yinc, yvals.data());
116 gridDefYvals(gridID_uv_destag, yvals.data());
117
118 if (cdoDebugExt)
119 {
120 cdo_print("%s():", __func__);
121 cdoPrintGrid(gridID_uv_destag, 1);
122 }
123
124 return gridID_uv_destag;
125 }
126
127 /*
128 @Function cdo_define_sample_grid
129 @Title Define a sampled grid of another grid
130
131 @Prototype int cdo_define_sample_grid(int gridSrcID, int sampleFactor)
132 @Parameter
133 @Item gridSrcID Source grid
134 @Item sampleFactor sampleFactor; typically 2,3,4 ...
135
136 @Description
137 The function @func{cdo_define_sample_grid} defines a sampled grid of another
138 grid
139
140 @EndFunction
141 */
142 int
cdo_define_sample_grid(int gridSrcID,int sampleFactor)143 cdo_define_sample_grid(int gridSrcID, int sampleFactor)
144 {
145 /* Example of horizontal grids (Harmonie HARM36_L25):
146 #
147 # gridID 2
148 #
149 gridtype = projection
150 gridsize = 622521
151 xsize = 789
152 ysize = 789
153 xunits = "m"
154 yunits = "m"
155 xfirst = 0
156 xinc = 2500
157 yfirst = 0
158 yinc = 2500
159 grid_mapping = Lambert_Conformal
160 grid_mapping_name = lambert_conformal_conic
161 standard_parallel = 52.5
162 longitude_of_central_meridian = 0.
163 latitude_of_projection_origin = 52.5
164 longitudeOfFirstGridPointInDegrees = -7.89
165 latitudeOfFirstGridPointInDegrees = 42.935
166 => RESULT:
167 #
168 # gridID 2
169 #
170 gridtype = projection
171 gridsize = 156025
172 xsize = 395
173 ysize = 395
174 xunits = "m"
175 yunits = "m"
176 xfirst = 0
177 xinc = 5000
178 yfirst = 0
179 yinc = 5000
180 grid_mapping = Lambert_Conformal
181 grid_mapping_name = lambert_conformal_conic
182 standard_parallel = 52.5
183 longitude_of_central_meridian = 0.
184 latitude_of_projection_origin = 52.5
185 longitudeOfFirstGridPointInDegrees = -7.89
186 latitudeOfFirstGridPointInDegrees = 42.935
187 */
188 Debug(cdoDebugExt, "%s(gridSrcID=%d, sampleFactor=%d) ...", __func__, gridSrcID, sampleFactor);
189
190 int gridtype = gridInqType(gridSrcID);
191 if (!(gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT || gridtype == GRID_PROJECTION || gridtype == GRID_CURVILINEAR
192 || gridtype == GRID_GENERIC))
193 cdo_abort("Unsupported gridtype: %s", gridNamePtr(gridtype));
194
195 int gridXsize = gridInqXsize(gridSrcID);
196 int gridYsize = gridInqYsize(gridSrcID);
197
198 if ((sampleFactor < 1) || (gridXsize < 1) || (gridYsize < 1) || (sampleFactor > (gridXsize / 4))
199 || (sampleFactor > (gridYsize / 4)))
200 cdo_abort("%s(): Unsupported sampleFactor (%d)! Note that: gridXsize = %d, gridYsize = %d", __func__, sampleFactor, gridXsize,
201 gridYsize);
202
203 if (cdoDebugExt > 20) cdoPrintGrid(gridSrcID, 1);
204
205 int xsize = (gridXsize + (sampleFactor - 1)) / sampleFactor; // HARM36_L25: (789 + 2-1) / 2 = 395
206 int ysize = (gridYsize + (sampleFactor - 1)) / sampleFactor;
207
208 int gridID_sampled = gridCreate(gridtype, xsize * ysize);
209
210 gridDefXsize(gridID_sampled, xsize);
211 gridDefYsize(gridID_sampled, ysize);
212
213 gridDefNP(gridID_sampled, gridInqNP(gridSrcID));
214 int datatype = CDI_UNDEFID;
215 cdiInqKeyInt(gridSrcID, CDI_GLOBAL, CDI_KEY_DATATYPE, &datatype);
216 cdiDefKeyInt(gridID_sampled, CDI_GLOBAL, CDI_KEY_DATATYPE, datatype);
217
218 grid_copy_keys(gridSrcID, gridID_sampled);
219
220 if (gridtype == GRID_PROJECTION) grid_copy_mapping(gridSrcID, gridID_sampled);
221
222 if (gridHasCoordinates(gridSrcID))
223 {
224 if (gridtype == GRID_CURVILINEAR)
225 {
226 std::vector<double> vals(gridXsize * gridYsize);
227 gridInqXvals(gridSrcID, vals.data());
228 double *pvals = vals.data();
229 for (int j = 0; j < gridYsize; j += sampleFactor)
230 for (int i = 0; i < gridXsize; i += sampleFactor) *pvals++ = vals[j * gridXsize + i];
231 gridDefXvals(gridID_sampled, vals.data());
232
233 gridInqYvals(gridSrcID, vals.data());
234 pvals = vals.data();
235 for (int j = 0; j < gridYsize; j += sampleFactor)
236 for (int i = 0; i < gridXsize; i += sampleFactor) *pvals++ = vals[j * gridXsize + i];
237 gridDefYvals(gridID_sampled, vals.data());
238 }
239 else
240 {
241 std::vector<double> xvals(gridXsize);
242 gridInqXvals(gridSrcID, xvals.data());
243 for (int i = 0, j = 0; i < gridXsize; i += sampleFactor) xvals[j++] = xvals[i];
244 gridDefXvals(gridID_sampled, xvals.data());
245
246 std::vector<double> yvals(gridYsize);
247 gridInqYvals(gridSrcID, yvals.data());
248 for (int i = 0, j = 0; i < gridYsize; i += sampleFactor) yvals[j++] = yvals[i];
249 gridDefYvals(gridID_sampled, yvals.data());
250 }
251 }
252
253 if (cdoDebugExt > 20)
254 {
255 cdo_print("cdo SampleGrid: define_sample_grid(): ");
256 cdoPrintGrid(gridID_sampled, 1);
257 }
258
259 return gridID_sampled;
260 }
261
262 /*
263 @Function cdo_define_subgrid_grid
264 @Title Define a sub-grid of another grid (LCC)
265
266 @Prototype int cdo_define_subgrid_grid(int gridIDsrc, int subI0, int subI1, int
267 subJ0, int subJ1)
268 @Parameter
269 @Item gridSrcID Source grid
270 @Item subI0,subI1, subJ0, subJ1 Sub-grid indices
271
272 @Description
273 The function @func{cdo_define_subgrid_grid} defines a sub-grid of another grid
274 (LCC)
275
276 @EndFunction
277 */
278 int
cdo_define_subgrid_grid(int gridSrcID,int subI0,int subI1,int subJ0,int subJ1)279 cdo_define_subgrid_grid(int gridSrcID, int subI0, int subI1, int subJ0, int subJ1)
280 {
281 /* Example of horizontal grids (Harmonie HARM36_L25):
282 #
283 # gridID 2
284 #
285 gridtype = projection
286 gridsize = 622521
287 xsize = 789
288 ysize = 789
289 xunits = "m"
290 yunits = "m"
291 xfirst = 0
292 xinc = 2500
293 yfirst = 0
294 yinc = 2500
295 grid_mapping = Lambert_Conformal
296 grid_mapping_name = lambert_conformal_conic
297 standard_parallel = 52.5
298 longitude_of_central_meridian = 0.
299 latitude_of_projection_origin = 52.5
300 longitudeOfFirstGridPointInDegrees = -7.89
301 latitudeOfFirstGridPointInDegrees = 42.935
302 => RESULT:
303 #
304 # gridID 2
305 #
306 gridtype = projection
307 gridsize = 156025
308 xsize = 350
309 ysize = 350
310 xunits = "m"
311 yunits = "m"
312 xfirst = 0
313 xinc = 2500
314 yfirst = 0
315 yinc = 2500
316 grid_mapping = Lambert_Conformal
317 grid_mapping_name = lambert_conformal_conic
318 standard_parallel = 52.5
319 longitude_of_central_meridian = 0.
320 latitude_of_projection_origin = 52.5
321 longitudeOfFirstGridPointInDegrees = ...
322 latitudeOfFirstGridPointInDegrees = ...
323 */
324 if (cdoDebugExt)
325 cdo_print("%s(gridSrcID=%d, (subI0,subI1,subJ0,subJ1) = (%d,%d,%d,%d) ...", __func__, gridSrcID, subI0, subI1, subJ0, subJ1);
326
327 int gridXsize = gridInqXsize(gridSrcID);
328 int gridYsize = gridInqYsize(gridSrcID);
329 int maxIndexI = gridXsize - 1;
330 int maxIndexJ = gridYsize - 1;
331
332 if ((subI0 < 0) || (subI0 > maxIndexI) || (subI1 <= subI0) || (subI1 > maxIndexI) || (subJ0 < 0) || (subJ0 > maxIndexJ)
333 || (subJ1 <= subJ0) || (subJ1 > maxIndexJ))
334 cdo_abort("%s() Incorrect subgrid specified! (subI0,subI1,subJ0,subJ1) =(%d,%d,%d,%d) Note that: gridXsize=%d, gridYsize=%d",
335 __func__, subI0, subI1, subJ0, subJ1, gridXsize, gridYsize);
336
337 auto gridtype = gridInqType(gridSrcID);
338 if (!(gridtype == GRID_PROJECTION && gridInqProjType(gridSrcID) == CDI_PROJ_LCC))
339 cdo_abort("%s() Error; Only LCC grid is supported; use selindexbox!", __func__);
340
341 struct CDI_GridProjParams gpp;
342 gridInqParamsLCC(gridSrcID, &gpp);
343 gpp.x_0 = gpp.mv;
344 gpp.y_0 = gpp.mv;
345
346 if (cdoDebugExt > 20) cdoPrintGrid(gridSrcID, 1);
347
348 if (cdoDebugExt)
349 {
350 cdo_print("%s() Original LCC grid:", __func__);
351 cdo_print("grid Xsize %d, grid Ysize %d", gridXsize, gridYsize);
352 cdo_print("xval_0 %4.3f, yval_0 %4.3f", gpp.xval_0, gpp.yval_0);
353 }
354
355 const auto gridIDcurvl = gridToCurvilinear(gridSrcID, 1);
356
357 gpp.xval_0 = gridInqXval(gridIDcurvl, 0);
358 gpp.yval_0 = gridInqYval(gridIDcurvl, 0);
359
360 if (cdoDebugExt)
361 {
362 cdo_print("%s() Original LCC grid as curvilinear (with lats-lons computed):", __func__);
363 cdo_print("grid Xsize %zu, grid Ysize %zu", gridInqXsize(gridIDcurvl), gridInqYsize(gridIDcurvl));
364 cdo_print("grid Xfirst %4.3f, grid Yfirst %4.3f", gridInqXval(gridIDcurvl, 0), gridInqYval(gridIDcurvl, 0));
365 cdo_print("grid Xlast %4.3f, grid Ylast %4.3f", gridInqXval(gridIDcurvl, gridInqSize(gridIDcurvl) - 1),
366 gridInqYval(gridIDcurvl, gridInqSize(gridIDcurvl) - 1));
367 cdo_print("xval_0 %4.3f, yval_0 %4.3f", gpp.xval_0, gpp.yval_0);
368 }
369
370 int xsize = subI1 - subI0 + 1;
371 int ysize = subJ1 - subJ0 + 1;
372
373 const auto gridID_sampled = gridCreate(gridtype, xsize * ysize);
374
375 gridDefXsize(gridID_sampled, xsize);
376 gridDefYsize(gridID_sampled, ysize);
377
378 if (gridHasCoordinates(gridSrcID))
379 {
380 std::vector<double> xvals(gridXsize), yvals(gridYsize);
381 gridInqXvals(gridSrcID, xvals.data());
382 gridInqYvals(gridSrcID, yvals.data());
383 gridDefXvals(gridID_sampled, xvals.data());
384 gridDefYvals(gridID_sampled, yvals.data());
385 }
386
387 gridDefNP(gridID_sampled, gridInqNP(gridSrcID));
388 int datatype = CDI_UNDEFID;
389 cdiInqKeyInt(gridSrcID, CDI_GLOBAL, CDI_KEY_DATATYPE, &datatype);
390 cdiDefKeyInt(gridID_sampled, CDI_GLOBAL, CDI_KEY_DATATYPE, datatype);
391
392 grid_copy_keys(gridSrcID, gridID_sampled);
393
394 gpp.xval_0 = gridInqXval(gridIDcurvl, subJ0 * gridXsize + subI0);
395 gpp.yval_0 = gridInqYval(gridIDcurvl, subJ0 * gridXsize + subI0);
396
397 if (cdoDebugExt)
398 {
399 cdo_print("%s() Sub-grid:", __func__);
400 cdo_print("grid Xsize %zu, grid Ysize %zu", gridInqXsize(gridID_sampled), gridInqYsize(gridID_sampled));
401 cdo_print("xval_0 %4.3f, yval_0 %4.3f", gpp.xval_0, gpp.yval_0);
402 }
403
404 gridDefParamsLCC(gridID_sampled, gpp);
405
406 gridDestroy(gridIDcurvl);
407
408 if (cdoDebugExt > 20)
409 {
410 cdo_print("%s(): ", __func__);
411 cdoPrintGrid(gridID_sampled, 1);
412 }
413
414 return gridID_sampled;
415 }
416