bes  Updated for version 3.20.10
HDFEOS2ArraySwathGeoMultiDimMapField.cc
1 // Retrieves the latitude and longitude of the HDF-EOS2 Swath without using dimension maps
3 //
4 // Authors: MuQun Yang <myang6@hdfgroup.org>
5 // Copyright (c) 2010-2020 The HDF Group
7 //For the swath using the multiple dimension maps,
8 // Latitude/longitude will be interpolated accordingly.
9 
10 #ifdef USE_HDFEOS2_LIB
11 
12 #include "HDFEOS2ArraySwathGeoMultiDimMapField.h"
13 #include <iostream>
14 #include <sstream>
15 #include <cassert>
16 #include <libdap/debug.h>
17 #include <libdap/InternalErr.h>
18 #include "BESDebug.h"
19 #include "HDF4RequestHandler.h"
20 
21 using namespace std;
22 using namespace libdap;
23 #define SIGNED_BYTE_TO_INT32 1
24 
25 bool
26 HDFEOS2ArraySwathGeoMultiDimMapField::read ()
27 {
28 
29  BESDEBUG("h4","Coming to HDFEOS2ArraySwathGeoMultiDimMapField read "<<endl);
30 
31  if(length() == 0)
32  return true;
33 
34  if(rank !=2)
35  throw InternalErr (__FILE__, __LINE__, "The field rank must be 2 for swath multi-dimension map reading.");
36 
37 #if 0
38  string check_pass_fileid_key_str="H4.EnablePassFileID";
39  bool check_pass_fileid_key = false;
40  check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
41 #endif
42 
43  bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
44 
45  // Declare offset, count and step
46  vector<int>offset;
47  offset.resize(rank);
48  vector<int>count;
49  count.resize(rank);
50  vector<int>step;
51  step.resize(rank);
52 
53  // Obtain offset,step and count from the client expression constraint
54  int nelms = format_constraint (&offset[0], &step[0], &count[0]);
55 
56  // Just declare offset,count and step in the int32 type.
57  vector<int32>offset32;
58  offset32.resize(rank);
59  vector<int32>count32;
60  count32.resize(rank);
61  vector<int32>step32;
62  step32.resize(rank);
63 
64  // Just obtain the offset,count and step in the datatype of int32.
65  for (int i = 0; i < rank; i++) {
66  offset32[i] = (int32) offset[i];
67  count32[i] = (int32) count[i];
68  step32[i] = (int32) step[i];
69  }
70 
71  int32 (*openfunc) (char *, intn);
72  int32 (*attachfunc) (int32, char *);
73  intn (*detachfunc) (int32);
74  intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
75  intn (*readfieldfunc) (int32, char *, int32 *, int32 *, int32 *, void *);
76 
77  // Define function pointers to handle the swath
78  string datasetname;
79  openfunc = SWopen;
80  attachfunc = SWattach;
81  detachfunc = SWdetach;
82  fieldinfofunc = SWfieldinfo;
83  readfieldfunc = SWreadfield;
84  datasetname = swathname;
85 
86  // Declare dimension map size, offset and inc
87  vector<int>dm_dimsizes;
88  dm_dimsizes.resize(rank);
89  vector<int>dm_offsets;
90  dm_offsets.resize(rank);
91  vector<int>dm_incs;
92  dm_incs.resize(rank);
93  bool no_interpolation = true;
94 
95  if(dim0inc != 0 || dim1inc !=0 || dim0offset != 0 || dim1offset != 0) {
96  dm_dimsizes[0] = dim0size;
97  dm_dimsizes[1] = dim1size;
98  dm_offsets[0] = dim0offset;
99  dm_offsets[1] = dim1offset;
100  dm_incs[0] = dim0inc;
101  dm_incs[1] = dim1inc;
102  no_interpolation = false;
103  }
104 
105  // We may eventually combine the following code with other code, so
106  // we don't add many comments from here to the end of the file.
107  // The jira ticket about combining code is HFRHANDLER-166.
108  int32 sfid = -1;
109  int32 swathid = -1;
110 
111  if(false == check_pass_fileid_key) {
112  sfid = openfunc (const_cast < char *>(filename.c_str ()), DFACC_READ);
113  if (sfid < 0) {
114  ostringstream eherr;
115  eherr << "File " << filename.c_str () << " cannot be open.";
116  throw InternalErr (__FILE__, __LINE__, eherr.str ());
117  }
118  }
119  else
120  sfid = swathfd;
121 
122 #if 0
123 cerr<<"swath name is "<<datasetname <<endl;
124 cerr<<"swath rank is "<<rank <<endl;
125 cerr<<"swath field name is "<<fieldname <<endl;
126 cerr<<"swath field new name is "<<name() <<endl;
127 
128 cerr<<"datadim0size "<<dim0size <<endl;
129 cerr<<"datadim1size "<<dim1size <<endl;
130 cerr<<"datadim0offset "<<dim0offset <<endl;
131 cerr<<"datadim1offset "<<dim1offset <<endl;
132 cerr<<"datadim0inc "<<dim0inc <<endl;
133 cerr<<"datadim1inc "<<dim1inc <<endl;
134 #endif
135 
136  swathid = attachfunc (sfid, const_cast < char *>(datasetname.c_str ()));
137 
138  if (swathid < 0) {
139  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
140  ostringstream eherr;
141  eherr << "Swath " << datasetname.c_str () << " cannot be attached.";
142  throw InternalErr (__FILE__, __LINE__, eherr.str ());
143  }
144 
145 
146  int32 tmp_rank = -1;
147  vector<int32>tmp_dims;
148  tmp_dims.resize(rank);
149  char tmp_dimlist[1024];
150  int32 type =-1;
151  intn r = -1;
152  r = fieldinfofunc (swathid, const_cast < char *>(fieldname.c_str ()),
153  &tmp_rank, &tmp_dims[0], &type, tmp_dimlist);
154  if (r != 0) {
155  detachfunc (swathid);
156  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
157  ostringstream eherr;
158  eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
159  throw InternalErr (__FILE__, __LINE__, eherr.str ());
160  }
161 
162  vector<int32> newdims;
163  newdims.resize(tmp_rank);
164 
165 
166  switch (type) {
167  case DFNT_INT8:
168  {
169  if(true == no_interpolation) {
170  vector<int8>val;
171  val.resize(nelms);
172  r = readfieldfunc (swathid, const_cast < char *>(fieldname.c_str ()), &offset32[0], &step32[0], &count32[0], &val[0]);
173  if (r != 0) {
174  detachfunc (swathid);
175  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
176  ostringstream eherr;
177  eherr << "field " << fieldname.c_str () << "cannot be read.";
178  throw InternalErr (__FILE__, __LINE__, eherr.str ());
179  }
180 #ifndef SIGNED_BYTE_TO_INT32
181  set_value ((dods_byte *) &val[0], nelms);
182 #else
183  vector<int32>newval;
184  newval.resize(nelms);
185 
186  for (int counter = 0; counter < nelms; counter++)
187  newval[counter] = (int32) (val[counter]);
188  set_value ((dods_int32 *) &newval[0], nelms);
189 #endif
190  }
191  else {
192  // Obtaining the total value and interpolating the data
193  // according to dimension map
194  vector <int8> total_val8;
195  r = GetFieldValue (swathid, fieldname, dm_dimsizes,dm_offsets,dm_incs, total_val8, newdims);
196 
197  if (r != 0) {
198  detachfunc (swathid);
199  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
200  ostringstream eherr;
201  eherr << "field " << fieldname.c_str () << "cannot be read.";
202  throw InternalErr (__FILE__, __LINE__, eherr.str ());
203 
204  }
205  vector<int8>val8;
206  val8.resize(nelms);
207  Field2DSubset (&val8[0], newdims[0], newdims[1],&total_val8[0],
208  &offset32[0], &count32[0], &step32[0]);
209 
210 #ifndef SIGNED_BYTE_TO_INT32
211  set_value ((dods_byte *) &val8[0], nelms);
212 #else
213  vector<int32>newval;
214  newval.resize(nelms);
215 
216  for (int counter = 0; counter < nelms; counter++)
217  newval[counter] = (int32) (val8[counter]);
218  set_value ((dods_int32 *) &newval[0], nelms);
219 #endif
220  }
221  }
222  break;
223  case DFNT_UINT8:
224  case DFNT_UCHAR8:
225  {
226  if(no_interpolation == false) {
227  vector <uint8> total_val_uint8;
228  r = GetFieldValue (swathid, fieldname, dm_dimsizes,dm_offsets,dm_incs, total_val_uint8, newdims);
229 
230  if (r != 0) {
231  detachfunc (swathid);
232  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
233  ostringstream eherr;
234  eherr << "field " << fieldname.c_str () << "cannot be read.";
235  throw InternalErr (__FILE__, __LINE__, eherr.str ());
236 
237  }
238  vector<uint8>val_uint8;
239  val_uint8.resize(nelms);
240  Field2DSubset (&val_uint8[0], newdims[0], newdims[1],&total_val_uint8[0],
241  &offset32[0], &count32[0], &step32[0]);
242 
243  set_value ((dods_byte *) &val_uint8[0], nelms);
244  }
245  else {
246 
247  vector<uint8>val;
248  val.resize(nelms);
249  r = readfieldfunc (swathid, const_cast < char *>(fieldname.c_str ()), &offset32[0], &step32[0], &count32[0], &val[0]);
250  if (r != 0) {
251  detachfunc (swathid);
252  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
253  ostringstream eherr;
254  eherr << "field " << fieldname.c_str () << "cannot be read.";
255  throw InternalErr (__FILE__, __LINE__, eherr.str ());
256  }
257 
258  set_value ((dods_byte *) &val[0], nelms);
259  }
260  }
261  break;
262 
263  case DFNT_INT16:
264  {
265  if(no_interpolation == false) {
266  vector <int16> total_val_int16;
267  r = GetFieldValue (swathid, fieldname, dm_dimsizes,dm_offsets,dm_incs, total_val_int16, newdims);
268 
269  if (r != 0) {
270  detachfunc (swathid);
271  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
272  ostringstream eherr;
273  eherr << "field " << fieldname.c_str () << "cannot be read.";
274  throw InternalErr (__FILE__, __LINE__, eherr.str ());
275 
276  }
277  vector<int16>val_int16;
278  val_int16.resize(nelms);
279  Field2DSubset (&val_int16[0], newdims[0], newdims[1],&total_val_int16[0],
280  &offset32[0], &count32[0], &step32[0]);
281 
282  set_value ((dods_int16 *) &val_int16[0], nelms);
283  }
284 
285  else {
286  vector<int16>val;
287  val.resize(nelms);
288  r = readfieldfunc (swathid, const_cast < char *>(fieldname.c_str ()), &offset32[0], &step32[0], &count32[0], &val[0]);
289  if (r != 0) {
290  detachfunc (swathid);
291  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
292  ostringstream eherr;
293  eherr << "field " << fieldname.c_str () << "cannot be read.";
294  throw InternalErr (__FILE__, __LINE__, eherr.str ());
295  }
296  set_value ((dods_int16 *) &val[0], nelms);
297  }
298  }
299  break;
300 
301  case DFNT_UINT16:
302  {
303  if(no_interpolation == false) {
304  vector <uint16> total_val_uint16;
305  r = GetFieldValue (swathid, fieldname, dm_dimsizes,dm_offsets,dm_incs, total_val_uint16, newdims);
306 
307  if (r != 0) {
308  detachfunc (swathid);
309  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
310  ostringstream eherr;
311  eherr << "field " << fieldname.c_str () << "cannot be read.";
312  throw InternalErr (__FILE__, __LINE__, eherr.str ());
313 
314  }
315  vector<uint16>val_uint16;
316  val_uint16.resize(nelms);
317  Field2DSubset (&val_uint16[0], newdims[0], newdims[1],&total_val_uint16[0],
318  &offset32[0], &count32[0], &step32[0]);
319 
320  set_value ((dods_uint16 *) &val_uint16[0], nelms);
321  }
322 
323  else {
324  vector<uint16>val;
325  val.resize(nelms);
326  r = readfieldfunc (swathid, const_cast < char *>(fieldname.c_str ()), &offset32[0], &step32[0], &count32[0], &val[0]);
327  if (r != 0) {
328  detachfunc (swathid);
329  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
330  ostringstream eherr;
331  eherr << "field " << fieldname.c_str () << "cannot be read.";
332  throw InternalErr (__FILE__, __LINE__, eherr.str ());
333  }
334 
335  set_value ((dods_uint16 *) &val[0], nelms);
336  }
337  }
338  break;
339 
340  case DFNT_INT32:
341  {
342  if(no_interpolation == false) {
343  vector <int32> total_val_int32;
344  r = GetFieldValue (swathid, fieldname, dm_dimsizes,dm_offsets,dm_incs, total_val_int32, newdims);
345 
346  if (r != 0) {
347  detachfunc (swathid);
348  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
349  ostringstream eherr;
350  eherr << "field " << fieldname.c_str () << "cannot be read.";
351  throw InternalErr (__FILE__, __LINE__, eherr.str ());
352 
353  }
354  vector<int32>val_int32;
355  val_int32.resize(nelms);
356  Field2DSubset (&val_int32[0], newdims[0], newdims[1],&total_val_int32[0],
357  &offset32[0], &count32[0], &step32[0]);
358 
359  set_value ((dods_int32 *) &val_int32[0], nelms);
360  }
361  else {
362  vector<int32>val;
363  val.resize(nelms);
364  r = readfieldfunc (swathid, const_cast < char *>(fieldname.c_str ()), &offset32[0], &step32[0], &count32[0], &val[0]);
365  if (r != 0) {
366  detachfunc (swathid);
367  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
368  ostringstream eherr;
369  eherr << "field " << fieldname.c_str () << "cannot be read.";
370  throw InternalErr (__FILE__, __LINE__, eherr.str ());
371  }
372 
373  set_value ((dods_int32 *) &val[0], nelms);
374  }
375  }
376  break;
377 
378  case DFNT_UINT32:
379  {
380  if(no_interpolation == false) {
381  vector <uint32> total_val_uint32;
382  r = GetFieldValue (swathid, fieldname, dm_dimsizes,dm_offsets,dm_incs, total_val_uint32, newdims);
383 
384  if (r != 0) {
385  detachfunc (swathid);
386  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
387  ostringstream eherr;
388  eherr << "field " << fieldname.c_str () << "cannot be read.";
389  throw InternalErr (__FILE__, __LINE__, eherr.str ());
390 
391  }
392  vector<uint32>val_uint32;
393  val_uint32.resize(nelms);
394  Field2DSubset (&val_uint32[0], newdims[0], newdims[1],&total_val_uint32[0],
395  &offset32[0], &count32[0], &step32[0]);
396 
397  set_value ((dods_uint32 *) &val_uint32[0], nelms);
398  }
399  else {
400  vector<uint32>val;
401  val.resize(nelms);
402  r = readfieldfunc (swathid, const_cast < char *>(fieldname.c_str ()), &offset32[0], &step32[0], &count32[0], &val[0]);
403  if (r != 0) {
404  detachfunc (swathid);
405  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
406  ostringstream eherr;
407  eherr << "field " << fieldname.c_str () << "cannot be read.";
408  throw InternalErr (__FILE__, __LINE__, eherr.str ());
409  }
410 
411  set_value ((dods_uint32 *) &val[0], nelms);
412  }
413  }
414  break;
415  case DFNT_FLOAT32:
416  {
417  if(no_interpolation == false) {
418  vector <float32> total_val_f32;
419  r = GetFieldValue (swathid, fieldname, dm_dimsizes,dm_offsets,dm_incs, total_val_f32, newdims);
420 
421  if (r != 0) {
422  detachfunc (swathid);
423  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
424  ostringstream eherr;
425  eherr << "field " << fieldname.c_str () << "cannot be read.";
426  throw InternalErr (__FILE__, __LINE__, eherr.str ());
427 
428  }
429  vector<float32>val_f32;
430  val_f32.resize(nelms);
431  Field2DSubset (&val_f32[0], newdims[0], newdims[1],&total_val_f32[0],
432  &offset32[0], &count32[0], &step32[0]);
433 
434  set_value ((dods_float32 *) &val_f32[0], nelms);
435  }
436  else {
437  vector<float32>val;
438  val.resize(nelms);
439  r = readfieldfunc (swathid, const_cast < char *>(fieldname.c_str ()), &offset32[0], &step32[0], &count32[0], &val[0]);
440  if (r != 0) {
441  detachfunc (swathid);
442  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
443  ostringstream eherr;
444  eherr << "field " << fieldname.c_str () << "cannot be read.";
445  throw InternalErr (__FILE__, __LINE__, eherr.str ());
446  }
447  set_value ((dods_float32 *) &val[0], nelms);
448  }
449  }
450  break;
451  case DFNT_FLOAT64:
452  {
453  if(no_interpolation == false) {
454  vector <float64> total_val_f64;
455  r = GetFieldValue (swathid, fieldname, dm_dimsizes,dm_offsets,dm_incs, total_val_f64, newdims);
456 
457  if (r != 0) {
458  detachfunc (swathid);
459  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
460  ostringstream eherr;
461  eherr << "field " << fieldname.c_str () << "cannot be read.";
462  throw InternalErr (__FILE__, __LINE__, eherr.str ());
463 
464  }
465  vector<float64>val_f64;
466  val_f64.resize(nelms);
467  Field2DSubset (&val_f64[0], newdims[0], newdims[1],&total_val_f64[0],
468  &offset32[0], &count32[0], &step32[0]);
469 
470  set_value ((dods_float64 *) &val_f64[0], nelms);
471  }
472  else {
473  vector<float64>val;
474  val.resize(nelms);
475  r = readfieldfunc (swathid, const_cast < char *>(fieldname.c_str ()), &offset32[0], &step32[0], &count32[0], &val[0]);
476  if (r != 0) {
477  detachfunc (swathid);
478  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
479  ostringstream eherr;
480  eherr << "field " << fieldname.c_str () << "cannot be read.";
481  throw InternalErr (__FILE__, __LINE__, eherr.str ());
482  }
483 
484  set_value ((dods_float64 *) &val[0], nelms);
485  }
486  }
487  break;
488  default:
489  detachfunc (swathid);
490  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
491  throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
492  }
493 
494  r = detachfunc (swathid);
495  if (r != 0) {
496  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
497  ostringstream eherr;
498  eherr << "Swath " << datasetname.c_str () << " cannot be detached.";
499  throw InternalErr (__FILE__, __LINE__, eherr.str ());
500  }
501 
502  HDFCFUtil::close_fileid(-1,-1,-1,sfid,check_pass_fileid_key);
503 
504 #if 0
505  r = closefunc (sfid);
506  if (r != 0) {
507  ostringstream eherr;
508  eherr << "Swath " << filename.c_str () << " cannot be closed.";
509  throw InternalErr (__FILE__, __LINE__, eherr.str ());
510  }
511 #endif
512 
513  return false;
514 }
515 
516 // Standard way of DAP handlers to pass the coordinates of the subsetted region to the handlers
517 // Return the number of elements to read.
518 int
519 HDFEOS2ArraySwathGeoMultiDimMapField::format_constraint (int *offset, int *step, int *count)
520 {
521  long nels = 1;
522  int id = 0;
523 
524  Dim_iter p = dim_begin ();
525  while (p != dim_end ()) {
526 
527  int start = dimension_start (p, true);
528  int stride = dimension_stride (p, true);
529  int stop = dimension_stop (p, true);
530 
531  // Check for illegal constraint
532  if (start > stop) {
533  ostringstream oss;
534  oss << "Array/Grid hyperslab start point "<< start <<
535  " is greater than stop point " << stop <<".";
536  throw Error(malformed_expr, oss.str());
537  }
538 
539  offset[id] = start;
540  step[id] = stride;
541  count[id] = ((stop - start) / stride) + 1; // count of elements
542  nels *= count[id]; // total number of values for variable
543 
544  BESDEBUG ("h4",
545  "=format_constraint():"
546  << "id=" << id << " offset=" << offset[id]
547  << " step=" << step[id]
548  << " count=" << count[id]
549  << endl);
550 
551  id++;
552  p++;
553  }// while (p != dim_end ())
554 
555  return nels;
556 }
557 
558 // Subset of latitude and longitude to follow the parameters
559 // from the DAP expression constraint
560 template < class T >
561 bool HDFEOS2ArraySwathGeoMultiDimMapField::Field2DSubset (T * outlatlon,
562  const int /*majordim //unused SBL 2/7/20 */,
563  const int minordim,
564  T * latlon,
565  int32 * offset,
566  int32 * count,
567  int32 * step)
568 {
569 #if 0
570  T (*templatlonptr)[majordim][minordim] = (T *[][]) latlon;
571 #endif
572  int i = 0;
573  int j = 0;
574 
575  // do subsetting
576  // Find the correct index
577  int dim0count = count[0];
578  int dim1count = count[1];
579 
580  int dim0index[dim0count];
581  int dim1index[dim1count];
582 
583  for (i = 0; i < count[0]; i++) // count[0] is the least changing dimension
584  dim0index[i] = offset[0] + i * step[0];
585 
586 
587  for (j = 0; j < count[1]; j++)
588  dim1index[j] = offset[1] + j * step[1];
589 
590  // Now assign the subsetting data
591  int k = 0;
592 
593  for (i = 0; i < count[0]; i++) {
594  for (j = 0; j < count[1]; j++) {
595 #if 0
596  outlatlon[k] = (*templatlonptr)[dim0index[i]][dim1index[j]];
597 #endif
598  outlatlon[k] = *(latlon + (dim0index[i] * minordim) + dim1index[j]);
599  k++;
600  }
601  }
602  return true;
603 }
604 
605 // Get latitude and longitude fields.
606 // It will call expand_dimmap_field to interpolate latitude and longitude.
607 template < class T > int
608 HDFEOS2ArraySwathGeoMultiDimMapField::
609 GetFieldValue (int32 swathid, const string & geofieldname,
610  const vector <int>&dimsizes,const vector<int>&offset,const vector<int>&inc,
611  vector < T > &vals, vector<int32>&newdims)
612 {
613 
614  int32 ret = -1;
615  int32 size = -1;
616  int32 sw_rank = -1;
617  int32 dims[130];
618  int32 type = -1;
619 
620  // Two dimensions for lat/lon; each dimension name is < 64 characters,
621  // The dimension names are separated by a comma.
622  char dimlist[130];
623  ret = SWfieldinfo (swathid, const_cast < char *>(geofieldname.c_str ()),
624  &sw_rank, dims, &type, dimlist);
625  if (ret != 0)
626  return -1;
627 
628  if(sw_rank !=2)
629  return -1;
630 
631  size = 1;
632  for (int i = 0; i <sw_rank; i++)
633  size *= dims[i];
634 
635  vals.resize (size);
636 
637  ret = SWreadfield (swathid, const_cast < char *>(geofieldname.c_str ()),
638  NULL, NULL, NULL, (void *) &vals[0]);
639  if (ret != 0)
640  return -1;
641 
642  vector < string > dimname;
643  HDFCFUtil::Split (dimlist, ',', dimname);
644 
645  for (int i = 0; i < sw_rank; i++) {
646 
647 #if 0
648 cerr<<"dimsizes["<<i<<"]: " <<dimsizes[i]<<endl;
649 cerr<<"offset["<<i<<"]: " <<offset[i]<<endl;
650 cerr<<"inc["<<i<<"]: " <<inc[i]<<endl;
651 #endif
652  int r;
653  r = _expand_dimmap_field (&vals, sw_rank, dims, i, dimsizes[i], offset[i], inc[i]);
654  if (r != 0)
655  return -1;
656  }
657 
658  // dims[] are expanded already.
659  for (int i = 0; i < sw_rank; i++) {
660  //cerr<<"i "<< i << " "<< dims[i] <<endl;
661  if (dims[i] < 0)
662  return -1;
663  newdims[i] = dims[i];
664  }
665 
666  return 0;
667 }
668 
669 // expand the dimension map field.
670 template < class T > int
671 HDFEOS2ArraySwathGeoMultiDimMapField::_expand_dimmap_field (vector < T >
672  *pvals, int32 sw_rank,
673  int32 dimsa[],
674  int dimindex,
675  int32 ddimsize,
676  int32 offset,
677  int32 inc)
678 {
679  vector < T > orig = *pvals;
680  vector < int32 > pos;
681  vector < int32 > dims;
682  vector < int32 > newdims;
683  pos.resize (sw_rank);
684  dims.resize (sw_rank);
685 
686  for (int i = 0; i < sw_rank; i++) {
687  pos[i] = 0;
688  dims[i] = dimsa[i];
689  }
690  newdims = dims;
691  newdims[dimindex] = ddimsize;
692  dimsa[dimindex] = ddimsize;
693 
694  int newsize = 1;
695 
696  for (int i = 0; i < sw_rank; i++) {
697  newsize *= newdims[i];
698  }
699  pvals->clear ();
700  pvals->resize (newsize);
701 
702  for (;;) {
703  // if end
704  if (pos[0] == dims[0]) {
705  // we past then end
706  break;
707  }
708  else if (pos[dimindex] == 0) {
709  // extract 1D values
710  vector < T > v;
711  for (int i = 0; i < dims[dimindex]; i++) {
712  pos[dimindex] = i;
713  v.push_back (orig[INDEX_nD_TO_1D (dims, pos)]);
714  }
715  // expand them
716 
717  vector < T > w;
718  for (int32 j = 0; j < ddimsize; j++) {
719  int32 i = (j - offset) / inc;
720  T f;
721 
722  if (i * inc + offset == j) // perfect match
723  {
724  f = (v[i]);
725  }
726  else {
727  int32 i1 = 0;
728  int32 i2 = (i<=0)?1:0;
729  int32 j1 = 0;
730  int32 j2 = 0;
731 
732 #if 0
733  if (i <= 0) {
734  //i1 = 0;
735  i2 = 1;
736  }
737 #endif
738  if ((unsigned int) i + 1 >= v.size ()) {
739  i1 = v.size () - 2;
740  i2 = v.size () - 1;
741  }
742  else {
743  i1 = i;
744  i2 = i + 1;
745  }
746  j1 = i1 * inc + offset;
747  j2 = i2 * inc + offset;
748  f = (((j - j1) * v[i2] + (j2 - j) * v[i1]) / (j2 - j1));
749  }
750  w.push_back (f);
751  pos[dimindex] = j;
752  (*pvals)[INDEX_nD_TO_1D (newdims, pos)] = f;
753  }
754  pos[dimindex] = 0;
755  }
756  // next pos
757  pos[sw_rank - 1]++;
758  for (int i = sw_rank - 1; i > 0; i--) {
759  if (pos[i] == dims[i]) {
760  pos[i] = 0;
761  pos[i - 1]++;
762  }
763  }
764  }
765 
766  return 0;
767 }
768 
769 
770 #endif
static void Split(const char *s, int len, char sep, std::vector< std::string > &names)
Definition: HDFCFUtil.cc:80
static void close_fileid(int32 sdfd, int32 file_id, int32 gridfd, int32 swathfd, bool pass_fileid_key)
Definition: HDFCFUtil.cc:3669