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