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 /*
9 This module contains the following operators:
10
11 Selgridcell selgridcell Select grid cells by indices
12 */
13
14 #include <cdi.h>
15
16 #include "cdo_options.h"
17 #include "process_int.h"
18 #include "cdo_cdi_wrapper.h"
19 #include <mpim_grid.h>
20 #include "gridreference.h"
21 #include "util_files.h"
22 #include "param_conversion.h"
23
24 int gengridcell(const int gridID1, const size_t gridsize2, const std::vector<long> &cellidx);
25
26 static int
genindexgrid(int gridID1,const size_t gridsize2,const std::vector<long> & cellidx)27 genindexgrid(int gridID1, const size_t gridsize2, const std::vector<long> &cellidx)
28 {
29 const auto gridID0 = gridID1;
30 auto gridtype1 = gridInqType(gridID1);
31
32 if (gridtype1 == GRID_LONLAT || gridtype1 == GRID_GAUSSIAN || gridtype1 == GRID_PROJECTION)
33 {
34 gridID1 = gridToCurvilinear(gridID1, 0);
35 gridtype1 = GRID_CURVILINEAR;
36 }
37 else if (gridtype1 == GRID_UNSTRUCTURED && !gridHasCoordinates(gridID1))
38 {
39 const auto reference = dereferenceGrid(gridID1);
40 if (reference.isValid) gridID1 = reference.gridID;
41 if (reference.notFound) cdo_warning("Reference to source grid not found!");
42 }
43
44 int gridID2 = -1;
45 if (gridtype1 == GRID_UNSTRUCTURED || gridtype1 == GRID_CURVILINEAR)
46 gridID2 = gengridcell(gridID1, gridsize2, cellidx);
47 else if (gridtype1 == GRID_GENERIC && gridInqYsize(gridID1) == 0)
48 gridID2 = gengridcell(gridID1, gridsize2, cellidx);
49
50 if (gridID0 != gridID1) gridDestroy(gridID1);
51
52 return gridID2;
53 }
54
55 static void
selectIndex(double * array1,double * array2,long nind,const std::vector<long> & cellidx)56 selectIndex(double *array1, double *array2, long nind, const std::vector<long> &cellidx)
57 {
58 for (long i = 0; i < nind; ++i) array2[i] = array1[cellidx[i]];
59 }
60
61 void *
Selgridcell(void * process)62 Selgridcell(void *process)
63 {
64 int gridID1 = -1, gridID2;
65 int index, gridtype = -1;
66 struct sindex_t
67 {
68 int gridID1, gridID2;
69 };
70 std::vector<int> indarr;
71
72 cdo_initialize(process);
73
74 cdo_operator_add("selgridcell", 0, 0, "grid cell indices (1-N)");
75 const auto DELGRIDCELL = cdo_operator_add("delgridcell", 0, 0, "grid cell indices (1-N)");
76
77 operator_input_arg(cdo_operator_enter(0));
78
79 const auto operatorID = cdo_operator_id();
80
81 if (cdo_operator_argc() < 1) cdo_abort("Too few arguments!");
82
83 int nind = 0;
84 if (cdo_operator_argc() == 1)
85 {
86 bool maskfile = true, indexfile = false;
87 const char *filename = cdo_operator_argv(0).c_str();
88 if (strncmp(filename, "index=", 6) == 0)
89 {
90 filename += 6;
91 indexfile = true;
92 }
93 else if (strncmp(filename, "mask=", 5) == 0)
94 {
95 filename += 5;
96 maskfile = true;
97 }
98
99 if (FileUtils::file_exists(filename))
100 {
101 if (indexfile)
102 {
103 size_t cdo_read_index(const char *indexfile, std::vector<int> &indarr);
104 const auto n = cdo_read_index(filename, indarr);
105 nind = n;
106 if (nind == 0) cdo_abort("Index file %s generates no input!", cdo_operator_argv(0).c_str());
107 }
108 else if (maskfile)
109 {
110 size_t cdo_read_mask(const char *maskfile, std::vector<bool> &imask);
111 std::vector<bool> mask;
112 const auto n = cdo_read_mask(filename, mask);
113 nind = 0;
114 for (size_t i = 0; i < n; ++i)
115 if (mask[i]) nind++;
116 if (nind == 0) cdo_abort("Mask is empty!");
117
118 indarr.resize(nind);
119 nind = 0;
120 for (size_t i = 0; i < n; ++i)
121 if (mask[i]) indarr[nind++] = i;
122
123 if (nind == 0) cdo_abort("Mask file %s generates no input!", cdo_operator_argv(0).c_str());
124 }
125 }
126 }
127
128 if (nind == 0)
129 {
130 indarr = cdo_argv_to_int(cdo_get_oper_argv());
131 nind = indarr.size();
132
133 if (Options::cdoVerbose)
134 for (int i = 0; i < nind; i++) cdo_print("int %d = %d", i + 1, indarr[i]);
135
136 for (int i = 0; i < nind; i++) indarr[i] -= 1;
137 }
138
139 if (nind == 0) cdo_abort("Argument %s generates no input!", cdo_operator_argv(0).c_str());
140
141 auto indmin = indarr[0];
142 auto indmax = indarr[0];
143 for (int i = 1; i < nind; i++)
144 {
145 indmin = std::min(indmin, indarr[i]);
146 indmax = std::max(indmax, indarr[i]);
147 }
148
149 if (indmin < 0) cdo_abort("Index < 1 not allowed!");
150
151 const auto streamID1 = cdo_open_read(0);
152
153 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
154 const auto vlistID2 = vlistDuplicate(vlistID1);
155
156 const auto taxisID1 = vlistInqTaxis(vlistID1);
157 const auto taxisID2 = taxisDuplicate(taxisID1);
158 vlistDefTaxis(vlistID2, taxisID2);
159
160 const auto nvars = vlistNvars(vlistID1);
161 std::vector<bool> vars(nvars);
162 for (int varID = 0; varID < nvars; varID++) vars[varID] = false;
163
164 const auto ngrids = vlistNgrids(vlistID1);
165 std::vector<sindex_t> sindex(ngrids);
166
167 long ncells = nind;
168 std::vector<long> cellidx;
169 if (operatorID == DELGRIDCELL)
170 {
171 const auto gridsize = vlistGridsizeMax(vlistID1);
172 ncells = gridsize - nind;
173 cellidx.resize(gridsize);
174 for (size_t i = 0; i < gridsize; ++i) cellidx[i] = 1;
175 for (long i = 0; i < nind; ++i) cellidx[indarr[i]] = 0;
176 long j = 0;
177 for (size_t i = 0; i < gridsize; ++i)
178 if (cellidx[i] == 1) cellidx[j++] = i;
179 if (j != ncells) cdo_abort("Internal error; number of cells differ");
180 }
181 else
182 {
183 cellidx.resize(nind);
184 for (int i = 0; i < nind; ++i) cellidx[i] = indarr[i];
185 }
186
187 if (ncells == 0) cdo_abort("Mask is empty!");
188
189 for (index = 0; index < ngrids; index++)
190 {
191 gridID1 = vlistGrid(vlistID1, index);
192 gridtype = gridInqType(gridID1);
193
194 const auto gridsize = gridInqSize(gridID1);
195 if (gridsize == 1) continue;
196 if (indmax >= (int) gridsize)
197 {
198 cdo_warning("Max grid index is greater than grid size, skipped grid %d!", index + 1);
199 continue;
200 }
201
202 gridID2 = genindexgrid(gridID1, ncells, cellidx);
203
204 if (gridID2 == -1)
205 {
206 cdo_warning("Unsupported grid type >%s<, skipped grid %d!", gridNamePtr(gridtype), index + 1);
207 continue;
208 }
209
210 sindex[index].gridID1 = gridID1;
211 sindex[index].gridID2 = gridID2;
212
213 vlistChangeGridIndex(vlistID2, index, gridID2);
214
215 for (int varID = 0; varID < nvars; varID++)
216 if (gridID1 == vlistInqVarGrid(vlistID1, varID)) vars[varID] = true;
217 }
218
219 {
220 int varID;
221 for (varID = 0; varID < nvars; varID++)
222 if (vars[varID]) break;
223
224 if (varID >= nvars) cdo_abort("No variables selected!");
225 }
226
227 const auto streamID2 = cdo_open_write(1);
228 cdo_def_vlist(streamID2, vlistID2);
229
230 auto gridsize = vlistGridsizeMax(vlistID1);
231 if (vlistNumber(vlistID1) != CDI_REAL) gridsize *= 2;
232 Varray<double> array1(gridsize);
233
234 auto gridsizemax = vlistGridsizeMax(vlistID2);
235 if (vlistNumber(vlistID2) != CDI_REAL) gridsizemax *= 2;
236 Varray<double> array2(gridsizemax);
237
238 int tsID = 0;
239 while (true)
240 {
241 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
242 if (nrecs == 0) break;
243
244 cdo_taxis_copy_timestep(taxisID2, taxisID1);
245 cdo_def_timestep(streamID2, tsID);
246
247 for (int recID = 0; recID < nrecs; recID++)
248 {
249 int varID, levelID;
250 cdo_inq_record(streamID1, &varID, &levelID);
251 size_t nmiss;
252 cdo_read_record(streamID1, array1.data(), &nmiss);
253
254 cdo_def_record(streamID2, varID, levelID);
255
256 if (vars[varID])
257 {
258 gridID1 = vlistInqVarGrid(vlistID1, varID);
259
260 for (index = 0; index < ngrids; index++)
261 if (gridID1 == sindex[index].gridID1) break;
262
263 if (index == ngrids) cdo_abort("Internal problem, grid not found!");
264
265 const auto gridsize2 = gridInqSize(sindex[index].gridID2);
266
267 selectIndex(array1.data(), array2.data(), ncells, cellidx);
268
269 if (nmiss)
270 {
271 const auto missval = vlistInqVarMissval(vlistID2, varID);
272 nmiss = varray_num_mv(gridsize2, array2, missval);
273 }
274
275 cdo_write_record(streamID2, array2.data(), nmiss);
276 }
277 else
278 {
279 cdo_write_record(streamID2, array1.data(), nmiss);
280 }
281 }
282
283 tsID++;
284 }
285
286 cdo_stream_close(streamID2);
287 cdo_stream_close(streamID1);
288
289 vlistDestroy(vlistID2);
290
291 cdo_finish();
292
293 return nullptr;
294 }
295