bes Updated for version 3.21.1
The Backend Server (BES) is the lower two tiers of the Hyrax data server
HDFEOS2ArrayGridGeoField.cc
1
2// retrieves the latitude and longitude of the HDF-EOS2 Grid
3// Authors: Kent Yang <myang6@hdfgroup.org>
4// Copyright (c) The HDF Group
6#ifdef USE_HDFEOS2_LIB
7
8#include "HDFEOS2ArrayGridGeoField.h"
9
10#include <stdio.h>
11#include <stdlib.h>
12#include <sys/stat.h>
13#include <iostream>
14#include <sstream>
15#include <cassert>
16#include <libdap/debug.h>
17#include "HDFEOS2.h"
18#include <libdap/InternalErr.h>
19#include <BESDebug.h>
20#include "HDFCFUtil.h"
21
22#include "misrproj.h"
23#include "errormacros.h"
24#include <proj.h>
25#include <sys/types.h>
26#include <fcntl.h>
27#include <unistd.h>
28
29#include "BESH4MCache.h"
30#include "HDF4RequestHandler.h"
31
32using namespace std;
33using namespace libdap;
34
35#define SIGNED_BYTE_TO_INT32 1
36
37// These two functions are used to handle MISR products with the SOM projections.
38extern "C" {
39 int inv_init(int insys, int inzone, double *inparm, int indatum, char *fn27, char *fn83, int *iflg, int (*inv_trans[])(double, double, double*, double*));
40 int sominv(double y, double x, double *lon, double *lat);
41}
42
43bool
44HDFEOS2ArrayGridGeoField::read ()
45{
46 BESDEBUG("h4","Coming to HDFEOS2ArrayGridGeoField read "<<endl);
47 if (length() == 0)
48 return true;
49
50 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
51
52 // Currently The latitude and longitude rank from HDF-EOS2 grid must be either 1-D or 2-D.
53 // However, For SOM projection the final rank will become 3.
54 if (rank < 1 || rank > 2) {
55 throw InternalErr (__FILE__, __LINE__, "The rank of geo field is greater than 2, currently we don't support 3-D lat/lon cases.");
56 }
57
58 // MISR SOM file's final rank is 3. So declare a new variable.
59 int final_rank = -1;
60
61 if (true == condenseddim)
62 final_rank = 1;
63 else if (4 == specialformat)// For the SOM projection, the final output of latitude/longitude rank should be 3.
64 final_rank = 3;
65 else
66 final_rank = rank;
67
68 vector<int> offset;
69 offset.resize(final_rank);
70 vector<int> count;
71 count.resize(final_rank);
72 vector<int> step;
73 step.resize(final_rank);
74
75 int nelms = -1;
76
77 // Obtain the number of the subsetted elements
78 nelms = format_constraint (offset.data(), step.data(), count.data());
79
80 // Define function pointers to handle both grid and swath Note: in
81 // this code, we only handle grid, implementing this way is to
82 // keep the same style as the read functions in other files.
83 int32 (*openfunc) (char *, intn);
84 int32 (*attachfunc) (int32, char *);
85 intn (*detachfunc) (int32);
86 intn (*fieldinfofunc) (int32, char *, int32 *, int32 *, int32 *, char *);
87 intn (*readfieldfunc) (int32, char *, int32 *, int32 *, int32 *, void *);
88
89 string datasetname;
90 openfunc = GDopen;
91 attachfunc = GDattach;
92 detachfunc = GDdetach;
93 fieldinfofunc = GDfieldinfo;
94 readfieldfunc = GDreadfield;
95 datasetname = gridname;
96
97 int32 gfid = -1;
98 int32 gridid = -1;
99
100 /* Declare projection code, zone, etc grid parameters. */
101 int32 projcode = -1;
102 int32 zone = -1;
103 int32 sphere = -1;
104 float64 params[16];
105
106 int32 xdim = 0;
107 int32 ydim = 0;
108
109 float64 upleft[2];
110 float64 lowright[2];
111
112 string cache_fpath="";
113 bool use_cache = false;
114
115 // Check if passing file IDs to data
116 if (true == check_pass_fileid_key)
117 gfid = gridfd;
118 else {
119 gfid = openfunc (const_cast < char *>(filename.c_str ()), DFACC_READ);
120 if (gfid < 0) {
121 ostringstream eherr;
122 eherr << "File " << filename.c_str () << " cannot be open.";
123 throw InternalErr (__FILE__, __LINE__, eherr.str ());
124 }
125 }
126
127 // Attach the grid id; make the grid valid.
128 gridid = attachfunc (gfid, const_cast < char *>(datasetname.c_str ()));
129 if (gridid < 0) {
130 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
131 ostringstream eherr;
132 eherr << "Grid " << datasetname.c_str () << " cannot be attached.";
133 throw InternalErr (__FILE__, __LINE__, eherr.str ());
134 }
135
136 if (false == llflag) {
137
138 // Cache
139 // Check if a BES key H4.EnableEOSGeoCacheFile is true, if yes, we will check
140 // if a lat/lon cache file exists and if we can read lat/lon from this file.
141
142 if (true == HDF4RequestHandler::get_enable_eosgeo_cachefile()) {
143
144 use_cache = true;
145 BESH4Cache *llcache = BESH4Cache::get_instance();
146
147 // Here we have a sanity check for the cached parameters:Cached directory,file prefix and cached directory size.
148 // Supposedly Hyrax BES cache feature should check this and the code exists. However, the
149 // current hyrax 1.9.7 doesn't provide this feature. KY 2014-10-24
150 string bescachedir = HDF4RequestHandler::get_cache_latlon_path();
151 string bescacheprefix = HDF4RequestHandler::get_cache_latlon_prefix();
152 long cachesize = HDF4RequestHandler::get_cache_latlon_size();
153
154 if (("" == bescachedir)||(""==bescacheprefix)||(cachesize <=0)){
155 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
156 throw InternalErr (__FILE__, __LINE__, "Either the cached dir is empty or the prefix is nullptr or the cache size is not set.");
157 }
158 else {
159 struct stat sb;
160 if (stat(bescachedir.c_str(),&sb) !=0) {
161 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
162 string err_mesg="The cached directory " + bescachedir;
163 err_mesg = err_mesg + " doesn't exist. ";
164 throw InternalErr(__FILE__,__LINE__,err_mesg);
165
166 }
167 else {
168 if (true == S_ISDIR(sb.st_mode)) {
169 if (access(bescachedir.c_str(),R_OK|W_OK|X_OK) == -1) {
170 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
171 string err_mesg="The cached directory " + bescachedir;
172 err_mesg = err_mesg + " can NOT be read,written or executable.";
173 throw InternalErr(__FILE__,__LINE__,err_mesg);
174 }
175
176 }
177 else {
178 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
179 string err_mesg="The cached directory " + bescachedir;
180 err_mesg = err_mesg + " is not a directory.";
181 throw InternalErr(__FILE__,__LINE__,err_mesg);
182
183 }
184 }
185 }
186
187 string cache_fname=HDF4RequestHandler::get_cache_latlon_prefix();
188
189 intn r = -1;
190 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
191 if (r!=0) {
192 detachfunc(gridid);
193 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
194 throw InternalErr (__FILE__, __LINE__, "GDprojinfo failed");
195 }
196
197 // Retrieve dimensions and X-Y coordinates of corners
198 if (GDgridinfo(gridid, &xdim, &ydim, upleft,
199 lowright) == -1) {
200 detachfunc(gridid);
201 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
202 throw InternalErr (__FILE__, __LINE__, "GDgridinfo failed");
203 }
204
205 // Retrieve pixel registration information
206 int32 pixreg = 0;
207 r = GDpixreginfo (gridid, &pixreg);
208 if (r != 0) {
209 detachfunc(gridid);
210 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
211 ostringstream eherr;
212 eherr << "cannot obtain grid pixel registration info.";
213 throw InternalErr (__FILE__, __LINE__, eherr.str ());
214 }
215
216 //Retrieve grid pixel origin
217 int32 origin = 0;
218 r = GDorigininfo (gridid, &origin);
219 if (r != 0) {
220 detachfunc(gridid);
221 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
222 ostringstream eherr;
223 eherr << "cannot obtain grid origin info.";
224 throw InternalErr (__FILE__, __LINE__, eherr.str ());
225 }
226
227
228 // Projection code,zone,sphere,pix,origin
229 cache_fname +=HDFCFUtil::get_int_str(projcode);
230 cache_fname +=HDFCFUtil::get_int_str(zone);
231 cache_fname +=HDFCFUtil::get_int_str(sphere);
232 cache_fname +=HDFCFUtil::get_int_str(pixreg);
233 cache_fname +=HDFCFUtil::get_int_str(origin);
234
235
236 // Xdimsize and ydimsize. Although it is rare, need to consider dim major.
237 // Whether latlon is ydim,xdim or xdim,ydim.
238 if (ydimmajor) {
239 cache_fname +=HDFCFUtil::get_int_str(ydim);
240 cache_fname +=HDFCFUtil::get_int_str(xdim);
241
242 }
243 else {
244 cache_fname +=HDFCFUtil::get_int_str(xdim);
245 cache_fname +=HDFCFUtil::get_int_str(ydim);
246 }
247
248 // upleft,lowright
249 // HDF-EOS upleft,lowright,params use DDDDMMMSSS.6 digits. So choose %17.6f.
250 cache_fname +=HDFCFUtil::get_double_str(upleft[0],17,6);
251 cache_fname +=HDFCFUtil::get_double_str(upleft[1],17,6);
252 cache_fname +=HDFCFUtil::get_double_str(lowright[0],17,6);
253 cache_fname +=HDFCFUtil::get_double_str(lowright[1],17,6);
254
255 // According to HDF-EOS2 document, only 13 parameters are used.
256 for(int ipar = 0; ipar<13;ipar++)
257 cache_fname+=HDFCFUtil::get_double_str(params[ipar],17,6);
258
259
260 cache_fpath = bescachedir + "/"+ cache_fname;
261
262 try {
263 do { // do while(0) is a trick to handle break; so ignore solarcloud's warning.
264 int expected_file_size = 0;
265 if (GCTP_CEA == projcode || GCTP_GEO == projcode)
266 expected_file_size = (xdim+ydim)*sizeof(double);
267 else if (GCTP_SOM == projcode)
268 expected_file_size = xdim*ydim*NBLOCK*2*sizeof(double);
269 else
270 expected_file_size = xdim*ydim*2*sizeof(double);
271
272 int fd = 0;
273 bool latlon_from_cache = llcache->get_data_from_cache(cache_fpath, expected_file_size,fd);
274 if (false == latlon_from_cache)
275 break;
276
277 // Get the offset of lat/lon in the cached file. Since lat is stored first and lon is stored second,
278 // so offset_1d for lat/lon is different.
279 // We still need to consider different projections. 1D,2D,3D reading.Need also to consider dim major and special format.
280 size_t offset_1d = 0;
281
282 // Get the count of the lat/lon from the cached file.
283 // Notice the data is read continuously. So starting from the offset point, we have to read all elements until the
284 // last points. The total count for the data point is bigger than the production of count and step.
285 int count_1d = 1;
286
287 if (GCTP_CEA == projcode|| GCTP_GEO== projcode) {
288
289 // It seems that for 1-D lat/lon, regardless of xdimmajor or ydimmajor. It is always Lat[YDim],Lon[XDim], check getCorrectSubset
290 // So we don't need to consider the dimension major case.
291 offset_1d = (fieldtype == 1) ?offset[0] :(ydim+offset[0]);
292 count_1d = 1+(count[0]-1)*step[0];
293 }
294 else if (GCTP_SOM == projcode) {
295
296 if (true == ydimmajor) {
297 offset_1d = (fieldtype == 1)?(offset[0]*xdim*ydim+offset[1]*xdim+offset[2])
298 :(offset[0]*xdim*ydim+offset[1]*xdim+offset[2]+expected_file_size/2/sizeof(double));
299 }
300 else {
301 offset_1d = (fieldtype == 1)?(offset[0]*xdim*ydim+offset[1]*ydim+offset[2])
302 :(offset[0]*xdim*ydim+offset[1]*ydim+offset[2]+expected_file_size/2/sizeof(double));
303 }
304
305 int total_count_dim0 = (count[0]-1)*step[0];
306 int total_count_dim1 = (count[1]-1)*step[1];
307 int total_count_dim2 = (count[2]-1)*step[2];
308 int total_dim1 = (true ==ydimmajor)?ydim:xdim;
309 int total_dim2 = (true ==ydimmajor)?xdim:ydim;
310
311 // Flatten the 3-D index to 1-D
312 // This calculation can be generalized from nD to 1D
313 // but since we only use it here. Just keep it this way.
314 count_1d = 1 + total_count_dim0*total_dim1*total_dim2 + total_count_dim1*total_dim2 + total_count_dim2;
315
316 }
317 else {// 2-D lat/lon case
318 if (true == ydimmajor)
319 offset_1d = (fieldtype == 1) ?(offset[0] * xdim + offset[1]):(expected_file_size/2/sizeof(double)+offset[0]*xdim+offset[1]);
320 else
321 offset_1d = (fieldtype == 1) ?(offset[0] * ydim + offset[1]):(expected_file_size/2/sizeof(double)+offset[0]*ydim+offset[1]);
322
323 // Flatten the 2-D index to 1-D
324 int total_count_dim0 = (count[0]-1)*step[0];
325 int total_count_dim1 = (count[1]-1)*step[1];
326 int total_dim1 = (true ==ydimmajor)?xdim:ydim;
327
328 count_1d = 1 + total_count_dim0*total_dim1 + total_count_dim1;
329 }
330
331 // Assign a vector to store lat/lon
332 vector<double> latlon_1d;
333 latlon_1d.resize(count_1d);
334
335 // Read lat/lon from the file.
336 //int fd;
337 //fd = open(cache_fpath.c_str(),O_RDONLY,0666);
338 // TODO: Use BESLog to report that the cached file cannot be read.
339 off_t fpos = lseek(fd,sizeof(double)*offset_1d,SEEK_SET);
340 if (-1 == fpos) {
341 llcache->unlock_and_close(cache_fpath);
342 llcache->purge_file(cache_fpath);
343 break;
344 }
345 ssize_t read_size = HDFCFUtil::read_vector_from_file(fd,latlon_1d,sizeof(double));
346 llcache->unlock_and_close(cache_fpath);
347 if ((-1 == read_size) || ((size_t)read_size != count_1d*sizeof(double))) {
348 llcache->purge_file(cache_fpath);
349 break;
350 }
351
352 int total_count = 1;
353 for (int i_rank = 0; i_rank<final_rank;i_rank++)
354 total_count = total_count*count[i_rank];
355
356 // We will see if there is a shortcut that the lat/lon is accessed with
357 // one-big-block. Actually this is the most common case. If we find
358 // such a case, we simply read the whole data into the latlon buffer and
359 // send it to BES.
360 if (total_count == count_1d) {
361 set_value((dods_float64*)latlon_1d.data(),nelms);
362 detachfunc(gridid);
363 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
364 return false;
365 }
366
367 vector<double>latlon;
368 latlon.resize(total_count);
369
370 // Retrieve latlon according to the projection
371 if (GCTP_CEA == projcode|| GCTP_GEO== projcode) {
372 for (int i = 0; i <total_count;i++)
373 latlon[i] = latlon_1d[i*step[0]];
374
375 }
376 else if (GCTP_SOM == projcode) {
377
378 for (int i =0; i<count[0];i++)
379 for(int j =0;j<count[1];j++)
380 for(int k=0;k<count[2];k++)
381 latlon[i*count[1]*count[2]+j*count[2]+k]=(true == ydimmajor)
382 ?latlon_1d[i*ydim*xdim*step[0]+j*xdim*step[1]+k*step[2]]
383 :latlon_1d[i*ydim*xdim*step[0]+j*ydim*step[1]+k*step[2]];
384 }
385 else {
386 for (int i =0; i<count[0];i++)
387 for(int j =0;j<count[1];j++)
388 latlon[i*count[1]+j]=(true == ydimmajor)
389 ?latlon_1d[i*xdim*step[0]+j*step[1]]
390 :latlon_1d[i*ydim*step[0]+j*step[1]];
391
392 }
393
394 set_value((dods_float64*)latlon.data(),nelms);
395 detachfunc(gridid);
396 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
397 return false;
398
399 } while (0);
400
401 }
402 catch(...) {
403 detachfunc(gridid);
404 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
405 throw;
406 }
407 }
408 }
409
410
411 // In this step, if use_cache is true, we always need to write the lat/lon into the cache.
412 // SOM projection should be calculated differently. If turning on the lat/lon cache feature, it also needs to be handled differently.
413 if (specialformat == 4) {// SOM projection
414 try {
415 CalculateSOMLatLon(gridid, offset.data(), count.data(), step.data(), nelms,cache_fpath,use_cache);
416 }
417 catch(...) {
418 detachfunc(gridid);
419 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
420 throw;
421 }
422 detachfunc(gridid);
423 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
424 return false;
425 }
426
427 // We define offset,count and step in int32 datatype.
428 vector<int32>offset32;
429 offset32.resize(rank);
430
431 vector<int32>count32;
432 count32.resize(rank);
433
434 vector<int32>step32;
435 step32.resize(rank);
436
437
438 // Obtain offset32 with the correct rank, the rank of lat/lon of
439 // GEO and CEA projections in the file may be 2 instead of 1.
440 try {
441 getCorrectSubset (offset.data(), count.data(), step.data(), offset32.data(), count32.data(), step32.data(), condenseddim, ydimmajor, fieldtype, rank);
442 }
443 catch(...) {
444 detachfunc(gridid);
445 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
446 throw;
447 }
448
449 // The following case handles when the lat/lon is not provided.
450 if (llflag == false) { // We have to calculate the lat/lon
451
452 vector<float64>latlon;
453 latlon.resize(nelms);
454
455 // If projection code etc. is not retrieved, retrieve them.
456 // When specialformat is 3, the file is a file of which the project code is set to -1, we need to skip it. KY 2014-09-11
457 if (projcode == -1 && specialformat !=3) {
458
459
460 intn r = 0;
461 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
462 if (r!=0) {
463 detachfunc(gridid);
464 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
465 throw InternalErr (__FILE__, __LINE__, "GDprojinfo failed");
466 }
467
468 // Retrieve dimensions and X-Y coordinates of corners
469 if (GDgridinfo(gridid, &xdim, &ydim, upleft,
470 lowright) == -1) {
471 detachfunc(gridid);
472 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
473 throw InternalErr (__FILE__, __LINE__, "GDgridinfo failed");
474 }
475 }
476
477
478 // Handle LAMAZ projection first.
479 if (GCTP_LAMAZ == projcode) {
480 try {
481 vector<double>latlon_all;
482 latlon_all.resize(xdim*ydim*2);
483
484 CalculateLAMAZLatLon(gridid, fieldtype, latlon.data(), latlon_all.data(),offset.data(), count.data(), step.data(), use_cache);
485 if (true == use_cache) {
486
487 BESH4Cache *llcache = BESH4Cache::get_instance();
488 llcache->write_cached_data(cache_fpath,xdim*ydim*2*sizeof(double),latlon_all);
489
490 }
491 }
492 catch(...) {
493 detachfunc(gridid);
494 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
495 throw;
496 }
497 set_value ((dods_float64 *) latlon.data(), nelms);
498 detachfunc(gridid);
499 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
500 return false;
501 }
502
503 // Aim to handle large MCD Grid such as 21600*43200 lat,lon
504 if (specialformat == 1) {
505
506 try {
507 vector<double>latlon_all;
508 latlon_all.resize(xdim+ydim);
509
510 CalculateLargeGeoLatLon(gridid, fieldtype,latlon.data(), latlon_all.data(),offset.data(), count.data(), step.data(), nelms,use_cache);
511 if (true == use_cache) {
512
513 BESH4Cache *llcache = BESH4Cache::get_instance();
514 llcache->write_cached_data(cache_fpath,(xdim+ydim)*sizeof(double),latlon_all);
515
516 }
517
518 }
519 catch(...) {
520 detachfunc(gridid);
521 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
522 throw;
523 }
524 set_value((dods_float64 *)latlon.data(),nelms);
525 detachfunc(gridid);
526 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
527
528 return false;
529 }
530
531 // Now handle other cases,note the values will be written after the if-block
532 else if (specialformat == 3) {// Have to provide latitude and longitude by ourselves
533 try {
534 CalculateSpeLatLon (gridid, fieldtype, latlon.data(), offset32.data(), count32.data(), step32.data());
535 }
536 catch(...) {
537 detachfunc(gridid);
538 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
539 throw;
540
541 }
542 detachfunc(gridid);
543 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
544 }
545 else {// This is mostly general case, it will calculate lat/lon with GDij2ll.
546
547 // Cache: check the flag and decide whether to calculate the lat/lon.
548 vector<double>latlon_all;
549
550 if (GCTP_GEO == projcode || GCTP_CEA == projcode)
551 latlon_all.resize(xdim+ydim);
552 else
553 latlon_all.resize(xdim*ydim*2);
554
555 CalculateLatLon (gridid, fieldtype, specialformat, latlon.data(),latlon_all.data(),
556 offset32.data(), count32.data(), step32.data(), nelms,use_cache);
557
558 if (true == use_cache) {
559 size_t num_item_expected = 0;
560 if (GCTP_GEO == projcode || GCTP_CEA == projcode)
561 num_item_expected = xdim + ydim;
562 else
563 num_item_expected = xdim*ydim*2;
564
565 BESH4Cache *llcache = BESH4Cache::get_instance();
566 llcache->write_cached_data(cache_fpath,num_item_expected*sizeof(double),latlon_all);
567
568 }
569
570 // The longitude values changed in the cache file is implemented in CalculateLatLon.
571 // Some longitude values need to be corrected.
572 if (speciallon && fieldtype == 2)
573 CorSpeLon(latlon.data(), nelms);
574 detachfunc(gridid);
575 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
576 }
577
578 set_value ((dods_float64 *) latlon.data(), nelms);
579
580 return false;
581 }
582
583
584 // Now lat and lon are stored as HDF-EOS2 fields. We need to read the lat and lon values from the fields.
585 int32 tmp_rank = -1;
586 vector<int32> tmp_dims;
587 tmp_dims.resize(rank);
588
589 char tmp_dimlist[1024];
590 int32 type = -1;
591 intn r = -1;
592
593 // Obtain field info.
594 r = fieldinfofunc (gridid, const_cast < char *>(fieldname.c_str ()),
595 &tmp_rank, tmp_dims.data(), &type, tmp_dimlist);
596
597 if (r != 0) {
598 detachfunc(gridid);
599 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
600 ostringstream eherr;
601 eherr << "Field " << fieldname.c_str () << " information cannot be obtained.";
602 throw InternalErr (__FILE__, __LINE__, eherr.str ());
603 }
604
605 // Retrieve dimensions and X-Y coordinates of corners
606 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
607 if (r != 0) {
608 detachfunc(gridid);
609 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
610 ostringstream eherr;
611 eherr << "Grid " << datasetname.c_str () << " information cannot be obtained.";
612 throw InternalErr (__FILE__, __LINE__, eherr.str ());
613 }
614
615 // Retrieve all GCTP projection information
616 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
617 if (r != 0) {
618 detachfunc(gridid);
619 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
620 ostringstream eherr;
621 eherr << "Grid " << datasetname.c_str () << " projection info. cannot be obtained.";
622 throw InternalErr (__FILE__, __LINE__, eherr.str ());
623 }
624
625 if (projcode != GCTP_GEO) { // Just retrieve the data like other fields
626 // We have to loop through all datatype and read the lat/lon out.
627 switch (type) {
628 case DFNT_INT8:
629 {
630 vector<int8> val;
631 val.resize(nelms);
632 r = readfieldfunc (gridid,
633 const_cast < char *>(fieldname.c_str ()),
634 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
635 if (r != 0) {
636 detachfunc(gridid);
637 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
638 ostringstream eherr;
639 eherr << "field " << fieldname.c_str () << "cannot be read.";
640 throw InternalErr (__FILE__, __LINE__, eherr.str ());
641 }
642
643 // DAP2 requires the map of SIGNED_BYTE to INT32 if
644 // SIGNED_BYTE_TO_INT32 is defined.
645#ifndef SIGNED_BYTE_TO_INT32
646 set_value ((dods_byte *) val.data(), nelms);
647#else
648 vector<int32>newval;
649 newval.resize(nelms);
650
651 for (int counter = 0; counter < nelms; counter++)
652 newval[counter] = (int32) (val[counter]);
653
654 set_value ((dods_int32 *) newval.data(), nelms);
655#endif
656
657 }
658 break;
659 case DFNT_UINT8:
660 case DFNT_UCHAR8:
661
662 {
663 vector<uint8> val;
664 val.resize(nelms);
665 r = readfieldfunc (gridid,
666 const_cast < char *>(fieldname.c_str ()),
667 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
668 if (r != 0) {
669 detachfunc(gridid);
670 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
671 ostringstream eherr;
672 eherr << "field " << fieldname.c_str () << "cannot be read.";
673 throw InternalErr (__FILE__, __LINE__, eherr.str ());
674 }
675 set_value ((dods_byte *) val.data(), nelms);
676
677 }
678 break;
679
680 case DFNT_INT16:
681 {
682 vector<int16> val;
683 val.resize(nelms);
684
685 r = readfieldfunc (gridid,
686 const_cast < char *>(fieldname.c_str ()),
687 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
688 if (r != 0) {
689 detachfunc(gridid);
690 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
691 ostringstream eherr;
692 eherr << "field " << fieldname.c_str () << "cannot be read.";
693 throw InternalErr (__FILE__, __LINE__, eherr.str ());
694 }
695
696 set_value ((dods_int16 *) val.data(), nelms);
697
698 }
699 break;
700
701 case DFNT_UINT16:
702 {
703 vector<uint16> val;
704 val.resize(nelms);
705
706 r = readfieldfunc (gridid,
707 const_cast < char *>(fieldname.c_str ()),
708 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
709 if (r != 0) {
710 detachfunc(gridid);
711 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
712 ostringstream eherr;
713 eherr << "field " << fieldname.c_str () << "cannot be read.";
714 throw InternalErr (__FILE__, __LINE__, eherr.str ());
715 }
716
717 set_value ((dods_uint16 *) val.data(), nelms);
718 }
719 break;
720
721 case DFNT_INT32:
722 {
723 vector<int32> val;
724 val.resize(nelms);
725
726 r = readfieldfunc (gridid,
727 const_cast < char *>(fieldname.c_str ()),
728 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
729 if (r != 0) {
730 detachfunc(gridid);
731 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
732 ostringstream eherr;
733 eherr << "field " << fieldname.c_str () << "cannot be read.";
734 throw InternalErr (__FILE__, __LINE__, eherr.str ());
735 }
736
737 set_value ((dods_int32 *) val.data(), nelms);
738 }
739 break;
740
741 case DFNT_UINT32:
742 {
743 vector<uint32> val;
744 val.resize(nelms);
745
746 r = readfieldfunc (gridid,
747 const_cast < char *>(fieldname.c_str ()),
748 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
749 if (r != 0) {
750 detachfunc(gridid);
751 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
752 ostringstream eherr;
753 eherr << "field " << fieldname.c_str () << "cannot be read.";
754 throw InternalErr (__FILE__, __LINE__, eherr.str ());
755 }
756 set_value ((dods_uint32 *) val.data(), nelms);
757 }
758 break;
759
760 case DFNT_FLOAT32:
761 {
762 vector<float32> val;
763 val.resize(nelms);
764
765 r = readfieldfunc (gridid,
766 const_cast < char *>(fieldname.c_str ()),
767 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
768 if (r != 0) {
769 detachfunc(gridid);
770 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
771 ostringstream eherr;
772 eherr << "field " << fieldname.c_str () << "cannot be read.";
773 throw InternalErr (__FILE__, __LINE__, eherr.str ());
774 }
775
776 set_value ((dods_float32 *) val.data(), nelms);
777 }
778 break;
779
780 case DFNT_FLOAT64:
781
782 {
783 vector<float64> val;
784 val.resize(nelms);
785
786 r = readfieldfunc (gridid,
787 const_cast < char *>(fieldname.c_str ()),
788 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
789 if (r != 0) {
790 detachfunc(gridid);
791 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
792 ostringstream eherr;
793 eherr << "field " << fieldname.c_str () << "cannot be read.";
794 throw InternalErr (__FILE__, __LINE__, eherr.str ());
795 }
796
797 set_value ((dods_float64 *) val.data(), nelms);
798 }
799 break;
800
801 default:
802 {
803 detachfunc(gridid);
804 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
805 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
806 }
807
808 }
809 }
810 else {// Only handle special cases for the Geographic Projection
811 // We find that lat/lon of the geographic projection in some
812 // files include fill values. So we recalculate lat/lon based
813 // on starting value,step values and number of steps.
814 // GDgetfillvalue will return 0 if having fill values.
815 // The other returned value indicates no fillvalue is found inside the lat or lon.
816 switch (type) {
817 case DFNT_INT8:
818 {
819 vector<int8> val;
820 val.resize(nelms);
821
822 int8 fillvalue = 0;
823
824 r = GDgetfillvalue (gridid,
825 const_cast < char *>(fieldname.c_str ()),
826 &fillvalue);
827 if (r == 0) {
828 int ifillvalue = fillvalue;
829
830 vector <int8> temp_total_val;
831 temp_total_val.resize(xdim*ydim);
832
833 r = readfieldfunc(gridid,
834 const_cast < char *>(fieldname.c_str ()),
835 nullptr, nullptr, nullptr, (void *)(temp_total_val.data()));
836
837 if (r != 0) {
838 detachfunc(gridid);
839 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
840 ostringstream eherr;
841 eherr << "field " << fieldname.c_str () << "cannot be read.";
842 throw InternalErr (__FILE__, __LINE__, eherr.str ());
843 }
844
845 try {
846 // Recalculate lat/lon for the geographic projection lat/lon that has fill values
847 HandleFillLatLon(temp_total_val, val.data(),ydimmajor,fieldtype,xdim,ydim,offset32.data(),count32.data(),step32.data(),ifillvalue);
848 }
849 catch(...) {
850 detachfunc(gridid);
851 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
852 throw;
853 }
854
855 }
856
857 else {
858
859 r = readfieldfunc (gridid,
860 const_cast < char *>(fieldname.c_str ()),
861 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
862 if (r != 0) {
863 detachfunc(gridid);
864 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
865 ostringstream eherr;
866 eherr << "field " << fieldname.c_str () << "cannot be read.";
867 throw InternalErr (__FILE__, __LINE__, eherr.str ());
868 }
869 }
870
871 if (speciallon && fieldtype == 2)
872 CorSpeLon (val.data(), nelms);
873
874
875#ifndef SIGNED_BYTE_TO_INT32
876 set_value ((dods_byte *) val.data(), nelms);
877#else
878 vector<int32>newval;
879 newval.resize(nelms);
880
881 for (int counter = 0; counter < nelms; counter++)
882 newval[counter] = (int32) (val[counter]);
883
884 set_value ((dods_int32 *) newval.data(), nelms);
885
886#endif
887
888 }
889 break;
890
891 case DFNT_UINT8:
892 case DFNT_UCHAR8:
893 {
894 vector<uint8> val;
895 val.resize(nelms);
896
897 uint8 fillvalue = 0;
898
899 r = GDgetfillvalue (gridid,
900 const_cast < char *>(fieldname.c_str ()),
901 &fillvalue);
902
903 if (r == 0) {
904
905 int ifillvalue = fillvalue;
906 vector <uint8> temp_total_val;
907 temp_total_val.resize(xdim*ydim);
908
909 r = readfieldfunc(gridid,
910 const_cast < char *>(fieldname.c_str ()),
911 nullptr, nullptr, nullptr, (void *)(temp_total_val.data()));
912
913 if (r != 0) {
914 detachfunc(gridid);
915 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
916 ostringstream eherr;
917 eherr << "field " << fieldname.c_str () << "cannot be read.";
918 throw InternalErr (__FILE__, __LINE__, eherr.str ());
919 }
920
921 try {
922 HandleFillLatLon(temp_total_val, val.data(),ydimmajor,fieldtype,xdim,ydim,offset32.data(),count32.data(),step32.data(),ifillvalue);
923 }
924 catch(...) {
925 detachfunc(gridid);
926 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
927 throw;
928 }
929
930 }
931
932 else {
933
934 r = readfieldfunc (gridid,
935 const_cast < char *>(fieldname.c_str ()),
936 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
937 if (r != 0) {
938 detachfunc(gridid);
939 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
940 ostringstream eherr;
941 eherr << "field " << fieldname.c_str () << "cannot be read.";
942 throw InternalErr (__FILE__, __LINE__, eherr.str ());
943 }
944 }
945
946 if (speciallon && fieldtype == 2)
947 CorSpeLon (val.data(), nelms);
948 set_value ((dods_byte *) val.data(), nelms);
949
950 }
951 break;
952
953 case DFNT_INT16:
954 {
955 vector<int16> val;
956 val.resize(nelms);
957
958 int16 fillvalue = 0;
959
960 r = GDgetfillvalue (gridid,
961 const_cast < char *>(fieldname.c_str ()),
962 &fillvalue);
963 if (r == 0) {
964
965 int ifillvalue = fillvalue;
966 vector <int16> temp_total_val;
967 temp_total_val.resize(xdim*ydim);
968
969 r = readfieldfunc(gridid,
970 const_cast < char *>(fieldname.c_str ()),
971 nullptr, nullptr, nullptr, (void *)(temp_total_val.data()));
972
973 if (r != 0) {
974 detachfunc(gridid);
975 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
976 ostringstream eherr;
977 eherr << "field " << fieldname.c_str () << "cannot be read.";
978 throw InternalErr (__FILE__, __LINE__, eherr.str ());
979 }
980
981 try {
982 HandleFillLatLon(temp_total_val, val.data(),ydimmajor,fieldtype,xdim,ydim,offset32.data(),count32.data(),step32.data(),ifillvalue);
983 }
984 catch(...) {
985 detachfunc(gridid);
986 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
987 throw;
988 }
989
990 }
991
992 else {
993
994 r = readfieldfunc (gridid,
995 const_cast < char *>(fieldname.c_str ()),
996 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
997 if (r != 0) {
998 detachfunc(gridid);
999 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1000 ostringstream eherr;
1001 eherr << "field " << fieldname.c_str () << "cannot be read.";
1002 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1003 }
1004 }
1005
1006
1007 if (speciallon && fieldtype == 2)
1008 CorSpeLon (val.data(), nelms);
1009
1010 set_value ((dods_int16 *) val.data(), nelms);
1011 }
1012 break;
1013
1014 case DFNT_UINT16:
1015 {
1016 uint16 fillvalue = 0;
1017 vector<uint16> val;
1018 val.resize(nelms);
1019
1020 r = GDgetfillvalue (gridid,
1021 const_cast < char *>(fieldname.c_str ()),
1022 &fillvalue);
1023
1024 if (r == 0) {
1025
1026 int ifillvalue = fillvalue;
1027
1028 vector <uint16> temp_total_val;
1029 temp_total_val.resize(xdim*ydim);
1030
1031 r = readfieldfunc(gridid,
1032 const_cast < char *>(fieldname.c_str ()),
1033 nullptr, nullptr, nullptr, (void *)(temp_total_val.data()));
1034
1035 if (r != 0) {
1036 detachfunc(gridid);
1037 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1038 ostringstream eherr;
1039 eherr << "field " << fieldname.c_str () << "cannot be read.";
1040 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1041 }
1042
1043 try {
1044 HandleFillLatLon(temp_total_val, val.data(),ydimmajor,fieldtype,xdim,ydim,offset32.data(),count32.data(),step32.data(),ifillvalue);
1045 }
1046 catch(...) {
1047 detachfunc(gridid);
1048 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1049 throw;
1050 }
1051 }
1052
1053 else {
1054
1055 r = readfieldfunc (gridid,
1056 const_cast < char *>(fieldname.c_str ()),
1057 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
1058 if (r != 0) {
1059 detachfunc(gridid);
1060 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1061 ostringstream eherr;
1062 eherr << "field " << fieldname.c_str () << "cannot be read.";
1063 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1064 }
1065 }
1066
1067 if (speciallon && fieldtype == 2)
1068 CorSpeLon (val.data(), nelms);
1069
1070 set_value ((dods_uint16 *) val.data(), nelms);
1071
1072 }
1073 break;
1074
1075 case DFNT_INT32:
1076 {
1077 vector<int32> val;
1078 val.resize(nelms);
1079
1080 int32 fillvalue = 0;
1081
1082 r = GDgetfillvalue (gridid,
1083 const_cast < char *>(fieldname.c_str ()),
1084 &fillvalue);
1085 if (r == 0) {
1086
1087 int ifillvalue = fillvalue;
1088
1089 vector <int32> temp_total_val;
1090 temp_total_val.resize(xdim*ydim);
1091
1092 r = readfieldfunc(gridid,
1093 const_cast < char *>(fieldname.c_str ()),
1094 nullptr, nullptr, nullptr, (void *)(temp_total_val.data()));
1095
1096 if (r != 0) {
1097 ostringstream eherr;
1098 detachfunc(gridid);
1099 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1100 eherr << "field " << fieldname.c_str () << "cannot be read.";
1101 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1102 }
1103
1104 try {
1105 HandleFillLatLon(temp_total_val, val.data(),ydimmajor,fieldtype,xdim,ydim,offset32.data(),count32.data(),step32.data(),ifillvalue);
1106 }
1107 catch(...) {
1108 detachfunc(gridid);
1109 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1110 throw;
1111 }
1112
1113 }
1114
1115 else {
1116
1117 r = readfieldfunc (gridid,
1118 const_cast < char *>(fieldname.c_str ()),
1119 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
1120 if (r != 0) {
1121 detachfunc(gridid);
1122 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1123 ostringstream eherr;
1124 eherr << "field " << fieldname.c_str () << "cannot be read.";
1125 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1126 }
1127
1128 }
1129 if (speciallon && fieldtype == 2)
1130 CorSpeLon (val.data(), nelms);
1131
1132 set_value ((dods_int32 *) val.data(), nelms);
1133
1134 }
1135 break;
1136
1137 case DFNT_UINT32:
1138 {
1139 vector<uint32> val;
1140 val.resize(nelms);
1141
1142 uint32 fillvalue = 0;
1143
1144 r = GDgetfillvalue (gridid,
1145 const_cast < char *>(fieldname.c_str ()),
1146 &fillvalue);
1147 if (r == 0) {
1148
1149 // this may cause overflow. Although we don't find the overflow in the NASA HDF products, may still fix it later. KY 2012-8-20
1150 int ifillvalue = (int)fillvalue;
1151 vector <uint32> temp_total_val;
1152 temp_total_val.resize(xdim*ydim);
1153 r = readfieldfunc(gridid,
1154 const_cast < char *>(fieldname.c_str ()),
1155 nullptr, nullptr, nullptr, (void *)(temp_total_val.data()));
1156
1157 if (r != 0) {
1158 detachfunc(gridid);
1159 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1160 ostringstream eherr;
1161 eherr << "field " << fieldname.c_str () << "cannot be read.";
1162 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1163 }
1164
1165 try {
1166 HandleFillLatLon(temp_total_val, val.data(),ydimmajor,fieldtype,xdim,ydim,offset32.data(),count32.data(),step32.data(),ifillvalue);
1167
1168 }
1169 catch(...) {
1170 detachfunc(gridid);
1171 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1172 throw;
1173 }
1174 }
1175
1176 else {
1177
1178 r = readfieldfunc (gridid,
1179 const_cast < char *>(fieldname.c_str ()),
1180 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
1181 if (r != 0) {
1182 detachfunc(gridid);
1183 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1184 ostringstream eherr;
1185 eherr << "field " << fieldname.c_str () << "cannot be read.";
1186 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1187 }
1188
1189 }
1190 if (speciallon && fieldtype == 2)
1191 CorSpeLon (val.data(), nelms);
1192
1193 set_value ((dods_uint32 *) val.data(), nelms);
1194
1195 }
1196 break;
1197
1198 case DFNT_FLOAT32:
1199 {
1200 vector<float32> val;
1201 val.resize(nelms);
1202
1203 float32 fillvalue =0;
1204 r = GDgetfillvalue (gridid,
1205 const_cast < char *>(fieldname.c_str ()),
1206 &fillvalue);
1207
1208
1209 if (r == 0) {
1210 // May cause overflow,not find this happen in NASA HDF files, may still need to handle later.
1211 // KY 2012-08-20
1212 auto ifillvalue =(int)fillvalue;
1213
1214 vector <float32> temp_total_val;
1215 temp_total_val.resize(xdim*ydim);
1216
1217 r = readfieldfunc(gridid,
1218 const_cast < char *>(fieldname.c_str ()),
1219 nullptr, nullptr, nullptr, (void *)(temp_total_val.data()));
1220
1221 if (r != 0) {
1222 detachfunc(gridid);
1223 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1224 ostringstream eherr;
1225 eherr << "field " << fieldname.c_str () << "cannot be read.";
1226 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1227 }
1228
1229 try {
1230 HandleFillLatLon(temp_total_val, val.data(),ydimmajor,fieldtype,xdim,ydim,offset32.data(),count32.data(),step32.data(),ifillvalue);
1231 }
1232 catch(...) {
1233 detachfunc(gridid);
1234 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1235 throw;
1236 }
1237
1238 }
1239 else {
1240
1241 r = readfieldfunc (gridid,
1242 const_cast < char *>(fieldname.c_str ()),
1243 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
1244 if (r != 0) {
1245 detachfunc(gridid);
1246 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1247 ostringstream eherr;
1248 eherr << "field " << fieldname.c_str () << "cannot be read.";
1249 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1250 }
1251
1252 }
1253 if (speciallon && fieldtype == 2)
1254 CorSpeLon (val.data(), nelms);
1255
1256 set_value ((dods_float32 *) val.data(), nelms);
1257
1258 }
1259 break;
1260
1261 case DFNT_FLOAT64:
1262
1263 {
1264 vector<float64> val;
1265 val.resize(nelms);
1266
1267 float64 fillvalue = 0;
1268 r = GDgetfillvalue (gridid,
1269 const_cast < char *>(fieldname.c_str ()),
1270 &fillvalue);
1271 if (r == 0) {
1272
1273 // May cause overflow,not find this happen in NASA HDF files, may still need to handle later.
1274 // KY 2012-08-20
1275 auto ifillvalue = (int)fillvalue;
1276 vector <float64> temp_total_val;
1277 temp_total_val.resize(xdim*ydim);
1278 r = readfieldfunc(gridid,
1279 const_cast < char *>(fieldname.c_str ()),
1280 nullptr, nullptr, nullptr, (void *)(temp_total_val.data()));
1281
1282 if (r != 0) {
1283 detachfunc(gridid);
1284 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1285 ostringstream eherr;
1286 eherr << "field " << fieldname.c_str () << "cannot be read.";
1287 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1288 }
1289
1290 try {
1291 HandleFillLatLon(temp_total_val, val.data(),ydimmajor,fieldtype,xdim,ydim,offset32.data(),count32.data(),step32.data(),ifillvalue);
1292 }
1293 catch(...) {
1294 detachfunc(gridid);
1295 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1296 throw;
1297 }
1298
1299 }
1300
1301 else {
1302
1303 r = readfieldfunc (gridid,
1304 const_cast < char *>(fieldname.c_str ()),
1305 offset32.data(), step32.data(), count32.data(), (void*)(val.data()));
1306 if (r != 0) {
1307 detachfunc(gridid);
1308 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1309 ostringstream eherr;
1310 eherr << "field " << fieldname.c_str () << "cannot be read.";
1311 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1312 }
1313
1314 }
1315 if (speciallon && fieldtype == 2)
1316 CorSpeLon (val.data(), nelms);
1317
1318 set_value ((dods_float64 *) val.data(), nelms);
1319
1320 }
1321 break;
1322 default:
1323 detachfunc(gridid);
1324 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1325 throw InternalErr (__FILE__, __LINE__, "unsupported data type.");
1326 }
1327
1328 }
1329
1330 r = detachfunc (gridid);
1331 if (r != 0) {
1332 ostringstream eherr;
1333 eherr << "Grid " << datasetname.c_str () << " cannot be detached.";
1334 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1335 }
1336
1337
1338 HDFCFUtil::close_fileid(-1,-1,gfid,-1,check_pass_fileid_key);
1339
1340 return false;
1341}
1342
1343// Standard way of DAP handlers to pass the coordinates of the subsetted region to the handlers
1344// Return the number of elements to read.
1345int
1346HDFEOS2ArrayGridGeoField::format_constraint (int *offset, int *step,
1347 int *count)
1348{
1349
1350 long nels = 1;
1351 int id = 0;
1352
1353 Dim_iter p = dim_begin ();
1354 while (p != dim_end ()) {
1355
1356 int start = dimension_start (p, true);
1357 int stride = dimension_stride (p, true);
1358 int stop = dimension_stop (p, true);
1359
1360 // Check for illegal constraint
1361 if (start > stop) {
1362 ostringstream oss;
1363 oss << "Array/Grid hyperslab start point "<< start <<
1364 " is greater than stop point " << stop <<".";
1365 throw Error(malformed_expr, oss.str());
1366 }
1367
1368 offset[id] = start;
1369 step[id] = stride;
1370 count[id] = ((stop - start) / stride) + 1; // count of elements
1371 nels *= count[id]; // total number of values for variable
1372
1373 BESDEBUG ("h4",
1374 "=format_constraint():"
1375 << "id=" << id << " offset=" << offset[id]
1376 << " step=" << step[id]
1377 << " count=" << count[id]
1378 << endl);
1379
1380 id++;
1381 p++;
1382 }// end while
1383
1384 return (int)nels;
1385}
1386
1387
1388// Calculate lat/lon based on HDF-EOS2 APIs.
1389void
1390HDFEOS2ArrayGridGeoField::CalculateLatLon (int32 gridid, int g_fieldtype,
1391 int g_specialformat,
1392 float64 * outlatlon,float64* latlon_all,
1393 const int32 * offset, const int32 * count,
1394 const int32 * step, int nelms,bool write_latlon_cache)
1395{
1396
1397 // Retrieve dimensions and X-Y coordinates of corners
1398 int32 xdim = 0;
1399 int32 ydim = 0;
1400 int r = -1;
1401 float64 upleft[2];
1402 float64 lowright[2];
1403
1404 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
1405 if (r != 0) {
1406 ostringstream eherr;
1407 eherr << "cannot obtain grid information.";
1408 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1409 }
1410
1411 // The coordinate values(MCD products) are set to -180.0, -90.0, etc.
1412 // We have to change them to DDDMMMSSS.SS format, so
1413 // we have to multiply them by 1000000.
1414 if (g_specialformat == 1) {
1415 upleft[0] = upleft[0] * 1000000;
1416 upleft[1] = upleft[1] * 1000000;
1417 lowright[0] = lowright[0] * 1000000;
1418 lowright[1] = lowright[1] * 1000000;
1419 }
1420
1421 // The coordinate values(CERES TRMM) are set to default,which are zeros.
1422 // Based on the grid names and size, we find it covers the whole global.
1423 // So we set the corner coordinates to (-180000000.00,90000000.00) and
1424 // (180000000.00,-90000000.00).
1425 if (g_specialformat == 2) {
1426 upleft[0] = 0.0;
1427 upleft[1] = 90000000.0;
1428 lowright[0] = 360000000.0;
1429 lowright[1] = -90000000.0;
1430 }
1431
1432 // Retrieve all GCTP projection information
1433 int32 projcode = 0;
1434 int32 zone = 0;
1435 int32 sphere = 0;
1436 float64 params[16];
1437
1438 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
1439 if (r != 0) {
1440 ostringstream eherr;
1441 eherr << "cannot obtain grid projection information";
1442 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1443 }
1444
1445 // Retrieve pixel registration information
1446 int32 pixreg = 0;
1447
1448 r = GDpixreginfo (gridid, &pixreg);
1449 if (r != 0) {
1450 ostringstream eherr;
1451 eherr << "cannot obtain grid pixel registration info.";
1452 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1453 }
1454
1455
1456 //Retrieve grid pixel origin
1457 int32 origin = 0;
1458
1459 r = GDorigininfo (gridid, &origin);
1460 if (r != 0) {
1461 ostringstream eherr;
1462 eherr << "cannot obtain grid origin info.";
1463 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1464 }
1465
1466 vector<int32>rows;
1467 vector<int32>cols;
1468 vector<float64>lon;
1469 vector<float64>lat;
1470 rows.resize(xdim*ydim);
1471 cols.resize(xdim*ydim);
1472 lon.resize(xdim*ydim);
1473 lat.resize(xdim*ydim);
1474
1475
1476 int i = 0;
1477 int j = 0;
1478 int k = 0;
1479
1480 if (ydimmajor) {
1481 /* Fill two arguments, rows and columns */
1482 // rows cols
1483 // /- xdim -/ /- xdim -/
1484 // 0 0 0 ... 0 0 1 2 ... x
1485 // 1 1 1 ... 1 0 1 2 ... x
1486 // ... ...
1487 // y y y ... y 0 1 2 ... x
1488
1489 for (k = j = 0; j < ydim; ++j) {
1490 for (i = 0; i < xdim; ++i) {
1491 rows[k] = j;
1492 cols[k] = i;
1493 ++k;
1494 }
1495 }
1496 }
1497 else {
1498 // rows cols
1499 // /- ydim -/ /- ydim -/
1500 // 0 1 2 ... y 0 0 0 ... y
1501 // 0 1 2 ... y 1 1 1 ... y
1502 // ... ...
1503 // 0 1 2 ... y 2 2 2 ... y
1504
1505 for (k = j = 0; j < xdim; ++j) {
1506 for (i = 0; i < ydim; ++i) {
1507 rows[k] = i;
1508 cols[k] = j;
1509 ++k;
1510 }
1511 }
1512 }
1513
1514
1515 r = GDij2ll (projcode, zone, params, sphere, xdim, ydim, upleft, lowright,
1516 xdim * ydim, rows.data(), cols.data(), lon.data(), lat.data(), pixreg, origin);
1517
1518 if (r != 0) {
1519 ostringstream eherr;
1520 eherr << "cannot calculate grid latitude and longitude";
1521 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1522 }
1523
1524 // ADDING CACHE file routine,save lon and lat to a cached file. lat first, lon second.
1525 if (true == write_latlon_cache) {
1526 if (GCTP_CEA == projcode || GCTP_GEO == projcode) {
1527 vector<double>temp_lat;
1528 vector<double>temp_lon;
1529 int32 temp_offset[2];
1530 int32 temp_count[2];
1531 int32 temp_step[2];
1532 temp_offset[0] = 0;
1533 temp_offset[1] = 0;
1534 temp_step[0] = 1;
1535 temp_step[1] = 1;
1536 if (ydimmajor) {
1537 // Latitude
1538 temp_count[0] = ydim;
1539 temp_count[1] = 1;
1540 temp_lat.resize(ydim);
1541 LatLon2DSubset(temp_lat.data(),ydim,xdim,lat.data(),temp_offset,temp_count,temp_step);
1542
1543 // Longitude
1544 temp_count[0] = 1;
1545 temp_count[1] = xdim;
1546 temp_lon.resize(xdim);
1547 LatLon2DSubset(temp_lon.data(),ydim,xdim,lon.data(),temp_offset,temp_count,temp_step);
1548
1549 for(i = 0; i<ydim;i++)
1550 latlon_all[i] = temp_lat[i];
1551
1552 // Longitude values for the simple projections are mapped to 0 to 360 and we need to map them between -180 and 180.
1553 // The routine need to be called before the latlon_all to make sure the longitude value is changed.
1554 // KY 2016-03-09, HFVHANDLER-301
1555 if (speciallon == true) {//Must also apply to the latitude case since the lat/lon is stored in one cached file
1556 CorSpeLon(temp_lon.data(),xdim);
1557 }
1558
1559 for(i = 0; i<xdim;i++)
1560 latlon_all[i+ydim] = temp_lon[i];
1561
1562 }
1563 else {
1564 // Latitude
1565 temp_count[1] = ydim;
1566 temp_count[0] = 1;
1567 temp_lat.resize(ydim);
1568 LatLon2DSubset(temp_lat.data(),xdim,ydim,lat.data(),temp_offset,temp_count,temp_step);
1569
1570 // Longitude
1571 temp_count[1] = 1;
1572 temp_count[0] = xdim;
1573 temp_lon.resize(xdim);
1574 LatLon2DSubset(temp_lon.data(),xdim,ydim,lon.data(),temp_offset,temp_count,temp_step);
1575
1576 for(i = 0; i<ydim;i++)
1577 latlon_all[i] = temp_lat[i];
1578
1579 // Longitude values for the simple projections are mapped to 0 to 360 and we need to map them between -180 and 180.
1580 // The routine need to be called before the latlon_all to make sure the longitude value is changed.
1581 // KY 2016-03-09, HFVHANDLER-301
1582 if (speciallon == true) //Must also apply to the latitude case since the lat/lon is stored in one cached file
1583 CorSpeLon(temp_lon.data(),xdim);
1584
1585 for(i = 0; i<xdim;i++)
1586 latlon_all[i+ydim] = temp_lon[i];
1587
1588 }
1589 }
1590 else {
1591 memcpy((char*)(&latlon_all[0]),lat.data(),xdim*ydim*sizeof(double));
1592 memcpy((char*)(&latlon_all[0])+xdim*ydim*sizeof(double),lon.data(),xdim*ydim*sizeof(double));
1593
1594 }
1595 }
1596
1597 // 2-D Lat/Lon, need to decompose the data for subsetting.
1598 if (nelms == (xdim * ydim)) { // no subsetting return all, for the performance reason.
1599 if (g_fieldtype == 1)
1600 memcpy (outlatlon, lat.data(), xdim * ydim * sizeof (double));
1601 else
1602 memcpy (outlatlon, lon.data(), xdim * ydim * sizeof (double));
1603 }
1604 else { // Messy subsetting case, needs to know the major dimension
1605 if (ydimmajor) {
1606 if (g_fieldtype == 1) // Lat
1607 LatLon2DSubset (outlatlon, ydim, xdim, lat.data(), offset, count,
1608 step);
1609 else // Lon
1610 LatLon2DSubset (outlatlon, ydim, xdim, lon.data(), offset, count,
1611 step);
1612 }
1613 else {
1614 if (g_fieldtype == 1) // Lat
1615 LatLon2DSubset (outlatlon, xdim, ydim, lat.data(), offset, count,
1616 step);
1617 else // Lon
1618 LatLon2DSubset (outlatlon, xdim, ydim, lon.data(), offset, count,
1619 step);
1620 }
1621 }
1622}
1623
1624
1625// Map the subset of the lat/lon buffer to the corresponding 2D array.
1626template<class T> void
1627HDFEOS2ArrayGridGeoField::LatLon2DSubset (T * outlatlon, int /*majordim */,
1628 int minordim, T * latlon,
1629 const int32 * offset, const int32 * count,
1630 const int32 * step) const
1631{
1632 int i = 0;
1633 int j = 0;
1634
1635 // do subsetting
1636 // Find the correct index
1637 int dim0count = count[0];
1638 int dim1count = count[1];
1639 int dim0index[dim0count];
1640 int dim1index[dim1count];
1641
1642 for (i = 0; i < count[0]; i++) // count[0] is the least changing dimension
1643 dim0index[i] = offset[0] + i * step[0];
1644
1645
1646 for (j = 0; j < count[1]; j++)
1647 dim1index[j] = offset[1] + j * step[1];
1648
1649 // Now assign the subsetting data
1650 int k = 0;
1651
1652 for (i = 0; i < count[0]; i++) {
1653 for (j = 0; j < count[1]; j++) {
1654
1655 outlatlon[k] = *(latlon + (dim0index[i] * minordim) + dim1index[j]);
1656 k++;
1657
1658 }
1659 }
1660}
1661
1662// Some HDF-EOS2 geographic projection lat/lon fields have fill values.
1663// This routine is used to replace those fill values by using the formula to calculate
1664// the lat/lon of the geographic projection.
1665template < class T > bool HDFEOS2ArrayGridGeoField::CorLatLon (T * latlon,
1666 int g_fieldtype,
1667 int elms,
1668 int fv)
1669{
1670
1671 // Since we only find the contiguous fill value of lat/lon from some position to the end
1672 // So to speed up the performance, the following algorithm is limited to that case.
1673
1674 // The first two values cannot be fill value.
1675 // We find a special case :the first latitude(index 0) is a special value.
1676 // So we need to have three non-fill values to calculate the increment.
1677
1678 if (elms < 3) {
1679 for (int i = 0; i < elms; i++)
1680 if ((int) (latlon[i]) == fv)
1681 return false;
1682 return true;
1683 }
1684
1685 // Number of elements is greater than 3.
1686
1687 for (int i = 0; i < 3; i++) // The first three elements should not include fill value.
1688 if ((int) (latlon[i]) == fv)
1689 return false;
1690
1691 if ((int) (latlon[elms - 1]) != fv)
1692 return true;
1693
1694 T increment = latlon[2] - latlon[1];
1695
1696 int index = 0;
1697
1698 // Find the first fill value
1699 index = findfirstfv (latlon, 0, elms - 1, fv);
1700 if (index < 2) {
1701 ostringstream eherr;
1702 eherr << "cannot calculate the fill value. ";
1703 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1704 }
1705
1706 for (int i = index; i < elms; i++) {
1707
1708 latlon[i] = latlon[i - 1] + increment;
1709
1710 // The latitude must be within (-90,90)
1711 if (i != (elms - 1) && (g_fieldtype == 1) &&
1712 ((float) (latlon[i]) < -90.0 || (float) (latlon[i]) > 90.0))
1713 return false;
1714
1715 // For longitude, since some files use (0,360)
1716 // some files use (-180,180), for simple check
1717 // we just choose (-180,360).
1718 // I haven't found longitude has missing values.
1719 if (i != (elms - 1) && (g_fieldtype == 2) &&
1720 ((float) (latlon[i]) < -180.0 || (float) (latlon[i]) > 360.0))
1721 return false;
1722 }
1723 if (g_fieldtype == 1) {
1724 if ((float) (latlon[elms - 1]) < -90.0)
1725 latlon[elms - 1] = (T)-90;
1726 if ((float) (latlon[elms - 1]) > 90.0)
1727 latlon[elms - 1] = (T)90;
1728 }
1729
1730 if (g_fieldtype == 2) {
1731 if ((float) (latlon[elms - 1]) < -180.0)
1732 latlon[elms - 1] = (T)-180.0;
1733 if ((float) (latlon[elms - 1]) > 360.0)
1734 latlon[elms - 1] = (T)360.0;
1735 }
1736 return true;
1737}
1738
1739// Make longitude (0-360) to (-180 - 180)
1740template < class T > void
1741HDFEOS2ArrayGridGeoField::CorSpeLon (T * lon, int xdim) const
1742{
1743 int i;
1744 float64 accuracy = 1e-3; // in case there is a lon value = 180.0 in the middle, make the error to be less than 1e-3.
1745 float64 temp = 0;
1746
1747 // Check if this lon. field falls to the (0-360) case.
1748 int speindex = -1;
1749
1750 for (i = 0; i < xdim; i++) {
1751 if ((double) lon[i] < 180.0)
1752 temp = 180.0 - (double) lon[i];
1753 if ((double) lon[i] > 180.0)
1754 temp = (double) lon[i] - 180.0;
1755
1756 if (temp < accuracy) {
1757 speindex = i;
1758 break;
1759 }
1760 else if ((static_cast < double >(lon[i]) < 180.0)
1761 &&(static_cast<double>(lon[i + 1]) > 180.0)) {
1762 speindex = i;
1763 break;
1764 }
1765 else
1766 continue;
1767 }
1768
1769 if (speindex != -1) {
1770 for (i = speindex + 1; i < xdim; i++) {
1771 lon[i] =
1772 static_cast < T > (static_cast < double >(lon[i]) - 360.0);
1773 }
1774 }
1775 return;
1776}
1777
1778// Get correct subsetting indexes. This is especially useful when 2D lat/lon can be condensed to 1D.
1779void
1780HDFEOS2ArrayGridGeoField::getCorrectSubset (const int *offset, const int *count,
1781 const int *step, int32 * offset32,
1782 int32 * count32, int32 * step32,
1783 bool gf_condenseddim, bool gf_ydimmajor,
1784 int gf_fieldtype, int gf_rank) const
1785{
1786
1787 if (gf_rank == 1) {
1788 offset32[0] = (int32) offset[0];
1789 count32[0] = (int32) count[0];
1790 step32[0] = (int32) step[0];
1791 }
1792 else if (gf_condenseddim) {
1793
1794 // Since offset,count and step for some dimensions will always
1795 // be 1, so first assign offset32,count32,step32 to 1.
1796 for (int i = 0; i < gf_rank; i++) {
1797 offset32[i] = 0;
1798 count32[i] = 1;
1799 step32[i] = 1;
1800 }
1801
1802 if (gf_ydimmajor && gf_fieldtype == 1) {// YDim major, User: Lat[YDim], File: Lat[YDim][XDim]
1803 offset32[0] = (int32) offset[0];
1804 count32[0] = (int32) count[0];
1805 step32[0] = (int32) step[0];
1806 }
1807 else if (gf_ydimmajor && gf_fieldtype == 2) { // YDim major, User: Lon[XDim],File: Lon[YDim][XDim]
1808 offset32[1] = (int32) offset[0];
1809 count32[1] = (int32) count[0];
1810 step32[1] = (int32) step[0];
1811 }
1812 else if (!gf_ydimmajor && gf_fieldtype == 1) {// XDim major, User: Lat[YDim], File: Lat[XDim][YDim]
1813 offset32[1] = (int32) offset[0];
1814 count32[1] = (int32) count[0];
1815 step32[1] = (int32) step[0];
1816 }
1817 else if (!gf_ydimmajor && gf_fieldtype == 2) {// XDim major, User: Lon[XDim], File: Lon[XDim][YDim]
1818 offset32[0] = (int32) offset[0];
1819 count32[0] = (int32) count[0];
1820 step32[0] = (int32) step[0];
1821 }
1822
1823 else {// errors
1824 throw InternalErr (__FILE__, __LINE__,
1825 "Lat/lon subset is wrong for condensed lat/lon");
1826 }
1827 }
1828 else {
1829 for (int i = 0; i < gf_rank; i++) {
1830 offset32[i] = (int32) offset[i];
1831 count32[i] = (int32) count[i];
1832 step32[i] = (int32) step[i];
1833 }
1834 }
1835}
1836
1837// Correct latitude and longitude that have fill values. Although I only found this
1838// happens for AIRS CO2 grids, I still implemented this as general as I can.
1839
1840template <class T> void
1841HDFEOS2ArrayGridGeoField::HandleFillLatLon(vector<T> total_latlon, T* latlon,bool gf_ydimmajor, int gf_fieldtype, int32 xdim , int32 ydim, const int32* offset, const int32* count, const int32* step, int fv) {
1842
1843 class vector <T> temp_lat;
1844 class vector <T> temp_lon;
1845
1846 if (true == gf_ydimmajor) {
1847
1848 if (1 == gf_fieldtype) {
1849 temp_lat.resize(ydim);
1850 for (int i = 0; i <(int)ydim; i++)
1851 temp_lat[i] = total_latlon[i*xdim];
1852
1853 if (false == CorLatLon(temp_lat.data(),gf_fieldtype,ydim,fv))
1854 throw InternalErr(__FILE__,__LINE__,"Cannot handle the fill values in lat/lon correctly");
1855
1856 for (int i = 0; i <(int)(count[0]); i++)
1857 latlon[i] = temp_lat[offset[0] + i* step[0]];
1858 }
1859 else {
1860
1861 temp_lon.resize(xdim);
1862 for (int i = 0; i <(int)xdim; i++)
1863 temp_lon[i] = total_latlon[i];
1864
1865
1866 if (false == CorLatLon(temp_lon.data(),gf_fieldtype,xdim,fv))
1867 throw InternalErr(__FILE__,__LINE__,"Cannot handle the fill values in lat/lon correctly");
1868
1869 for (int i = 0; i <(int)(count[1]); i++)
1870 latlon[i] = temp_lon[offset[1] + i* step[1]];
1871
1872 }
1873 }
1874 else {
1875
1876 if (1 == gf_fieldtype) {
1877 temp_lat.resize(xdim);
1878 for (int i = 0; i <(int)xdim; i++)
1879 temp_lat[i] = total_latlon[i];
1880
1881 if (false == CorLatLon(temp_lat.data(),gf_fieldtype,ydim,fv))
1882 throw InternalErr(__FILE__,__LINE__,"Cannot handle the fill values in lat/lon correctly");
1883
1884 for (int i = 0; i <(int)(count[1]); i++)
1885 latlon[i] = temp_lat[offset[1] + i* step[1]];
1886 }
1887 else {
1888
1889 temp_lon.resize(ydim);
1890 for (int i = 0; i <(int)ydim; i++)
1891 temp_lon[i] = total_latlon[i*xdim];
1892
1893
1894 if (false == CorLatLon(temp_lon.data(),gf_fieldtype,xdim,fv))
1895 throw InternalErr(__FILE__,__LINE__,"Cannot handle the fill values in lat/lon correctly");
1896
1897 for (int i = 0; i <(int)(count[0]); i++)
1898 latlon[i] = temp_lon[offset[0] + i* step[0]];
1899 }
1900
1901 }
1902}
1903
1904// A helper recursive function to find the first filled value index.
1905template < class T > int
1906HDFEOS2ArrayGridGeoField::findfirstfv (T * array, int start, int end,
1907 int fillvalue)
1908{
1909
1910 if (start == end || start == (end - 1)) {
1911 if (static_cast < int >(array[start]) == fillvalue)
1912 return start;
1913 else
1914 return end;
1915 }
1916 else {
1917 int current = (start + end) / 2;
1918
1919 if (static_cast < int >(array[current]) == fillvalue)
1920 return findfirstfv (array, start, current, fillvalue);
1921 else
1922 return findfirstfv (array, current, end, fillvalue);
1923 }
1924}
1925
1926// Calculate Special Latitude and Longitude.
1927//One MOD13C2 file doesn't provide projection code
1928// The upperleft and lowerright coordinates are all -1
1929// We have to calculate lat/lon by ourselves.
1930// Since it doesn't provide the project code, we double check their information
1931// and find that it covers the whole globe with 0.05 degree resolution.
1932// Lat. is from 90 to -90 and Lon is from -180 to 180.
1933void
1934HDFEOS2ArrayGridGeoField::CalculateSpeLatLon (int32 gridid, int gf_fieldtype,
1935 float64 * outlatlon,
1936 const int32 * offset32,
1937 const int32 * count32, const int32 * step32) const
1938{
1939
1940 // Retrieve dimensions and X-Y coordinates of corners
1941 int32 xdim = 0;
1942 int32 ydim = 0;
1943 int r = -1;
1944 float64 upleft[2];
1945 float64 lowright[2];
1946
1947 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
1948 if (r != 0) {
1949 ostringstream eherr;
1950 eherr << "cannot obtain grid information.";
1951 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1952 }
1953 //Since this is a special calcuation out of using the GDij2ll function,
1954 // the rank is always assumed to be 2 and we condense to 1. So the
1955 // count for longitude should be count[1] instead of count[0]. See function GetCorSubset
1956
1957 // Since the project parameters in StructMetadata are all set to be default, I will use
1958 // the default HDF-EOS2 cell center as the origin of the coordinate. See the HDF-EOS2 user's guide
1959 // for details. KY 2012-09-10
1960
1961 if (0 == xdim || 0 == ydim)
1962 throw InternalErr(__FILE__,__LINE__,"xdim or ydim cannot be zero");
1963
1964 if (gf_fieldtype == 1) {
1965 double latstep = 180.0 / ydim;
1966
1967 for (int i = 0; i < (int) (count32[0]); i++)
1968 outlatlon[i] = 90.0 -latstep/2 - latstep * (offset32[0] + i * step32[0]);
1969 }
1970 else {// Longitude should use count32[1] etc.
1971 double lonstep = 360.0 / xdim;
1972
1973 for (int i = 0; i < (int) (count32[1]); i++)
1974 outlatlon[i] = -180.0 + lonstep/2 + lonstep * (offset32[1] + i * step32[1]);
1975 }
1976}
1977
1978// Calculate latitude and longitude for the MISR SOM projection HDF-EOS2 product.
1979// since the latitude and longitude of the SOM projection are 3-D, so we need to handle this projection in a special way.
1980// Based on our current understanding, the third dimension size is always 180.
1981// This is according to the MISR Lat/lon calculation document
1982// at http://eosweb.larc.nasa.gov/PRODOCS/misr/DPS/DPS_v50_RevS.pdf
1983void
1984HDFEOS2ArrayGridGeoField::CalculateSOMLatLon(int32 gridid, const int *start, const int *count, const int *step, int nelms,const string & cache_fpath,bool write_latlon_cache)
1985{
1986 int32 projcode = -1;
1987 int32 zone = -1;
1988 int32 sphere = -1;
1989 float64 params[NPROJ];
1990 intn r = -1;
1991
1992 r = GDprojinfo (gridid, &projcode, &zone, &sphere, params);
1993 if (r!=0)
1994 throw InternalErr (__FILE__, __LINE__, "GDprojinfo doesn't return the correct values");
1995
1996 int MAXNDIM = 10;
1997 int32 dim[MAXNDIM];
1998 char dimlist[STRLEN];
1999 r = GDinqdims(gridid, dimlist, dim);
2000 // r is the number of dims. or 0.
2001 // So the valid returned value can be greater than 0. Only throw error when r is less than 0.
2002 if (r<0)
2003 throw InternalErr (__FILE__, __LINE__, "GDinqdims doesn't return the correct values");
2004
2005 bool is_block_180 = false;
2006 for(int i=0; i<MAXNDIM; i++)
2007 {
2008 if (dim[i]==NBLOCK)
2009 {
2010 is_block_180 = true;
2011 break;
2012 }
2013 }
2014 if (false == is_block_180) {
2015 ostringstream eherr;
2016 eherr <<"Number of Block is not " << NBLOCK ;
2017 throw InternalErr(__FILE__,__LINE__,eherr.str());
2018 }
2019
2020 int32 xdim = 0;
2021 int32 ydim = 0;
2022 float64 ulc[2];
2023 float64 lrc[2];
2024
2025 r = GDgridinfo (gridid, &xdim, &ydim, ulc, lrc);
2026 if (r!=0)
2027 throw InternalErr(__FILE__,__LINE__,"GDgridinfo doesn't return the correct values");
2028
2029
2030 float32 offset[NOFFSET];
2031 char som_rw_code[]="r";
2032 r = GDblkSOMoffset(gridid, offset, NOFFSET, som_rw_code);
2033 if (r!=0)
2034 throw InternalErr(__FILE__,__LINE__,"GDblkSOMoffset doesn't return the correct values");
2035
2036 int status = misr_init(NBLOCK, xdim, ydim, offset, ulc, lrc);
2037 if (status!=0)
2038 throw InternalErr(__FILE__,__LINE__,"misr_init doesn't return the correct values");
2039
2040 int iflg = 0;
2041 int (*inv_trans[MAXPROJ+1])(double, double, double*, double*);
2042 inv_init((int)projcode, (int)zone, (double*)params, (int)sphere, nullptr, nullptr, (int*)&iflg, inv_trans);
2043 if (iflg)
2044 throw InternalErr(__FILE__,__LINE__,"inv_init doesn't return correct values");
2045
2046 // Change to vector in the future. KY 2012-09-20
2047 double somx = 0.;
2048 double somy = 0.;
2049 double lat_r = 0.;
2050 double lon_r = 0.;
2051 int i = 0;
2052 int j = 0;
2053 int k = 0;
2054 int b = 0;
2055 int npts=0;
2056 float l = 0;
2057 float s = 0;
2058
2059 // Seems setting blockdim = 0 always, need to understand this more. KY 2012-09-20
2060 int blockdim=0; //20; //84.2115,84.2018, 84.192, ... //0 for all
2061 if (blockdim==0) //66.2263, 66.224, ....
2062 {
2063
2064 if (true == write_latlon_cache) {
2065 vector<double>latlon_all;
2066 latlon_all.resize(xdim*ydim*NBLOCK*2);
2067 for(i =1; i <NBLOCK+1;i++)
2068 for(j=0;j<xdim;j++)
2069 for(k=0;k<ydim;k++)
2070 {
2071 b = i;
2072 l = (float)j;
2073 s = (float)k;
2074 misrinv(b, l, s, &somx, &somy); /* (b,l.l,s.s) -> (X,Y) */
2075 sominv(somx, somy, &lon_r, &lat_r); /* (X,Y) -> (lat,lon) */
2076 latlon_all[npts] = lat_r*R2D;
2077 latlon_all[xdim*ydim*NBLOCK+npts] = lon_r*R2D;
2078 npts++;
2079
2080 }
2081
2082 BESH4Cache *llcache = BESH4Cache::get_instance();
2083 llcache->write_cached_data(cache_fpath,xdim*ydim*NBLOCK*2*sizeof(double),latlon_all);
2084
2085 // Send the subset of latlon to DAP.
2086 vector<double>latlon;
2087 latlon.resize(nelms);
2088
2089 // Leave the following #if 0 #endif for possible future references. KY 2022-11-28
2090#if 0
2091 //double[180*xdim*ydim];
2092 //int s1=start[0]+1, e1=s1+count[0]*step[0];
2093 //int s2=start[1], e2=s2+count[1]*step[1];
2094 //int s3=start[2], e3=s3+count[2]*step[2];
2095 //int s1=start[0]+1;
2096 //int s2=start[1];
2097 //int s3=start[2];
2098#endif
2099
2100 npts =0;
2101 for(i=0; i<count[0]; i++) //i = 1; i<180+1; i++)
2102 for(j=0; j<count[1]; j++)//j=0; j<xdim; j++)
2103 for(k=0; k<count[2]; k++)//k=0; k<ydim; k++)
2104 {
2105 if (fieldtype == 1) {
2106 latlon[npts] = latlon_all[start[0]*ydim*xdim+start[1]*ydim+start[2]+
2107 i*ydim*xdim*step[0]+j*ydim*step[1]+k*step[2]];
2108 }
2109 else {
2110 latlon[npts] = latlon_all[xdim*ydim*NBLOCK+start[0]*ydim*xdim+start[1]*ydim+start[2]+
2111 i*ydim*xdim*step[0]+j*ydim*step[1]+k*step[2]];
2112
2113 }
2114 npts++;
2115 }
2116
2117 set_value ((dods_float64 *) latlon.data(), nelms); //(180*xdim*ydim)); //nelms);
2118 }
2119 else {
2120 vector<double>latlon;
2121 latlon.resize(nelms); //double[180*xdim*ydim];
2122 int s1=start[0]+1;
2123 int e1=s1+count[0]*step[0];
2124 int s2=start[1];
2125 int e2=s2+count[1]*step[1];
2126 int s3=start[2];
2127 int e3=s3+count[2]*step[2];
2128 for(i=s1; i<e1; i+=step[0]) //i = 1; i<180+1; i++)
2129 for(j=s2; j<e2; j+=step[1])//j=0; j<xdim; j++)
2130 for(k=s3; k<e3; k+=step[2])//k=0; k<ydim; k++)
2131 {
2132 b = i;
2133 l = j;
2134 s = k;
2135 misrinv(b, l, s, &somx, &somy); /* (b,l.l,s.s) -> (X,Y) */
2136 sominv(somx, somy, &lon_r, &lat_r); /* (X,Y) -> (lat,lon) */
2137 if (fieldtype==1)
2138 latlon[npts] = lat_r*R2D;
2139 else
2140 latlon[npts] = lon_r*R2D;
2141 npts++;
2142 }
2143 set_value ((dods_float64 *) latlon.data(), nelms); //(180*xdim*ydim)); //nelms);
2144 }
2145 }
2146}
2147
2148// The following code aims to handle large MCD Grid(GCTP_GEO projection) such as 21600*43200 lat and lon.
2149// These MODIS MCD files don't follow standard format for lat/lon (DDDMMMSSS);
2150// they simply represent lat/lon as -180.0000000 or -90.000000.
2151// HDF-EOS2 library won't give the correct value based on these value.
2152// We need to calculate the latitude and longitude values.
2153void
2154HDFEOS2ArrayGridGeoField::CalculateLargeGeoLatLon(int32 gridid, int gf_fieldtype, float64* latlon, float64* latlon_all,const int *start, const int *count, const int *step, int nelms,bool write_latlon_cache) const
2155{
2156
2157 int32 xdim = 0;
2158 int32 ydim = 0;
2159 float64 upleft[2];
2160 float64 lowright[2];
2161 int r = 0;
2162 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
2163 if (r!=0) {
2164 throw InternalErr(__FILE__,__LINE__, "GDgridinfo failed");
2165 }
2166
2167 if (0 == xdim || 0 == ydim) {
2168 throw InternalErr(__FILE__,__LINE__, "xdim or ydim should not be zero. ");
2169 }
2170
2171 if (upleft[0]>180.0 || upleft[0] <-180.0 ||
2172 upleft[1]>90.0 || upleft[1] <-90.0 ||
2173 lowright[0] >180.0 || lowright[0] <-180.0 ||
2174 lowright[1] >90.0 || lowright[1] <-90.0) {
2175
2176 throw InternalErr(__FILE__,__LINE__, "lat/lon corner points are out of range. ");
2177 }
2178
2179 if (count[0] != nelms) {
2180 throw InternalErr(__FILE__,__LINE__, "rank is not 1 ");
2181 }
2182 float lat_step = (float)(lowright[1] - upleft[1])/ydim;
2183 float lon_step = (float)(lowright[0] - upleft[0])/xdim;
2184
2185 if (true == write_latlon_cache) {
2186
2187 for(int i = 0;i<ydim;i++)
2188 latlon_all[i] = upleft[1] + i*lat_step + lat_step/2;
2189
2190 for(int i = 0;i<xdim;i++)
2191 latlon_all[i+ydim] = upleft[0] + i*lon_step + lon_step/2;
2192
2193 }
2194
2195 // Treat the origin of the coordinate as the center of the cell.
2196 // This has been the setting of MCD43 data. KY 2012-09-10
2197 if (1 == gf_fieldtype) { //Latitude
2198 auto start_lat = (float)(upleft[1] + start[0] *lat_step + lat_step/2);
2199 float step_lat = lat_step *step[0];
2200 for (int i = 0; i < count[0]; i++)
2201 latlon[i] = start_lat +i *step_lat;
2202 }
2203 else { // Longitude
2204 auto start_lon = (float)(upleft[0] + start[0] *lon_step + lon_step/2);
2205 float step_lon = lon_step *step[0];
2206 for (int i = 0; i < count[0]; i++)
2207 latlon[i] = start_lon +i *step_lon;
2208 }
2209
2210}
2211
2212
2213// Calculate latitude and longitude for LAMAZ projection lat/lon products.
2214// GDij2ll returns infinite numbers over the north pole or the south pole.
2215void
2216HDFEOS2ArrayGridGeoField::CalculateLAMAZLatLon(int32 gridid, int gf_fieldtype, float64* latlon, float64* latlon_all, const int *start, const int *count, const int *step, bool write_latlon_cache)
2217{
2218 int32 xdim = 0;
2219 int32 ydim = 0;
2220 intn r = 0;
2221 float64 upleft[2];
2222 float64 lowright[2];
2223
2224 r = GDgridinfo (gridid, &xdim, &ydim, upleft, lowright);
2225 if (r != 0)
2226 throw InternalErr(__FILE__,__LINE__,"GDgridinfo failed");
2227
2228 vector<float64> tmp1;
2229 tmp1.resize(xdim*ydim);
2230 int32 tmp2[] = {0, 0};
2231 int32 tmp3[] = {xdim, ydim};
2232 int32 tmp4[] = {1, 1};
2233
2234 CalculateLatLon (gridid, gf_fieldtype, specialformat, tmp1.data(), latlon_all, tmp2, tmp3, tmp4, xdim*ydim,write_latlon_cache);
2235
2236 if (write_latlon_cache == true) {
2237
2238 vector<float64> temp_lat_all;
2239 vector<float64> lat_all;
2240 temp_lat_all.resize(xdim*ydim);
2241 lat_all.resize(xdim*ydim);
2242
2243 vector<float64> temp_lon_all;
2244 vector<float64> lon_all;
2245 temp_lon_all.resize(xdim*ydim);
2246 lon_all.resize(xdim*ydim);
2247
2248 for(int w=0; w < xdim*ydim; w++){
2249 temp_lat_all[w] = latlon_all[w];
2250 lat_all[w] = latlon_all[w];
2251 temp_lon_all[w] = latlon_all[w+xdim*ydim];
2252 lon_all[w] = latlon_all[w+xdim*ydim];
2253 }
2254
2255 // If we find infinite number among lat or lon values, we use the nearest neighbor method to calculate lat or lon.
2256 if (ydimmajor) {
2257 for(int i=0; i<ydim; i++)//Lat
2258 for(int j=0; j<xdim; j++)
2259 if (isundef_lat(lat_all[i*xdim+j]))
2260 lat_all[i*xdim+j]=nearestNeighborLatVal(temp_lat_all.data(), i, j, ydim, xdim);
2261 for(int i=0; i<ydim; i++)
2262 for(int j=0; j<xdim; j++)
2263 if (isundef_lon(lon_all[i*xdim+j]))
2264 lon_all[i*xdim+j]=nearestNeighborLonVal(temp_lon_all.data(), i, j, ydim, xdim);
2265 }
2266 else {
2267 for(int i=0; i<xdim; i++)
2268 for(int j=0; j<ydim; j++)
2269 if (isundef_lat(lat_all[i*ydim+j]))
2270 lat_all[i*ydim+j]=nearestNeighborLatVal(temp_lat_all.data(), i, j, xdim, ydim);
2271
2272 for(int i=0; i<xdim; i++)
2273 for(int j=0; j<ydim; j++)
2274 if (isundef_lon(lon_all[i*ydim+j]))
2275 lon_all[i*ydim+j]=nearestNeighborLonVal(temp_lon_all.data(), i, j, xdim, ydim);
2276
2277 }
2278
2279 for(int i = 0; i<xdim*ydim;i++) {
2280 latlon_all[i] = lat_all[i];
2281 latlon_all[i+xdim*ydim] = lon_all[i];
2282 }
2283
2284 }
2285
2286 // Need to optimize the access of LAMAZ subset
2287 vector<float64> tmp5;
2288 tmp5.resize(xdim*ydim);
2289
2290 for(int w=0; w < xdim*ydim; w++)
2291 tmp5[w] = tmp1[w];
2292
2293 // If we find infinite number among lat or lon values, we use the nearest neighbor method to calculate lat or lon.
2294 if (ydimmajor) {
2295 if (gf_fieldtype==1) {// Lat.
2296 for(int i=0; i<ydim; i++)
2297 for(int j=0; j<xdim; j++)
2298 if (isundef_lat(tmp1[i*xdim+j]))
2299 tmp1[i*xdim+j]=nearestNeighborLatVal(tmp5.data(), i, j, ydim, xdim);
2300 } else if (gf_fieldtype==2){ // Lon.
2301 for(int i=0; i<ydim; i++)
2302 for(int j=0; j<xdim; j++)
2303 if (isundef_lon(tmp1[i*xdim+j]))
2304 tmp1[i*xdim+j]=nearestNeighborLonVal(tmp5.data(), i, j, ydim, xdim);
2305 }
2306 } else { // end if (ydimmajor)
2307 if (gf_fieldtype==1) {
2308 for(int i=0; i<xdim; i++)
2309 for(int j=0; j<ydim; j++)
2310 if (isundef_lat(tmp1[i*ydim+j]))
2311 tmp1[i*ydim+j]=nearestNeighborLatVal(tmp5.data(), i, j, xdim, ydim);
2312 } else if (gf_fieldtype==2) {
2313 for(int i=0; i<xdim; i++)
2314 for(int j=0; j<ydim; j++)
2315 if (isundef_lon(tmp1[i*ydim+j]))
2316 tmp1[i*ydim+j]=nearestNeighborLonVal(tmp5.data(), i, j, xdim, ydim);
2317 }
2318 }
2319
2320 for(int i=start[0], k=0; i<start[0]+count[0]*step[0]; i+=step[0])
2321 for(int j=start[1]; j<start[1]+count[1]*step[1]; j+= step[1])
2322 latlon[k++] = tmp1[i*ydim+j];
2323
2324}
2325#endif
2326
BESFileLockingCache::unlock_and_close
virtual void unlock_and_close(const std::string &target)
Definition BESFileLockingCache.cc:745
BESFileLockingCache::purge_file
virtual void purge_file(const std::string &file)
Purge a single file from the cache.
Definition BESFileLockingCache.cc:1125
BESH4Cache::get_instance
static BESH4Cache * get_instance()
Definition BESH4MCache.cc:102
ostringstream
STL class.
vector
STL class.
HDFCFUtil::close_fileid
static void close_fileid(int32 sdfd, int32 file_id, int32 gridfd, int32 swathfd, bool pass_fileid_key)
Definition HDFCFUtil.cc:3559