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