11 #include "HDFSPArrayGeoField.h"
18 #include "InternalErr.h"
20 #include "HDFCFUtil.h"
21 #include "HDF4RequestHandler.h"
25 #define SIGNED_BYTE_TO_INT32 1
28 #define NUM_PIXEL_NAME "Pixels per Scan Line"
29 #define NUM_POINTS_LINE_NAME "Number of Pixel Control Points"
30 #define NUM_SCAN_LINE_NAME "Number of Scan Lines"
31 #define NUM_LAT_NAME "Number of Lines"
32 #define NUM_LON_NAME "Number of Columns"
33 #define LAT_STEP_NAME "Latitude Step"
34 #define LON_STEP_NAME "Longitude Step"
35 #define SWP_LAT_NAME "SW Point Latitude"
36 #define SWP_LON_NAME "SW Point Longitude"
39 bool HDFSPArrayGeoField::read ()
42 BESDEBUG(
"h4",
"Coming to HDFSPArrayGeoField read "<<endl);
57 nelms = format_constraint (&offset[0], &step[0], &count[0]);
60 vector<int32>offset32;
61 offset32.resize(rank);
69 for (
int i = 0; i < rank; i++) {
70 offset32[i] = (int32) offset[i];
71 count32[i] = (int32) count[i];
72 step32[i] = (int32) step[i];
82 readtrmml2_v6 (&offset32[0], &count32[0], &step32[0], nelms);
89 readtrmml3a_v6 (&offset32[0], &count32[0], &step32[0], nelms);
95 readtrmml3b_v6 (&offset32[0], &count32[0], &step32[0], nelms);
101 readtrmml3c_v6 (&offset32[0], &count32[0], &step32[0], nelms);
110 readtrmml3_v7 (&offset32[0], &step32[0], nelms);
116 readceravgsyn (&offset32[0], &count32[0], &step32[0], nelms);
123 readceres4ig (&offset32[0], &count32[0], &step32[0], nelms);
130 readcersavgid2 (&offset[0], &count[0], &step[0], nelms);
137 readceres4ig (&offset32[0], &count32[0], &step32[0], nelms);
146 readcersavgid1 (&offset[0], &count[0], &step[0], nelms);
149 else if (rank == 2) {
150 readcersavgid2 (&offset[0], &count[0], &step[0], nelms);
158 readceravgsyn (&offset32[0], &count32[0], &step32[0], nelms);
165 readcerzavg (&offset32[0], &count32[0], &step32[0], nelms);
172 readobpgl2 (&offset32[0], &count32[0], &step32[0], nelms);
179 readobpgl3 (&offset[0], &step[0], nelms);
186 throw InternalErr (__FILE__, __LINE__,
"Unsupported HDF files");
192 throw InternalErr (__FILE__, __LINE__,
"Unsupported HDF files");
207 HDFSPArrayGeoField::readtrmml2_v6 (int32 * offset32, int32 * count32,
208 int32 * step32,
int nelms)
212 string check_pass_fileid_key_str=
"H4.EnablePassFileID";
213 bool check_pass_fileid_key =
false;
214 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
217 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
221 if(
false == check_pass_fileid_key) {
222 sdid = SDstart (
const_cast < char *
>(filename.c_str ()), DFACC_READ);
225 eherr <<
"File " << filename.c_str () <<
" cannot be open.";
226 throw InternalErr (__FILE__, __LINE__, eherr.str ());
234 vector<int32>geooffset32;
235 geooffset32.resize(rank+1);
237 vector<int32>geocount32;
238 geocount32.resize(rank+1);
240 vector<int32>geostep32;
241 geostep32.resize(rank+1);
243 for (
int i = 0; i < rank; i++) {
244 geooffset32[i] = offset32[i];
245 geocount32[i] = count32[i];
246 geostep32[i] = step32[i];
249 if (fieldtype == 1) {
250 geooffset32[rank] = 0;
251 geocount32[rank] = 1;
255 if (fieldtype == 2) {
256 geooffset32[rank] = 1;
257 geocount32[rank] = 1;
261 int32 sdsindex = SDreftoindex (sdid, (int32) fieldref);
263 if (sdsindex == -1) {
266 eherr <<
"SDS index " << sdsindex <<
" is not right.";
267 throw InternalErr (__FILE__, __LINE__, eherr.str ());
270 sdsid = SDselect (sdid, sdsindex);
274 eherr <<
"SDselect failed.";
275 throw InternalErr (__FILE__, __LINE__, eherr.str ());
286 r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (
void*)(&val[0]));
291 eherr <<
"SDreaddata failed.";
292 throw InternalErr (__FILE__, __LINE__, eherr.str ());
295 #ifndef SIGNED_BYTE_TO_INT32
296 set_value ((dods_byte *) &val[0], nelms);
299 newval.resize(nelms);
301 for (
int counter = 0; counter < nelms; counter++)
302 newval[counter] = (int32) (val[counter]);
303 set_value ((dods_int32 *) &newval[0], nelms);
313 r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (
void*)(&val[0]));
318 eherr <<
"SDreaddata failed";
319 throw InternalErr (__FILE__, __LINE__, eherr.str ());
322 set_value ((dods_byte *) &val[0], nelms);
330 r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (
void*)(&val[0]));
335 eherr <<
"SDreaddata failed";
336 throw InternalErr (__FILE__, __LINE__, eherr.str ());
339 set_value ((dods_int16 *) &val[0], nelms);
347 r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (
void*)(&val[0]));
352 eherr <<
"SDreaddata failed";
353 throw InternalErr (__FILE__, __LINE__, eherr.str ());
356 set_value ((dods_uint16 *) &val[0], nelms);
363 r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (
void*)(&val[0]));
368 eherr <<
"SDreaddata failed";
369 throw InternalErr (__FILE__, __LINE__, eherr.str ());
372 set_value ((dods_int32 *) &val[0], nelms);
379 r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (
void*)(&val[0]));
384 eherr <<
"SDreaddata failed";
385 throw InternalErr (__FILE__, __LINE__, eherr.str ());
387 set_value ((dods_uint32 *) &val[0], nelms);
394 r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (
void*)(&val[0]));
399 eherr <<
"SDreaddata failed";
400 throw InternalErr (__FILE__, __LINE__, eherr.str ());
403 set_value ((dods_float32 *) &val[0], nelms);
410 r = SDreaddata (sdsid, &geooffset32[0], &geostep32[0], &geocount32[0], (
void*)(&val[0]));
415 eherr <<
"SDreaddata failed";
416 throw InternalErr (__FILE__, __LINE__, eherr.str ());
419 set_value ((dods_float64 *) &val[0], nelms);
425 throw InternalErr (__FILE__, __LINE__,
"unsupported data type.");
428 r = SDendaccess (sdsid);
432 eherr <<
"SDendaccess failed.";
433 throw InternalErr (__FILE__, __LINE__, eherr.str ());
443 HDFSPArrayGeoField::readtrmml3a_v6 (int32 * offset32, int32 * count32,
444 int32 * step32,
int nelms)
447 const float slat = 89.5;
448 const float slon = 0.5;
452 if (fieldtype == 1) {
455 float sval = slat - offset32[0];
457 while (icount < (
int) (count32[0])) {
458 val[icount] = sval - step32[0] * icount;
463 if (fieldtype == 2) {
465 float sval = slon + offset32[0];
467 while (icount < (
int) (count32[0])) {
468 val[icount] = sval + step32[0] * icount;
472 set_value ((dods_float32 *) &val[0], nelms);
479 HDFSPArrayGeoField::readtrmml3b_v6 (int32 * offset32, int32 * count32,
480 int32 * step32,
int nelms)
483 const float slat = -49.875;
484 const float slon = -179.875;
488 if (fieldtype == 1) {
491 float sval = slat + 0.25 * (int) (offset32[0]);
493 while (icount < (
int) (count32[0])) {
494 val[icount] = sval + 0.25 * (int) (step32[0]) * icount;
499 if (fieldtype == 2) {
501 float sval = slon + 0.25 * (int) (offset32[0]);
503 while (icount < (
int) (count32[0])) {
504 val[icount] = sval + 0.25 * (int) (step32[0]) * icount;
508 set_value ((dods_float32 *) &val[0], nelms);
514 HDFSPArrayGeoField::readtrmml3c_v6 (int32 * offset32, int32 * count32,
515 int32 * step32,
int nelms)
518 const float slat = -36.75;
519 const float slon = -179.75;
523 if (fieldtype == 1) {
526 float sval = slat + 0.5 * (int) (offset32[0]);
528 while (icount < (
int) (count32[0])) {
529 val[icount] = sval + 0.5 * (int) (step32[0]) * icount;
534 if (fieldtype == 2) {
536 float sval = slon + 0.5 * (int) (offset32[0]);
538 while (icount < (
int) (count32[0])) {
539 val[icount] = sval + 0.5 * (int) (step32[0]) * icount;
543 set_value ((dods_float32 *) &val[0], nelms);
548 HDFSPArrayGeoField::readtrmml3_v7 (int32 * offset32,
549 int32 * step32,
int nelms)
553 string check_pass_fileid_key_str=
"H4.EnablePassFileID";
554 bool check_pass_fileid_key =
false;
555 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
558 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
562 if(
false == check_pass_fileid_key) {
563 sdid = SDstart (
const_cast < char *
>(filename.c_str ()), DFACC_READ);
566 eherr <<
"File " << filename.c_str () <<
" cannot be open.";
567 throw InternalErr (__FILE__, __LINE__, eherr.str ());
574 string gridinfo_name =
"GridHeader";
580 throw InternalErr (__FILE__,__LINE__,
581 "The maximum number of grids to be supported in the current implementation is 9.");
585 ostringstream fieldref_str;
586 fieldref_str << fieldref;
587 gridinfo_name = gridinfo_name + fieldref_str.str();
591 int32 attr_index = 0;
592 attr_index = SDfindattr (sdid, gridinfo_name.c_str());
593 if (attr_index == FAIL) {
595 string err_mesg =
"SDfindattr failed,should find attribute "+gridinfo_name+
" .";
596 throw InternalErr (__FILE__, __LINE__, err_mesg);
599 int32 attr_dtype = 0;
602 char attr_name[H4_MAX_NC_NAME];
604 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
605 if (status == FAIL) {
607 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
610 vector<char> attr_value;
611 attr_value.resize(n_values * DFKNTsize(attr_dtype));
613 status = SDreadattr (sdid, attr_index, &attr_value[0]);
614 if (status == FAIL) {
616 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
619 float lat_start = 0;;
620 float lon_start = 0.;
627 HDFCFUtil::parser_trmm_v7_gridheader(attr_value,latsize,lonsize,
628 lat_start,lon_start,lat_res,lon_res,
false);
634 if(0 == latsize || 0 == lonsize) {
636 throw InternalErr (__FILE__, __LINE__,
"Either latitude or longitude size is 0. ");
644 for (
int i = 0; i < nelms; ++i)
645 val[i] = lat_start+offset32[0]*lat_res+lat_res/2 + i*lat_res*step32[0];
647 else if(fieldtype == 2) {
648 for (
int i = 0; i < nelms; ++i)
649 val[i] = lon_start+offset32[0]*lon_res+lon_res/2 + i*lon_res*step32[0];
652 set_value ((dods_float32 *) &val[0], nelms);
668 HDFSPArrayGeoField::readobpgl2 (int32 * offset32, int32 * count32,
669 int32 * step32,
int nelms)
672 string check_pass_fileid_key_str=
"H4.EnablePassFileID";
673 bool check_pass_fileid_key =
false;
674 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
677 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
681 if(
false == check_pass_fileid_key) {
682 sdid = SDstart (
const_cast < char *
>(filename.c_str ()), DFACC_READ);
685 eherr <<
"File " << filename.c_str () <<
" cannot be open.";
686 throw InternalErr (__FILE__, __LINE__, eherr.str ());
696 int32 attr_index = 0;
697 int32 attr_dtype = 0;
699 int32 num_pixel_data = 0;
700 int32 num_point_data = 0;
701 int32 num_scan_data = 0;
703 attr_index = SDfindattr (sdid, NUM_PIXEL_NAME);
704 if (attr_index == FAIL) {
706 string attr_name(NUM_PIXEL_NAME);
707 string err_mesg =
"SDfindattr failed,should find attribute "+attr_name+
" .";
708 throw InternalErr (__FILE__, __LINE__, err_mesg);
711 char attr_name[H4_MAX_NC_NAME];
713 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
714 if (status == FAIL) {
716 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
721 throw InternalErr (__FILE__, __LINE__,
722 "Only one value of number of scan line ");
725 status = SDreadattr (sdid, attr_index, &num_pixel_data);
726 if (status == FAIL) {
728 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
731 attr_index = SDfindattr (sdid, NUM_POINTS_LINE_NAME);
732 if (attr_index == FAIL) {
734 string attr_name(NUM_POINTS_LINE_NAME);
735 string err_mesg =
"SDfindattr failed,should find attribute "+attr_name+
" .";
736 throw InternalErr (__FILE__, __LINE__, err_mesg);
741 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
742 if (status == FAIL) {
744 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
749 throw InternalErr (__FILE__, __LINE__,
750 "Only one value of number of point ");
753 status = SDreadattr (sdid, attr_index, &num_point_data);
754 if (status == FAIL) {
756 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
759 attr_index = SDfindattr (sdid, NUM_SCAN_LINE_NAME);
760 if (attr_index == FAIL) {
762 string attr_name(NUM_SCAN_LINE_NAME);
763 string err_mesg =
"SDfindattr failed,should find attribute "+attr_name+
" .";
764 throw InternalErr (__FILE__, __LINE__, err_mesg);
769 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
770 if (status == FAIL) {
772 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
777 throw InternalErr (__FILE__, __LINE__,
"Only one value of number of point ");
780 status = SDreadattr (sdid, attr_index, &num_scan_data);
781 if (status == FAIL) {
783 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
787 if ( 0 == num_scan_data || 0 == num_point_data || 0 == num_pixel_data) {
789 throw InternalErr (__FILE__, __LINE__,
"num_scan or num_point or num_pixel should not be zero. ");
792 if ( 1 == num_point_data && num_pixel_data != 1) {
794 throw InternalErr (__FILE__, __LINE__,
795 "num_point is 1 and num_pixel is not 1, interpolation cannot be done ");
798 bool compmapflag =
false;
799 if (num_pixel_data == num_point_data)
802 int32 sdsindex = SDreftoindex (sdid, (int32) fieldref);
803 if (sdsindex == -1) {
806 eherr <<
"SDS index " << sdsindex <<
" is not right.";
807 throw InternalErr (__FILE__, __LINE__, eherr.str ());
810 sdsid = SDselect (sdid, sdsindex);
814 eherr <<
"SDselect failed.";
815 throw InternalErr (__FILE__, __LINE__, eherr.str ());
832 throw InternalErr (__FILE__, __LINE__,
"datatype is not float, unsupported.");
838 r = SDreaddata (sdsid, offset32, step32, count32, &val[0]);
843 eherr <<
"SDreaddata failed";
844 throw InternalErr (__FILE__, __LINE__, eherr.str ());
849 int total_elm = num_scan_data * num_point_data;
850 vector<float32>orival;
851 orival.resize(total_elm);
852 int32 all_start[2],all_edges[2];
856 all_edges[0] = num_scan_data;
857 all_edges[1] = num_point_data;
859 r = SDreaddata (sdsid, all_start, NULL, all_edges,
865 eherr <<
"SDreaddata failed";
866 throw InternalErr (__FILE__, __LINE__, eherr.str ());
868 int interpolate_elm = num_scan_data *num_pixel_data;
870 vector<float32> interp_val;
871 interp_val.resize(interpolate_elm);
877 if (num_pixel_data % num_point_data == 0)
878 tempseg = num_pixel_data / num_point_data;
880 tempseg = num_pixel_data / num_point_data + 1;
883 (num_pixel_data%num_point_data)?(num_pixel_data-1-(tempseg*(num_point_data-2))):tempseg;
885 if ( 0 == last_tempseg || 0 == tempseg) {
888 throw InternalErr(__FILE__,__LINE__,
"Segments cannot be zero");
891 int interp_val_index = 0;
893 for (
int i = 0; i <num_scan_data; i++) {
896 for(
int j =0; j <num_point_data -2; j ++) {
897 tempdiff = orival[i*num_point_data+j+1] - orival[i*num_point_data+j];
898 for (
int k = 0; k <tempseg; k++) {
899 interp_val[interp_val_index] = orival[i*num_point_data+j] +
906 tempdiff = orival[i*num_point_data+num_point_data-1]
907 -orival[i*num_point_data+num_point_data-2];
908 for (
int k = 0; k <last_tempseg; k++) {
909 interp_val[interp_val_index] = orival[i*num_point_data+num_point_data-2] +
910 tempdiff/last_tempseg *k;
914 interp_val[interp_val_index] = orival[i*num_point_data+num_point_data-1];
919 LatLon2DSubset(&val[0],num_scan_data,num_pixel_data,&interp_val[0],
920 offset32,count32,step32);
931 int32 realcount2 = oricount32[1] - 1;
933 for (i = 0; i < (int) count32[0]; i++) {
934 for (j = 0; j < (int) realcount2 - 1; j++) {
936 orival[i * (int) oricount32[1] + j + 1] -
937 orival[i * (
int) oricount32[1] + j];
938 for (k = 0; k < tempnewseg; k++) {
939 val[i * (int) count32[1] + j * tempnewseg + k] =
940 orival[i * (
int) oricount32[1] + j] +
941 tempdiff / tempnewseg * k;
945 orival[i * (int) oricount32[1] + j + 1] -
946 orival[i * (
int) oricount32[1] + j];
952 tempnewseg * (
int) (realcount2 - 1));
953 for (k2 = 0; k2 <= lastseg; k2++)
954 val[i * (
int) count32[1] + j * tempnewseg + k2] =
955 orival[i * (int) oricount32[1] + j] +
956 tempdiff / lastseg * k2;
959 delete[](float32 *) orival;
963 set_value ((dods_float32 *) &val[0], nelms);
969 throw InternalErr (__FILE__, __LINE__,
"unsupported data type.");
972 r = SDendaccess (sdsid);
976 eherr <<
"SDendaccess failed.";
977 throw InternalErr (__FILE__, __LINE__, eherr.str ());
987 HDFSPArrayGeoField::readobpgl3 (
int *offset,
int *step,
int nelms)
991 string check_pass_fileid_key_str=
"H4.EnablePassFileID";
992 bool check_pass_fileid_key =
false;
993 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
996 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
1000 if(
false == check_pass_fileid_key) {
1001 sdid = SDstart (
const_cast < char *
>(filename.c_str ()), DFACC_READ);
1003 ostringstream eherr;
1004 eherr <<
"File " << filename.c_str () <<
" cannot be open.";
1005 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1014 int32 attr_index = 0;
1015 int32 attr_dtype = 0;
1017 int32 num_lat_data = 0;
1018 int32 num_lon_data = 0;
1019 float32 lat_step = 0.;
1020 float32 lon_step = 0.;
1021 float32 swp_lat = 0.;
1022 float32 swp_lon = 0.;
1025 attr_index = SDfindattr (sdid, NUM_LAT_NAME);
1026 if (attr_index == FAIL) {
1028 string attr_name(NUM_LAT_NAME);
1029 string err_mesg =
"SDfindattr failed,should find attribute "+attr_name+
" .";
1030 throw InternalErr (__FILE__, __LINE__, err_mesg);
1034 char attr_name[H4_MAX_NC_NAME];
1036 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1037 if (status == FAIL) {
1039 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
1042 if (n_values != 1) {
1044 throw InternalErr (__FILE__, __LINE__,
"Only should have one value ");
1047 status = SDreadattr (sdid, attr_index, &num_lat_data);
1048 if (status == FAIL) {
1050 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
1054 attr_index = SDfindattr (sdid, NUM_LON_NAME);
1055 if (attr_index == FAIL) {
1057 string attr_name2(NUM_LON_NAME);
1058 string err_mesg =
"SDfindattr failed,should find attribute "+attr_name2+
" .";
1059 throw InternalErr (__FILE__, __LINE__, err_mesg);
1064 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1065 if (status == FAIL) {
1067 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
1070 if (n_values != 1) {
1072 throw InternalErr (__FILE__, __LINE__,
"Only should have one value ");
1075 status = SDreadattr (sdid, attr_index, &num_lon_data);
1076 if (status == FAIL) {
1078 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
1082 attr_index = SDfindattr (sdid, LAT_STEP_NAME);
1083 if (attr_index == FAIL) {
1085 string attr_name2(LAT_STEP_NAME);
1086 string err_mesg =
"SDfindattr failed,should find attribute "+attr_name2+
" .";
1087 throw InternalErr (__FILE__, __LINE__, err_mesg);
1092 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1093 if (status == FAIL) {
1095 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
1098 if (n_values != 1) {
1100 throw InternalErr (__FILE__, __LINE__,
"Only should have one value ");
1103 status = SDreadattr (sdid, attr_index, &lat_step);
1104 if (status == FAIL) {
1106 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
1110 attr_index = SDfindattr (sdid, LON_STEP_NAME);
1111 if (attr_index == FAIL) {
1113 string attr_name2(LON_STEP_NAME);
1114 string err_mesg =
"SDfindattr failed,should find attribute "+attr_name2+
" .";
1115 throw InternalErr (__FILE__, __LINE__, err_mesg);
1120 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1121 if (status == FAIL) {
1123 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
1126 if (n_values != 1) {
1128 throw InternalErr (__FILE__, __LINE__,
"Only should have one value ");
1131 status = SDreadattr (sdid, attr_index, &lon_step);
1132 if (status == FAIL) {
1134 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
1138 attr_index = SDfindattr (sdid, SWP_LAT_NAME);
1139 if (attr_index == FAIL) {
1141 string attr_name2(SWP_LAT_NAME);
1142 string err_mesg =
"SDfindattr failed,should find attribute "+attr_name2+
" .";
1143 throw InternalErr (__FILE__, __LINE__, err_mesg);
1148 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1149 if (status == FAIL) {
1151 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
1154 if (n_values != 1) {
1156 throw InternalErr (__FILE__, __LINE__,
"Only should have one value ");
1159 status = SDreadattr (sdid, attr_index, &swp_lat);
1160 if (status == FAIL) {
1162 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
1166 attr_index = SDfindattr (sdid, SWP_LON_NAME);
1167 if (attr_index == FAIL) {
1169 string attr_name2(SWP_LON_NAME);
1170 string err_mesg =
"SDfindattr failed,should find attribute "+attr_name2+
" .";
1171 throw InternalErr (__FILE__, __LINE__, err_mesg);
1176 SDattrinfo (sdid, attr_index, attr_name, &attr_dtype, &n_values);
1177 if (status == FAIL) {
1179 throw InternalErr (__FILE__, __LINE__,
"SDattrinfo failed ");
1182 if (n_values != 1) {
1184 throw InternalErr (__FILE__, __LINE__,
"Only should have one value ");
1187 status = SDreadattr (sdid, attr_index, &swp_lon);
1188 if (status == FAIL) {
1190 throw InternalErr (__FILE__, __LINE__,
"SDreadattr failed ");
1194 if (fieldtype == 1) {
1196 vector<float32> allval;
1197 allval.resize(num_lat_data);
1201 for (
int j = 0; j < num_lat_data; j++)
1202 allval[j] = (num_lat_data - j - 1) * lat_step + swp_lat;
1207 for (
int k = 0; k < nelms; k++)
1208 val[k] = allval[(
int) (offset[0] + step[0] * k)];
1210 set_value ((dods_float32 *) &val[0], nelms);
1213 if (fieldtype == 2) {
1215 vector<float32>allval;
1216 allval.resize(num_lon_data);
1217 for (
int j = 0; j < num_lon_data; j++)
1218 allval[j] = swp_lon + j * lon_step;
1222 for (
int k = 0; k < nelms; k++)
1223 val[k] = allval[(
int) (offset[0] + step[0] * k)];
1225 set_value ((dods_float32 *) &val[0], nelms);
1239 HDFSPArrayGeoField::readcersavgid2 (
int *offset,
int *count,
int *step,
1243 const int dimsize0 = 180;
1244 const int dimsize1 = 360;
1245 float32 val[count[0]][count[1]];
1246 float32 orival[dimsize0][dimsize1];
1251 if (fieldtype == 1) {
1252 for (
int i = 0; i < dimsize0; i++)
1253 for (
int j = 0; j < dimsize1; j++)
1254 orival[i][j] = 89.5 - i;
1256 for (
int i = 0; i < count[0]; i++) {
1257 for (
int j = 0; j < count[1]; j++) {
1259 orival[offset[0] + step[0] * i][offset[1] + step[1] * j];
1264 if (fieldtype == 2) {
1284 for (j = 0; j < dimsize1; j++) {
1285 orival[0][j] = -179.5;
1286 orival[dimsize0 - 1][j] = -179.5;
1295 for (i = 0; i < latrange; i++)
1296 for (j = 0; j < (dimsize1 / lonextent); j++)
1297 for (k = 0; k < lonextent; k++)
1298 orival[i + latindex_north][j * lonextent + k] =
1299 -179.5 + lonextent * j;
1302 latindex_south = dimsize0 - 1 - latrange;
1303 for (i = 0; i < latrange; i++)
1304 for (j = 0; j < (dimsize1 / lonextent); j++)
1305 for (k = 0; k < lonextent; k++)
1306 orival[i + latindex_south][j * lonextent + k] =
1307 -179.5 + lonextent * j;
1310 latindex_north = latindex_north + latrange;
1314 for (i = 0; i < latrange; i++)
1315 for (j = 0; j < (dimsize1 / lonextent); j++)
1316 for (k = 0; k < lonextent; k++)
1317 orival[i + latindex_north][j * lonextent + k] =
1318 -179.5 + lonextent * j;
1321 latindex_south = latindex_south - latrange;
1322 for (i = 0; i < latrange; i++)
1323 for (j = 0; j < (dimsize1 / lonextent); j++)
1324 for (k = 0; k < lonextent; k++)
1325 orival[i + latindex_south][j * lonextent + k] =
1326 -179.5 + lonextent * j;
1329 latindex_north = latindex_north + latrange;
1333 for (i = 0; i < latrange; i++)
1334 for (j = 0; j < (dimsize1 / lonextent); j++)
1335 for (k = 0; k < lonextent; k++)
1336 orival[i + latindex_north][j * lonextent + k] =
1337 -179.5 + lonextent * j;
1340 latindex_south = latindex_south - latrange;
1342 for (i = 0; i < latrange; i++)
1343 for (j = 0; j < (dimsize1 / lonextent); j++)
1344 for (k = 0; k < lonextent; k++)
1345 orival[i + latindex_south][j * lonextent + k] =
1346 -179.5 + lonextent * j;
1349 latindex_north = latindex_north + latrange;
1353 for (i = 0; i < latrange; i++)
1354 for (j = 0; j < (dimsize1 / lonextent); j++)
1355 for (k = 0; k < lonextent; k++)
1356 orival[i + latindex_north][j * lonextent + k] =
1357 -179.5 + lonextent * j;
1360 latindex_south = latindex_south - latrange;
1361 for (i = 0; i < latrange; i++)
1362 for (j = 0; j < (dimsize1 / lonextent); j++)
1363 for (k = 0; k < lonextent; k++)
1364 orival[i + latindex_south][j * lonextent + k] =
1365 -179.5 + lonextent * j;
1367 for (i = 0; i < count[0]; i++) {
1368 for (j = 0; j < count[1]; j++) {
1370 orival[offset[0] + step[0] * i][offset[1] + step[1] * j];
1374 set_value ((dods_float32 *) (&val[0][0]), nelms);
1380 HDFSPArrayGeoField::readcersavgid1 (
int *offset,
int *count,
int *step,
1387 if (fieldtype == 1) {
1388 const int dimsize0 = 180;
1389 float32 val[count[0]];
1390 float32 orival[dimsize0];
1392 for (
int i = 0; i < dimsize0; i++)
1393 orival[i] = 89.5 - i;
1395 for (
int i = 0; i < count[0]; i++) {
1396 val[i] = orival[offset[0] + step[0] * i];
1398 set_value ((dods_float32 *) (&val[0]), nelms);
1402 if (fieldtype == 2) {
1407 throw InternalErr (__FILE__, __LINE__,
1408 "the number of element must be 1");
1409 set_value ((dods_float32 *) (&val), nelms);
1418 HDFSPArrayGeoField::readceravgsyn (int32 * offset32, int32 * count32,
1419 int32 * step32,
int nelms)
1423 string check_pass_fileid_key_str=
"H4.EnablePassFileID";
1424 bool check_pass_fileid_key =
false;
1425 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
1427 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
1431 if(
false == check_pass_fileid_key) {
1432 sdid = SDstart (
const_cast < char *
>(filename.c_str ()), DFACC_READ);
1434 ostringstream eherr;
1435 eherr <<
"File " << filename.c_str () <<
" cannot be open.";
1436 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1445 int32 sdsindex = SDreftoindex (sdid, fieldref);
1447 if (sdsindex == -1) {
1449 ostringstream eherr;
1450 eherr <<
"SDS index " << sdsindex <<
" is not right.";
1451 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1454 sdsid = SDselect (sdid, sdsindex);
1457 ostringstream eherr;
1458 eherr <<
"SDselect failed.";
1459 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1473 SDendaccess (sdsid);
1475 throw InternalErr (__FILE__, __LINE__,
1476 "datatype is not float, unsupported.");
1481 r = SDreaddata (sdsid, offset32, step32, count32, &val[0]);
1483 SDendaccess (sdsid);
1485 ostringstream eherr;
1486 eherr <<
"SDreaddata failed";
1487 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1490 if (fieldtype == 1) {
1491 for (i = 0; i < nelms; i++)
1492 val[i] = 90 - val[i];
1494 if (fieldtype == 2) {
1495 for (i = 0; i < nelms; i++)
1497 val[i] = val[i] - 360.0;
1500 set_value ((dods_float32 *) &val[0], nelms);
1508 r = SDreaddata (sdsid, offset32, step32, count32, &val[0]);
1510 SDendaccess (sdsid);
1512 ostringstream eherr;
1513 eherr <<
"SDreaddata failed";
1514 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1517 if (fieldtype == 1) {
1518 for (i = 0; i < nelms; i++)
1519 val[i] = 90 - val[i];
1521 if (fieldtype == 2) {
1522 for (i = 0; i < nelms; i++)
1524 val[i] = val[i] - 360.0;
1526 set_value ((dods_float64 *) &val[0], nelms);
1530 SDendaccess (sdsid);
1532 throw InternalErr (__FILE__, __LINE__,
"unsupported data type.");
1535 r = SDendaccess (sdsid);
1537 ostringstream eherr;
1538 eherr <<
"SDendaccess failed.";
1539 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1548 HDFSPArrayGeoField::readceres4ig (int32 * offset32, int32 * count32,
1549 int32 * step32,
int nelms)
1552 string check_pass_fileid_key_str=
"H4.EnablePassFileID";
1553 bool check_pass_fileid_key =
false;
1554 check_pass_fileid_key = HDFCFUtil::check_beskeys(check_pass_fileid_key_str);
1557 bool check_pass_fileid_key = HDF4RequestHandler::get_pass_fileid();
1561 if(
false == check_pass_fileid_key) {
1562 sdid = SDstart (
const_cast < char *
>(filename.c_str ()), DFACC_READ);
1564 ostringstream eherr;
1565 eherr <<
"File " << filename.c_str () <<
" cannot be open.";
1566 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1576 int32 sdsindex = SDreftoindex (sdid, (int32) fieldref);
1577 if (sdsindex == -1) {
1579 ostringstream eherr;
1580 eherr <<
"SDS index " << sdsindex <<
" is not right.";
1581 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1584 sdsid = SDselect (sdid, sdsindex);
1587 ostringstream eherr;
1588 eherr <<
"SDselect failed.";
1589 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1593 int32 sds_dtype = 0;
1595 char sdsname[H4_MAX_NC_NAME];
1596 int32 dim_sizes[H4_MAX_VAR_DIMS];
1599 SDgetinfo (sdsid, sdsname, &sdsrank, dim_sizes, &sds_dtype, &n_attrs);
1601 SDendaccess (sdsid);
1603 ostringstream eherr;
1604 eherr <<
"SDgetinfo failed.";
1605 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1608 vector<int32> orioffset32;
1609 vector<int32> oricount32;
1610 vector<int32> oristep32;
1611 orioffset32.resize(sdsrank);
1612 oricount32.resize(sdsrank);
1613 oristep32.resize(sdsrank);
1617 switch (sds_dtype) {
1627 SDendaccess (sdsid);
1629 throw InternalErr (__FILE__, __LINE__,
1630 "datatype is not float, unsupported.");
1633 vector<float32> val;
1635 if (fieldtype == 1) {
1636 if (sptype == CER_CGEO) {
1638 SDendaccess (sdsid);
1640 throw InternalErr (__FILE__, __LINE__,
1641 "For CER_ISCCP-D2like-GEO case, lat/lon must be 3-D");
1646 orioffset32[1] = offset32[0];
1649 oricount32[1] = count32[0];
1652 oristep32[1] = step32[0];
1655 if (sptype == CER_ES4) {
1657 SDendaccess (sdsid);
1659 throw InternalErr (__FILE__, __LINE__,
1660 "For CER_ES4 case, lat/lon must be 2-D");
1663 orioffset32[0] = offset32[0];
1665 oricount32[0] = count32[0];
1667 oristep32[0] = step32[0];
1672 if (fieldtype == 2) {
1673 if (sptype == CER_CGEO) {
1675 SDendaccess (sdsid);
1677 throw InternalErr (__FILE__, __LINE__,
1678 "For CER_ISCCP-D2like-GEO case, lat/lon must be 3-D");
1681 orioffset32[2] = offset32[0];
1684 oricount32[2] = count32[0];
1687 oristep32[2] = step32[0];
1690 if (sptype == CER_ES4) {
1692 SDendaccess (sdsid);
1694 throw InternalErr (__FILE__, __LINE__,
1695 "For CER_ES4 case, lat/lon must be 2-D");
1697 orioffset32[1] = offset32[0];
1699 oricount32[1] = count32[0];
1701 oristep32[1] = step32[0];
1706 r = SDreaddata (sdsid, &orioffset32[0], &oristep32[0], &oricount32[0], &val[0]);
1708 SDendaccess (sdsid);
1710 ostringstream eherr;
1711 eherr <<
"SDreaddata failed";
1712 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1716 for (
int i = 0; i < nelms; i++)
1717 val[i] = 90 - val[i];
1718 if (fieldtype == 2) {
1723 for (
int i = 0; i < nelms; i++)
1725 val[i] = val[i] - 360.0;
1728 set_value ((dods_float32 *) &val[0], nelms);
1732 SDendaccess (sdsid);
1734 throw InternalErr (__FILE__, __LINE__,
"unsupported data type.");
1737 r = SDendaccess (sdsid);
1739 ostringstream eherr;
1740 eherr <<
"SDendaccess failed.";
1741 throw InternalErr (__FILE__, __LINE__, eherr.str ());
1749 HDFSPArrayGeoField::readcerzavg (int32 * offset32, int32 * count32,
1750 int32 * step32,
int nelms)
1752 if (fieldtype == 1) {
1757 float32 latstep = 1.0;
1759 for (
int i = 0; i < nelms; i++)
1761 89.5 - ((
int) (offset32[0]) +
1762 ((
int) (step32[0])) * i) * latstep;
1763 set_value ((dods_float32 *) &val[0], nelms);
1766 if (fieldtype == 2) {
1767 if (count32[0] != 1 || nelms != 1)
1768 throw InternalErr (__FILE__, __LINE__,
1769 "Longitude should only have one value for zonal mean");
1773 set_value ((dods_float32 *) & val, nelms);
1781 HDFSPArrayGeoField::format_constraint (
int *offset,
int *step,
int *count)
1786 Dim_iter p = dim_begin ();
1787 while (p != dim_end ()) {
1789 int start = dimension_start (p,
true);
1790 int stride = dimension_stride (p,
true);
1791 int stop = dimension_stop (p,
true);
1796 oss <<
"Array/Grid hyperslab start point "<< start <<
1797 " is greater than stop point " << stop <<
".";
1798 throw Error(malformed_expr, oss.str());
1803 count[id] = ((stop - start) / stride) + 1;
1807 "=format_constraint():"
1808 <<
"id=" <<
id <<
" offset=" << offset[
id]
1809 <<
" step=" << step[
id]
1810 <<
" count=" << count[
id]
1823 template <
typename T >
1824 void HDFSPArrayGeoField::LatLon2DSubset (T * outlatlon,
1835 T templatlonptr[majordim][minordim] (typeof templatlonptr) latlon;
1844 const int dim0count = count[0];
1845 const int dim1count = count[1];
1846 int dim0index[dim0count];
1847 int dim1index[dim1count];
1849 for (i = 0; i < count[0]; i++)
1850 dim0index[i] = offset[0] + i * step[0];
1852 for (j = 0; j < count[1]; j++)
1853 dim1index[j] = offset[1] + j * step[1];
1858 for (i = 0; i < count[0]; i++) {
1859 for (j = 0; j < count[1]; j++) {
1861 outlatlon[k] = templatlonptr[dim0index[i]][dim1index[j]];
1863 outlatlon[k] = *(latlon + (dim0index[i] * minordim) + dim1index[j]);