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 <algorithm> // sort
9
10 #include <cdi.h>
11
12 #include "cdo_options.h"
13 #include "dmemory.h"
14 #include "process_int.h"
15 #include "cdi_lockedIO.h"
16 #include <mpim_grid.h>
17 #include "grid_point_search.h"
18 #include "verifygrid.h"
19
20 constexpr int MAX_CHILDS = 9;
21
22 struct CellIndex
23 {
24 long ncells;
25 long *neighbor; // neighbor cell index
26 long *parent; // parent cell index
27 long *child; // child cell index
28 const char *filename;
29 };
30
31 static void
copy_data_to_index(long ncells,const Varray<double> & data,long * cellindex)32 copy_data_to_index(long ncells, const Varray<double> &data, long *cellindex)
33 {
34 for (long i = 0; i < ncells; ++i) cellindex[i] = std::lround(data[i]);
35 }
36
37 static void
free_cellindex(CellIndex * cellindex)38 free_cellindex(CellIndex *cellindex)
39 {
40 if (cellindex->neighbor) Free(cellindex->neighbor);
41 if (cellindex->parent) Free(cellindex->parent);
42 if (cellindex->child) Free(cellindex->child);
43 Free(cellindex);
44 }
45
46 static CellIndex *
read_cellindex(const std::string & filename)47 read_cellindex(const std::string &filename)
48 {
49 const auto streamID = stream_open_read_locked(filename.c_str());
50 const auto vlistID = streamInqVlist(streamID);
51 const auto ngrids = vlistNgrids(vlistID);
52 int gridID = -1;
53 for (int index = 0; index < ngrids; ++index)
54 {
55 gridID = vlistGrid(vlistID, index);
56 if (gridInqType(gridID) == GRID_UNSTRUCTURED && gridInqNvertex(gridID) == 3) break;
57 }
58
59 if (gridID == -1) cdo_abort("No ICON grid found in %s!", filename);
60
61 // int nid = CDI_UNDEFID;
62 int pid = CDI_UNDEFID;
63 // int cid = CDI_UNDEFID;
64 const auto nvars = vlistNvars(vlistID);
65 char varname[CDI_MAX_NAME];
66 for (int varID = 0; varID < nvars; ++varID)
67 {
68 vlistInqVarName(vlistID, varID, varname);
69 /*
70 if (strcmp(varname, "neighbor_cell_index") == 0)
71 {
72 nid = varID;
73 break;
74 }
75 */
76 if (strcmp(varname, "parent_cell_index") == 0)
77 {
78 pid = varID;
79 break;
80 }
81 /*
82 if ( strcmp(varname, "child_cell_index") == 0 )
83 {
84 cid = varID;
85 break;
86 }
87 */
88 }
89
90 // if (nid == CDI_UNDEFID) cdo_abort("neighbor_cell_index not found in %s!", filename);
91 // if (pid == CDI_UNDEFID) cdo_abort("parent_cell_index not found in %s!", filename);
92 // if (cid == CDI_UNDEFID) cdo_abort("child_cell_index not found in %s!", filename);
93
94 const long ncells = gridInqSize(gridID);
95
96 CellIndex *cellindex = (CellIndex *) Malloc(sizeof(CellIndex));
97 cellindex->ncells = ncells;
98
99 cellindex->neighbor = nullptr;
100 // cellindex->neighbor = (long*) Malloc(3*ncells*sizeof(long));
101 cellindex->parent = (long *) Malloc(ncells * sizeof(long));
102 cellindex->child = nullptr;
103 // cellindex->child = (cid != CDI_UNDEFID) ? (int*) Malloc(MAX_CHILDS*ncells*sizeof(int)) : nullptr;
104 Varray<double> data(ncells);
105
106 for (long i = 0; i < ncells; ++i) cellindex->parent[i] = 0;
107
108 const auto nrecs = streamInqTimestep(streamID, 0);
109 for (int recID = 0; recID < nrecs; recID++)
110 {
111 int varID, levelID;
112 size_t nmiss;
113 streamInqRecord(streamID, &varID, &levelID);
114 if (varID == pid /* || varID == nid || varID == cid */)
115 {
116 streamReadRecord(streamID, data.data(), &nmiss);
117 // if (varID == pid)
118 {
119 if (Options::cdoVerbose) cdo_print("Read parent_cell_index");
120 copy_data_to_index(ncells, data, cellindex->parent);
121 }
122 // else if ( varID == nid ) copy_data_to_index(ncells, data, cellindex->neighbor+levelID*ncells);
123 // else if ( varID == cid ) copy_data_to_index(ncells, data, cellindex->child+levelID*ncells);
124 }
125 }
126
127 // Fortran to C index
128 for (long i = 0; i < ncells; ++i) cellindex->parent[i] -= 1;
129 // for ( long i = 0; i < 3*ncells; ++i ) cellindex->neighbor[i] -= 1;
130
131 streamClose(streamID);
132
133 return cellindex;
134 }
135
136 static int
read_grid(const char * filename)137 read_grid(const char *filename)
138 {
139 const auto streamID = stream_open_read_locked(filename);
140 const auto vlistID = streamInqVlist(streamID);
141 const auto ngrids = vlistNgrids(vlistID);
142 int gridID = -1;
143 for (int index = 0; index < ngrids; ++index)
144 {
145 gridID = vlistGrid(vlistID, index);
146 if (gridInqType(gridID) == GRID_UNSTRUCTURED && gridInqNvertex(gridID) == 3) break;
147 }
148
149 if (gridID == -1) cdo_abort("No ICON grid found in %s!", filename);
150
151 const auto gridID2 = gridDuplicate(gridID);
152
153 streamClose(streamID);
154
155 return gridID2;
156 }
157
158 /**
159 * Return the first index of element x fits.
160 *
161 * If no interval can be found return -1.
162
163 * @param *array ascending sorted list
164 * @param n length of the sorted list
165 * @param search the element to find a position for
166 */
167 static long
find_index(int search,long n,const long * array)168 find_index(int search, long n, const long *array)
169 {
170 long first = 0;
171 long last = n - 1;
172 long middle = (first + last) / 2;
173
174 while (first <= last)
175 {
176 if (array[middle] < search)
177 first = middle + 1;
178 else if (array[middle] == search)
179 {
180 for (long i = middle; i >= 0; i--)
181 {
182 if (array[i] == search)
183 middle = i;
184 else
185 break;
186 }
187 return middle;
188 }
189 else
190 last = middle - 1;
191
192 middle = (first + last) / 2;
193 }
194
195 return -1;
196 }
197
198 struct SortInfo
199 {
200 int p, i;
201 };
202
203 static bool
cmpsinfo(const SortInfo & a,const SortInfo & b)204 cmpsinfo(const SortInfo &a, const SortInfo &b)
205 {
206 return a.p < b.p;
207 }
208
209 static void
compute_child_from_parent(CellIndex * cellindex1,CellIndex * cellindex2)210 compute_child_from_parent(CellIndex *cellindex1, CellIndex *cellindex2)
211 {
212 if (Options::cdoVerbose) cdo_print("%s", __func__);
213
214 const auto ncells1 = cellindex1->ncells;
215 long *parent1 = cellindex1->parent;
216
217 std::vector<long> idx1(ncells1);
218 for (long i = 0; i < ncells1; ++i) idx1[i] = i;
219 for (long i = 1; i < ncells1; ++i)
220 if (parent1[i] < parent1[i - 1])
221 {
222 if (Options::cdoVerbose) cdo_print("Sort parent index of %s!", cellindex1->filename);
223 std::vector<SortInfo> sinfo(ncells1);
224 for (long j = 0; j < ncells1; ++j)
225 {
226 sinfo[j].p = parent1[j];
227 sinfo[j].i = idx1[j];
228 }
229 std::sort(sinfo.begin(), sinfo.end(), cmpsinfo);
230 for (long j = 0; j < ncells1; ++j)
231 {
232 parent1[j] = sinfo[j].p;
233 idx1[j] = sinfo[j].i;
234 }
235 break;
236 }
237
238 const auto ncells2 = cellindex2->ncells;
239 long *child2 = (long *) Malloc(MAX_CHILDS * ncells2 * sizeof(long));
240 cellindex2->child = child2;
241 for (long i = 0; i < ncells2; ++i)
242 {
243 for (long k = 0; k < MAX_CHILDS; ++k) child2[i * MAX_CHILDS + k] = -1;
244 long j = find_index(i, ncells1, parent1);
245 if (j < 0) continue;
246 for (long k = 0; k < MAX_CHILDS; ++k)
247 {
248 if (i != parent1[j + k]) break;
249 // child2[i*MAX_CHILDS+k] = j+k;
250 child2[i * MAX_CHILDS + k] = idx1[j + k];
251 }
252 // if ( i%10000 == 0 ) printf("%d %d %d %d %d %d\n", i, j, parent1[j], parent1[j+1], parent1[j+2], parent1[j+3]);
253 }
254 }
255
256 static void
read_coordinates(const char * filename,long n,double * lon,double * lat,int nv,double * lon_bnds,double * lat_bnds)257 read_coordinates(const char *filename, long n, double *lon, double *lat, int nv, double *lon_bnds, double *lat_bnds)
258 {
259 const auto streamID = streamOpenRead(filename);
260 const auto vlistID = streamInqVlist(streamID);
261 const auto ngrids = vlistNgrids(vlistID);
262 int gridID = -1;
263 for (int index = 0; index < ngrids; ++index)
264 {
265 gridID = vlistGrid(vlistID, index);
266 if (gridInqType(gridID) == GRID_UNSTRUCTURED && (long) gridInqSize(gridID) == n && gridInqNvertex(gridID) == 3) break;
267 }
268
269 if (gridID == -1) cdo_abort("No ICON grid with %ld cells found in %s!", n, filename);
270
271 gridInqXvals(gridID, lon);
272 gridInqYvals(gridID, lat);
273
274 // Convert lat/lon units if required
275 cdo_grid_to_radian(gridID, CDI_XAXIS, n, lon, "grid center lon");
276 cdo_grid_to_radian(gridID, CDI_YAXIS, n, lat, "grid center lat");
277
278 if (nv == 3 && lon_bnds && lat_bnds)
279 {
280 gridInqXbounds(gridID, lon_bnds);
281 gridInqYbounds(gridID, lat_bnds);
282
283 cdo_grid_to_radian(gridID, CDI_XAXIS, n * 3, lon_bnds, "grid corner lon");
284 cdo_grid_to_radian(gridID, CDI_YAXIS, n * 3, lat_bnds, "grid corner lat");
285 }
286
287 streamClose(streamID);
288 }
289
290 constexpr int MAX_SEARCH = 128; // the triangles are distorted!
291
292 static void
compute_child_from_bounds(CellIndex * cellindex2,long ncells2,double * grid_center_lon2,double * grid_center_lat2,double * grid_corner_lon2,double * grid_corner_lat2,long ncells1,const Varray<double> & grid_center_lon1,const Varray<double> & grid_center_lat1)293 compute_child_from_bounds(CellIndex *cellindex2, long ncells2, double *grid_center_lon2, double *grid_center_lat2,
294 double *grid_corner_lon2, double *grid_corner_lat2, long ncells1, const Varray<double> &grid_center_lon1,
295 const Varray<double> &grid_center_lat1)
296 {
297 if (Options::cdoVerbose) cdo_print("%s", __func__);
298
299 bool xIsCyclic = false;
300 size_t dims[2];
301 dims[0] = ncells1;
302 dims[1] = 0;
303 GridPointSearch gps;
304 grid_point_search_create(gps, xIsCyclic, dims, ncells1, grid_center_lon1, grid_center_lat1);
305 knnWeightsType knnWeights(MAX_SEARCH);
306 size_t *nbr_addr = &knnWeights.m_addr[0];
307
308 int ncorner = 3;
309 Point3D center_point_xyz;
310 Varray<Point3D> cell_corners_xyz(4);
311 Point center_point_plane_projection;
312 Varray<Point> cell_corners_plane_projection(4);
313
314 long *child2 = (long *) Malloc(MAX_CHILDS * ncells2 * sizeof(long));
315 cellindex2->child = child2;
316 for (long cell_no2 = 0; cell_no2 < ncells2; ++cell_no2)
317 {
318 for (int k = 0; k < MAX_CHILDS; ++k) child2[cell_no2 * MAX_CHILDS + k] = -1;
319
320 set_cell_corners_xyz(ncorner, &grid_corner_lon2[cell_no2 * ncorner], &grid_corner_lat2[cell_no2 * ncorner], cell_corners_xyz);
321 cell_corners_xyz[ncorner] = cell_corners_xyz[0];
322
323 const auto coordinate_to_ignore = find_coordinate_to_ignore(cell_corners_xyz);
324
325 bool invert_result = false;
326 const auto cval = (coordinate_to_ignore == 1) ? cell_corners_xyz[0].X
327 : ((coordinate_to_ignore == 2) ? cell_corners_xyz[0].Y : cell_corners_xyz[0].Z);
328 if (cval < 0.0) invert_result = true;
329
330 set_cell_corners_plane_projection(coordinate_to_ignore, ncorner, cell_corners_xyz, cell_corners_plane_projection);
331
332 const auto polygon_area = calculate_the_polygon_area(cell_corners_plane_projection, ncorner + 1);
333 auto is_clockwise = are_polygon_vertices_arranged_in_clockwise_order(polygon_area);
334
335 if (invert_result) is_clockwise = !is_clockwise;
336 if (is_clockwise) continue;
337
338 grid_search_point(gps, grid_center_lon2[cell_no2], grid_center_lat2[cell_no2], knnWeights);
339 int k = 0;
340
341 double center_coordinates[3];
342 for (int i = 0; i < MAX_SEARCH; ++i)
343 {
344 const auto cell_no1 = nbr_addr[i];
345 if (cell_no1 < SIZE_MAX)
346 {
347 gcLLtoXYZ(grid_center_lon1[cell_no1], grid_center_lat1[cell_no1], center_coordinates);
348 center_point_xyz.X = center_coordinates[0];
349 center_point_xyz.Y = center_coordinates[1];
350 center_point_xyz.Z = center_coordinates[2];
351
352 set_center_point_plane_projection(coordinate_to_ignore, center_point_xyz, center_point_plane_projection);
353
354 auto winding_number
355 = winding_numbers_algorithm(cell_corners_plane_projection, ncorner + 1, center_point_plane_projection);
356 // printf("%d %g %g %g %g %d\n", cell_no2,
357 // RAD2DEG*grid_center_lon2[cell_no2],
358 // RAD2DEG*grid_center_lat2[cell_no2],
359 // RAD2DEG*grid_center_lon1[cell_no1],
360 // RAD2DEG*grid_center_lat1[cell_no1], winding_number);
361 if (winding_number != 0)
362 {
363 if (k >= MAX_CHILDS) cdo_abort("Internal problem, limit of MAX_CHILDS reached (limit=9).");
364 child2[cell_no2 * MAX_CHILDS + k++] = (long) cell_no1;
365 }
366 }
367 }
368 // printf("%d k = %d\n", cell_no2, k);
369 }
370
371 grid_point_search_delete(gps);
372 }
373
374 static void
compute_child_from_coordinates(CellIndex * cellindex1,CellIndex * cellindex2)375 compute_child_from_coordinates(CellIndex *cellindex1, CellIndex *cellindex2)
376 {
377 if (Options::cdoVerbose) cdo_print("%s", __func__);
378
379 const auto ncells1 = cellindex1->ncells;
380 const auto ncells2 = cellindex2->ncells;
381
382 Varray<double> lon1(ncells1), lat1(ncells1), lon2(ncells2), lat2(ncells2);
383 Varray<double> lon2_bnds(3 * ncells2), lat2_bnds(3 * ncells2);
384
385 read_coordinates(cellindex1->filename, ncells1, lon1.data(), lat1.data(), 0, nullptr, nullptr);
386 read_coordinates(cellindex2->filename, ncells2, lon2.data(), lat2.data(), 3, lon2_bnds.data(), lat2_bnds.data());
387
388 compute_child_from_bounds(cellindex2, ncells2, lon2.data(), lat2.data(), lon2_bnds.data(), lat2_bnds.data(), ncells1, lon1, lat1);
389 }
390
391 static void
compute_child(CellIndex * cellindex1,CellIndex * cellindex2)392 compute_child(CellIndex *cellindex1, CellIndex *cellindex2)
393 {
394 bool lparent = true;
395 const auto ncells1 = cellindex1->ncells;
396 long *parent1 = cellindex1->parent;
397 {
398 long i;
399 for (i = 0; i < ncells1; ++i)
400 if (parent1[i] >= 0) break;
401 if (i == ncells1) lparent = false;
402 }
403
404 if (lparent)
405 compute_child_from_parent(cellindex1, cellindex2);
406 else
407 {
408 compute_child_from_coordinates(cellindex1, cellindex2);
409 // cdo_abort("Missing parent index of %s!", cellindex1->filename);
410 }
411 }
412
413 static void
compute_sum(long i,long & n,double & sum,double & sumq,long kci,CellIndex ** cellindex,const Varray<double> & array)414 compute_sum(long i, long &n, double &sum, double &sumq, long kci, CellIndex **cellindex, const Varray<double> &array)
415 {
416 // printf("compute: i, kci %d %d\n", i, kci);
417 const auto ncells2 = cellindex[kci]->ncells;
418 if (i < 0 || i > ncells2) cdo_abort("Child grid cell index %ld out of bounds %ld!", i, ncells2);
419
420 for (int k = 0; k < MAX_CHILDS; ++k)
421 {
422 long index = cellindex[kci]->child[i * MAX_CHILDS + k];
423 if (index == -1) break;
424 if (kci == 1)
425 {
426 sum += array[index];
427 sumq += array[index] * array[index];
428 n += 1;
429 }
430 else
431 compute_sum(index, n, sum, sumq, kci - 1, cellindex, array);
432 }
433 }
434
435 static void
samplegrid(double missval,long nci,CellIndex ** cellindex,const Varray<double> & array1,Varray<double> & array2,Varray<double> & array3)436 samplegrid(double missval, long nci, CellIndex **cellindex, const Varray<double> &array1, Varray<double> &array2,
437 Varray<double> &array3)
438 {
439 static bool lstat = true;
440 long kci = nci - 1;
441 const auto ncells2 = cellindex[kci]->ncells;
442 long nx = 0;
443 double x = 0;
444 #ifdef _OPENMP
445 //#pragma omp parallel for default(none) shared(missval, ncells2, kci, cellindex, array1, array2, array3)
446 #endif
447 for (long i = 0; i < ncells2; ++i)
448 {
449 long n = 0;
450 double sum = 0, sumq = 0;
451 compute_sum(i, n, sum, sumq, kci, cellindex, array1);
452 array2[i] = n ? sum / n : missval; // mean
453 double var1 = (n * n > n) ? (sumq * n - sum * sum) / (n * n - n) : missval;
454 if (var1 < 0 && var1 > -1.e-5) var1 = 0;
455 array3[i] = var_to_std(var1, missval); // std1
456 if (lstat && n)
457 {
458 nx++;
459 x += n;
460 }
461 }
462 if (Options::cdoVerbose && lstat)
463 {
464 lstat = false;
465 cdo_print("Mean number of childs %g", nx ? x / nx : 0);
466 }
467 }
468
469 void *
Samplegridicon(void * process)470 Samplegridicon(void *process)
471 {
472 cdo_initialize(process);
473
474 cdo_operator_add("samplegridicon", 0, 0, "sample grids");
475
476 const auto nsamplegrids = cdo_operator_argc();
477 if (nsamplegrids < 2) cdo_abort("Parameter missing!");
478
479 std::vector<CellIndex *> cellindex(nsamplegrids);
480
481 for (int i = 0; i < nsamplegrids; ++i)
482 {
483 cellindex[i] = read_cellindex(cdo_operator_argv(i));
484 cellindex[i]->filename = cdo_operator_argv(i).c_str();
485 if (Options::cdoVerbose) cdo_print("Found %ld grid cells in %s", cellindex[i]->ncells, cellindex[i]->filename);
486 }
487
488 for (int i = 0; i < nsamplegrids - 1; ++i) compute_child(cellindex[i], cellindex[i + 1]);
489
490 const auto gridID2 = read_grid(cdo_operator_argv(nsamplegrids - 1).c_str());
491
492 const auto streamID1 = cdo_open_read(0);
493 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
494
495 long gridsize = vlistGridsizeMax(vlistID1);
496 if (Options::cdoVerbose) cdo_print("Source gridsize = %zu", gridsize);
497 if (gridsize != cellindex[0]->ncells)
498 cdo_abort("Gridsize (%ld) of input stream and first grid (%ld) differ!", gridsize, cellindex[0]->ncells);
499 if (vlistNumber(vlistID1) != CDI_REAL) gridsize *= 2;
500 Varray<double> array1(gridsize);
501
502 const auto vlistID2 = vlistDuplicate(vlistID1);
503 const auto vlistID3 = vlistDuplicate(vlistID1);
504
505 const auto taxisID1 = vlistInqTaxis(vlistID1);
506 const auto taxisID2 = taxisDuplicate(taxisID1);
507 const auto taxisID3 = taxisDuplicate(taxisID1);
508 vlistDefTaxis(vlistID2, taxisID2);
509 vlistDefTaxis(vlistID3, taxisID3);
510
511 const auto ngrids = vlistNgrids(vlistID1);
512 for (int index = 0; index < ngrids; ++index)
513 {
514 const auto gridID = vlistGrid(vlistID1, index);
515 const auto gridtype = gridInqType(gridID);
516 if (!(gridtype == GRID_UNSTRUCTURED && gridInqNvertex(gridID) == 3))
517 cdo_abort("Unsupported gridtype: %s with %d corners", gridNamePtr(gridtype), gridInqNvertex(gridID));
518
519 vlistChangeGridIndex(vlistID2, index, gridID2);
520 vlistChangeGridIndex(vlistID3, index, gridID2);
521 }
522
523 const auto streamID2 = cdo_open_write(1);
524 cdo_def_vlist(streamID2, vlistID2);
525
526 const auto streamID3 = cdo_open_write(2);
527 cdo_def_vlist(streamID3, vlistID3);
528
529 auto gridsize2 = gridInqSize(gridID2);
530 if (Options::cdoVerbose) cdo_print("Target gridsize = %ld", gridsize2);
531 if (vlistNumber(vlistID2) != CDI_REAL) gridsize2 *= 2;
532 Varray<double> array2(gridsize2), array3(gridsize2);
533
534 int tsID = 0;
535 while (true)
536 {
537 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
538 if (nrecs == 0) break;
539
540 cdo_taxis_copy_timestep(taxisID2, taxisID1);
541 cdo_taxis_copy_timestep(taxisID3, taxisID1);
542 cdo_def_timestep(streamID2, tsID);
543 cdo_def_timestep(streamID3, tsID);
544
545 for (int recID = 0; recID < nrecs; recID++)
546 {
547 int varID, levelID;
548 cdo_inq_record(streamID1, &varID, &levelID);
549 size_t nmiss;
550 cdo_read_record(streamID1, array1.data(), &nmiss);
551
552 const auto missval = vlistInqVarMissval(vlistID1, varID);
553
554 samplegrid(missval, nsamplegrids, &cellindex[0], array1, array2, array3);
555
556 nmiss = varray_num_mv(gridsize2, array2, missval);
557 cdo_def_record(streamID2, varID, levelID);
558 cdo_write_record(streamID2, array2.data(), nmiss);
559
560 nmiss = varray_num_mv(gridsize2, array3, missval);
561 cdo_def_record(streamID3, varID, levelID);
562 cdo_write_record(streamID3, array3.data(), nmiss);
563 }
564
565 tsID++;
566 }
567
568 cdo_stream_close(streamID3);
569 cdo_stream_close(streamID2);
570 cdo_stream_close(streamID1);
571
572 vlistDestroy(vlistID2);
573 gridDestroy(gridID2);
574
575 for (int i = 0; i < nsamplegrids; ++i) free_cellindex(cellindex[i]);
576
577 cdo_finish();
578
579 return nullptr;
580 }
581