bes Updated for version 3.21.1
The Backend Server (BES) is the lower two tiers of the Hyrax data server
HDF5MissLLArray.cc
1
2// This file is part of the hdf5 data handler for the OPeNDAP data server.
3// Currently, it provides the missing latitude,longitude fields for the HDF-EOS5 files for the default option.
5
6#include <sys/stat.h>
7#include <iostream>
8#include <sstream>
9#include <memory>
10#include <BESDebug.h>
11#include <libdap/InternalErr.h>
12
13#include "HDF5MissLLArray.h"
14#include "HDF5CFUtil.h"
15#if 0
16#include "HDF5RequestHandler.h"
17#endif
18
19using namespace std;
20using namespace libdap;
21
22BaseType *HDF5MissLLArray::ptr_duplicate()
23{
24 auto HDF5MissLLArray_unique = make_unique<HDF5MissLLArray>(*this);
25 return HDF5MissLLArray_unique.release();
26}
27
28bool HDF5MissLLArray::read()
29{
30
31 BESDEBUG("h5","Coming to HDF5MissLLArray read "<<endl);
32
33 if (g_info.projection == HE5_GCTP_GEO)
34 read_data_geo();
35 else
36 read_data_non_geo();
37
38 return true;
39}
40
41
42bool HDF5MissLLArray::read_data_non_geo() {
43
44 int64_t nelms = -1;
45 vector<int64_t>offset;
46 vector<int64_t>count;
47 vector<int64_t>step;
48
49 if (rank <= 0)
50 throw InternalErr (__FILE__, __LINE__,
51 "The number of dimension of this variable should be greater than 0");
52 else {
53 offset.resize(rank);
54 count.resize(rank);
55 step.resize(rank);
56 nelms = format_constraint (offset.data(), step.data(), count.data());
57 }
58
59 if (nelms <= 0)
60 throw InternalErr (__FILE__, __LINE__,
61 "The number of elments is negative.");
62
63 int64_t total_elms = g_info.xdim_size * g_info.ydim_size;
64 if (total_elms > DODS_INT_MAX)
65 throw InternalErr (__FILE__, __LINE__,
66 "Currently we cannot calculate lat/lon that is greater than 2G for HDF-EOS5.");
67
68 if(g_info.ydim_size <=0 || g_info.xdim_size <=0)
69 throw InternalErr (__FILE__, __LINE__,
70 "The number of elments at each dimension should be greater than 0.");
71
72 vector<size_t>pos(rank,0);
73 for (int i = 0; i< rank; i++)
74 pos[i] = offset[i];
75
76 vector<size_t>dimsizes;
77 dimsizes.push_back(g_info.ydim_size);
78 dimsizes.push_back(g_info.xdim_size);
79
80 double upleft[2];
81 double lowright[2];
82 vector<int>rows;
83 vector<int>cols;
84 vector<double>lon;
85 vector<double>lat;
86 rows.resize(g_info.xdim_size * g_info.ydim_size);
87 cols.resize(g_info.xdim_size * g_info.ydim_size);
88 lon.resize(g_info.xdim_size * g_info.ydim_size);
89 lat.resize(g_info.xdim_size * g_info.ydim_size);
90
91 upleft[0] = g_info.point_left;
92 upleft[1] = g_info.point_upper;
93 lowright[0] = g_info.point_right;
94 lowright[1] = g_info.point_lower;
95
96 int r = -1;
97
98 int k = 0;
99 for (int j = 0; j < g_info.ydim_size; ++j) {
100 for (int i = 0; i < g_info.xdim_size; ++i) {
101 rows[k] = j;
102 cols[k] = i;
103 ++k;
104 }
105 }
106
107 BESDEBUG("h5", " Before calling GDij2ll, check all projection parameters. " << endl);
108 BESDEBUG("h5", " eos5_projcode is " << g_info.projection <<endl);
109 BESDEBUG("h5", " eos5_zone is " << g_info.zone <<endl);
110 BESDEBUG("h5", " eos5_params[0] is " << g_info.param[0] <<endl);
111 BESDEBUG("h5", " eos5_params[1] is " << g_info.param[1] <<endl);
112 BESDEBUG("h5", " eos5_sphere is " << g_info.sphere <<endl);
113 BESDEBUG("h5", " xdimsize is " << g_info.xdim_size <<endl);
114 BESDEBUG("h5", " ydimsize is " << g_info.ydim_size <<endl);
115 BESDEBUG("h5", " eos5_pixelreg is " << g_info.pixelregistration <<endl);
116 BESDEBUG("h5", " eos5_origin is " << g_info.gridorigin <<endl);
117 BESDEBUG("h5", " upleft[0] is " << upleft[0] <<endl);
118 BESDEBUG("h5", " upleft[1] is " << upleft[1] <<endl);
119 BESDEBUG("h5", " lowright[0] is " << lowright[0] <<endl);
120 BESDEBUG("h5", " lowright[1] is " << lowright[1] <<endl);
121
122
123 // Calculate Lat/lon by using GCTP
124 r = GDij2ll (g_info.projection, g_info.zone, g_info.param,g_info.sphere, g_info.xdim_size, g_info.ydim_size, upleft, lowright,
125 g_info.xdim_size * g_info.ydim_size, rows.data(), cols.data(), lon.data(), lat.data(), g_info.pixelregistration, g_info.gridorigin);
126 if (r != 0) {
127 ostringstream eherr;
128 eherr << "cannot calculate grid latitude and longitude";
129 throw InternalErr (__FILE__, __LINE__, eherr.str ());
130 }
131
132
133 BESDEBUG("h5", " The first value of lon is " << lon[0] <<endl);
134 BESDEBUG("h5", " The first value of lat is " << lat[0] <<endl);
135
136
137 if(is_lat) {
138 if(total_elms == nelms)
139 set_value_ll(lat.data(),total_elms);
140 else {
141 vector<double>val;
142 subset<double>(
143 lat.data(),
144 rank,
145 dimsizes,
146 offset.data(),
147 step.data(),
148 count.data(),
149 &val,
150 pos,
151 0);
152 set_value_ll(val.data(),nelms);
153 }
154 }
155 else {
156
157 if(total_elms == nelms)
158 set_value_ll(lon.data(),total_elms);
159 else {
160 vector<double>val;
161 subset<double>(
162 lon.data(),
163 rank,
164 dimsizes,
165 offset.data(),
166 step.data(),
167 count.data(),
168 &val,
169 pos,
170 0);
171 set_value_ll(val.data(),nelms);
172 }
173 }
174 return true;
175
176}
177
178bool HDF5MissLLArray::read_data_geo(){
179
180 BESDEBUG("h5","HDF5MissLLArray: Coming to read_data_geo "<<endl);
181 int64_t nelms = -1;
182 vector<int64_t>offset;
183 vector<int64_t>count;
184 vector<int64_t>step;
185
186
187 if (rank <= 0)
188 throw InternalErr (__FILE__, __LINE__,
189 "The number of dimension of this variable should be greater than 0");
190
191 offset.resize(rank);
192 count.resize(rank);
193 step.resize(rank);
194 nelms = format_constraint (offset.data(), step.data(), count.data());
195
196
197 if (nelms <= 0 || nelms >DODS_INT_MAX)
198 throw InternalErr (__FILE__, __LINE__,
199 "The number of elements for geographic lat/lon is negative or greater than 2G.");
200
201 vector<float>val;
202 val.resize(nelms);
203
204#if 0
205 if (is_lat) {
206
207 if (HE5_HDFE_GD_UL == g_info.gridorigin || HE5_HDFE_GD_UR == g_info.gridorigin) {
208
209 start = (float)(g_info.point_upper);
210 end = (float)(g_info.point_lower);
211
212 }
213 else {// (gridorigin == HE5_HDFE_GD_LL || gridorigin == HE5_HDFE_GD_LR)
214
215 start = (float)(g_info.point_lower);
216 end = (float)(g_info.point_upper);
217 }
218
219 if(g_info.ydim_size <=0)
220 throw InternalErr (__FILE__, __LINE__,
221 "The number of elments should be greater than 0.");
222
223 float lat_step = (end - start) /(float)(g_info.ydim_size);
224
225 if ( HE5_HDFE_CENTER == g_info.pixelregistration ) {
226 for (int i = 0; i < nelms; i++)
227 val[i] = (((float)(offset[0]+i*step[0]) + 0.5F) * lat_step + start) / 1000000.0F;
228 }
229 else { // HE5_HDFE_CORNER
230 for (int i = 0; i < nelms; i++)
231 val[i] = ((float)(offset[0]+i * step[0])*lat_step + start) / 1000000.0F;
232 }
233 }
234 else {
235
236 if (HE5_HDFE_GD_UL == g_info.gridorigin || HE5_HDFE_GD_LL == g_info.gridorigin) {
237
238 start = (float)(g_info.point_left);
239 end = (float)(g_info.point_right);
240
241 }
242 else {// (gridorigin == HE5_HDFE_GD_UR || gridorigin == HE5_HDFE_GD_LR)
243
244 start = (float)(g_info.point_right);
245 end = (float)(g_info.point_left);
246 }
247
248 if(g_info.xdim_size <=0)
249 throw InternalErr (__FILE__, __LINE__,
250 "The number of elments should be greater than 0.");
251 float lon_step = (end - start) /(float)(g_info.xdim_size);
252
253 if (HE5_HDFE_CENTER == g_info.pixelregistration) {
254
255 for (int i = 0; i < nelms; i++)
256 val[i] = (((float)(offset[0] + i *step[0]) + 0.5F) * lon_step + start ) / 1000000.0F;
257
258 }
259 else { // HE5_HDFE_CORNER
260 for (int i = 0; i < nelms; i++)
261 val[i] = ((float)(offset[0]+i*step[0]) * lon_step + start) / 1000000.0F;
262 }
263 }
264#endif
265#if 0
266for (int i =0; i <nelms; i++)
267"h5","final data val "<< i <<" is " << val[i] <<endl;
268#endif
269
270 if (is_lat)
271 read_data_geo_lat(nelms, offset, step, val) ;
272 else
273 read_data_geo_lon(nelms, offset, step, val) ;
274 set_value_ll(val.data(), nelms);
275
276
277 return true;
278}
279
280void HDF5MissLLArray::read_data_geo_lat(int64_t nelms, const vector<int64_t> &offset,
281 const vector<int64_t> &step, vector<float> &val) const
282{
283
284 float start = 0.0;
285 float end = 0.0;
286
287 if (HE5_HDFE_GD_UL == g_info.gridorigin || HE5_HDFE_GD_UR == g_info.gridorigin) {
288
289 start = (float)(g_info.point_upper);
290 end = (float)(g_info.point_lower);
291
292 }
293 else {// (gridorigin == HE5_HDFE_GD_LL || gridorigin == HE5_HDFE_GD_LR)
294
295 start = (float)(g_info.point_lower);
296 end = (float)(g_info.point_upper);
297 }
298
299 if(g_info.ydim_size <=0)
300 throw InternalErr (__FILE__, __LINE__,
301 "The number of elements should be greater than 0.");
302
303 float lat_step = (end - start) /(float)(g_info.ydim_size);
304
305 if ( HE5_HDFE_CENTER == g_info.pixelregistration ) {
306 for (int i = 0; i < nelms; i++)
307 val[i] = (((float)(offset[0]+i*step[0]) + 0.5F) * lat_step + start) / 1000000.0F;
308 }
309 else { // HE5_HDFE_CORNER
310 for (int i = 0; i < nelms; i++)
311 val[i] = ((float)(offset[0]+i * step[0])*lat_step + start) / 1000000.0F;
312 }
313
314}
315
316void HDF5MissLLArray::read_data_geo_lon(int64_t nelms, const vector<int64_t> &offset, const vector<int64_t> &step,
317 vector<float> &val) const {
318
319 float start = 0.0;
320 float end = 0.0;
321
322 if (HE5_HDFE_GD_UL == g_info.gridorigin || HE5_HDFE_GD_LL == g_info.gridorigin) {
323
324 start = (float)(g_info.point_left);
325 end = (float)(g_info.point_right);
326
327 }
328 else {// (gridorigin == HE5_HDFE_GD_UR || gridorigin == HE5_HDFE_GD_LR)
329
330 start = (float)(g_info.point_right);
331 end = (float)(g_info.point_left);
332 }
333
334 if (g_info.xdim_size <=0)
335 throw InternalErr (__FILE__, __LINE__,
336 "The number of elements should be greater than 0.");
337
338 float lon_step = (end - start) /(float)(g_info.xdim_size);
339
340 if (HE5_HDFE_CENTER == g_info.pixelregistration) {
341
342 for (int i = 0; i < nelms; i++)
343 val[i] = (((float)(offset[0] + i *step[0]) + 0.5F) * lon_step + start ) / 1000000.0F;
344
345 }
346 else { // HE5_HDFE_CORNER
347 for (int i = 0; i < nelms; i++)
348 val[i] = ((float)(offset[0]+i*step[0]) * lon_step + start) / 1000000.0F;
349 }
350
351}
352
353// parse constraint expr. and make hdf5 coordinate point location.
354// return number of elements to read.
355int64_t
356HDF5MissLLArray::format_constraint (int64_t *offset, int64_t *step, int64_t *count)
357{
358 int64_t nels = 1;
359 int id = 0;
360
361 Dim_iter p = dim_begin ();
362
363 while (p != dim_end ()) {
364
365 int64_t start = dimension_start (p, true);
366 int64_t stride = dimension_stride (p, true);
367 int64_t stop = dimension_stop (p, true);
368
369 // Check for illegal constraint
370 if (start > stop) {
371 ostringstream oss;
372 oss << "Array/Grid hyperslab start point "<< start <<
373 " is greater than stop point " << stop <<".";
374 throw Error(malformed_expr, oss.str());
375 }
376
377 offset[id] = start;
378 step[id] = stride;
379 count[id] = ((stop - start) / stride) + 1; // count of elements
380 nels *= count[id]; // total number of values for variable
381
382 BESDEBUG ("h5",
383 "=format_constraint():"
384 << "id=" << id << " offset=" << offset[id]
385 << " step=" << step[id]
386 << " count=" << count[id]
387 << endl);
388
389 id++;
390 p++;
391 }
392
393 return nels;
394}
395
397//
398// \param input Input variable
399// \param dim dimension info of the input
400// \param start start indexes of each dim
401// \param stride stride of each dim
402// \param edge count of each dim
403// \param poutput output variable
404// \parrm index dimension index
405// \return 0 if successful. -1 otherwise.
406//
407template<typename T>
408int HDF5MissLLArray::subset(
409 void* input,
410 int s_rank,
411 const vector<size_t> & dim,
412 int64_t start[],
413 int64_t stride[],
414 int64_t edge[],
415 vector<T> *poutput,
416 vector<size_t>& pos,
417 int index)
418{
419 for(int k=0; k<edge[index]; k++)
420 {
421 pos[index] = start[index] + k*stride[index];
422 if(index+1<s_rank)
423 subset(input, s_rank, dim, start, stride, edge, poutput,pos,index+1);
424 if(index==s_rank-1)
425 {
426 size_t cur_pos = INDEX_nD_TO_1D( dim, pos);
427 auto tempbuf = (void*)((char*)input+cur_pos*sizeof(T));
428 poutput->push_back(*(static_cast<T*>(tempbuf)));
429 }
430 } // end of for
431 return 0;
432} // end of template<typename T> static int subset
433
434size_t HDF5MissLLArray::INDEX_nD_TO_1D (const std::vector < size_t > &dims,
435 const std::vector < size_t > &pos) const {
436 //
437 // "int a[10][20][30] // & a[1][2][3] == a + (20*30+1 + 30*2 + 1 *3)"
438 // "int b[10][2] // &b[1][1] == b + (2*1 + 1)"
439 //
440 if(dims.size () != pos.size ())
441 throw InternalErr(__FILE__,__LINE__,"dimension error in INDEX_nD_TO_1D routine.");
442 size_t sum = 0;
443 size_t start = 1;
444
445 for (const auto &p:pos) {
446 size_t m = 1;
447
448 for (size_t j = start; j < dims.size (); j++)
449 m *= dims[j];
450 sum += m *p;
451 start++;
452 }
453 return sum;
454}
455
HDF5CFUtil.h
This file includes several helper functions for translating HDF5 to CF-compliant.
HDF5RequestHandler.h
include the entry functions to execute the handlers