#include #include #include #include #include "grib_api.h" #include #include #include #include #include #include #include #include #include "test.h" #if 0 double maxRelativeError = 1e-4; /*1e-7;*/ double maxAbsoluteError = 1e-50; #else const double maxRelativeError = 1e-7; /*1e-7;*/ const double maxAbsoluteError = 0; #endif #define ISECTION_0 2 #define ISECTION_1 1024 /* beware of for ocean data */ #define ISECTION_2 3000 #define ISECTION_3 2 #define ISECTION_4 512 #define RSECTION_2 512 #define RSECTION_3 2 #define RSECTION_4 1 struct { int code; char *message; } gribex_errors[] = { {201 ," An invalid coding/decoding option was requested.",}, {202 ,"The number of bits per data value exceeds a word length.",}, {203 ,"Missing data indicated and data field contains non-zero values.",}, {204 ,"Number of bits per data value is invalid.",}, {205 ,"Number of bits per data value is same as the number of bits per computer word.",}, {301 ,"Error in encoding or decoding the letters GRIB.",}, {302 ,"Error decoding the length of the GRIB message.",}, {303 ,"Error encoding or decoding the GRIB edition number.",}, {304 ,"Error decoding octets 22 and 23, experimental edition check.",}, {305 ,"Input data is not GRIB or pseudo-GRIB.",}, {401 ,"Error encoding or decoding the length of section 1.",}, {402 ,"Error encoding or decoding the parameter version number.",}, {403 ,"Error encoding or decoding six header fields in section 1, from 'identification of centre' to 'indicator of type of level'.",}, {404 ,"Error encoding or decoding height, pressure, etc of levels.",}, {405 ,"Error encoding or decoding six fields in section 1, from 'year of century' to 'indicator of unit of time range'.",}, {406 ,"Error encoding or decoding period of time.",}, {407 ,"Error encoding or decoding time range indicator.",}, {408 ,"Error encoding or decoding number averaged.",}, {409 ,"Error encoding or decoding number missing from averages etc.",}, {410 ,"Error encoding or decoding century of data or reserved field.",}, {411 ,"Error encoding or decoding units decimal scale factor.",}, {412 ,"Error encoding or decoding ECMWF local data.",}, {413 ,"Grib edition not handled.",}, {499 ,"Error found when checking values for section 1 against expected valid GRIB values.",}, {501 ,"Error encoding or decoding length of section 2.",}, {502 ,"Error encoding or decoding number of vertical coordinate parameters.",}, {503 ,"Error encoding or decoding location of list of vertical coordinate parameters or list of numbers of points.",}, {504 ,"Error encoding or decoding data representation type.",}, {505 ,"Error encoding or decoding number of points along a parallel or meridian.",}, {506 ,"Error encoding or decoding latitude or longitude of the first grid point.",}, {507 ,"Error encoding or decoding the components flag.",}, {508 ,"Error encoding or decoding latitude or longitude of the last grid point.",}, {509 ,"Error encoding or decoding the i direction increment.",}, {510 ,"Error encoding or decoding the number of parallels between the pole and the Equator.",}, {511 ,"Error encoding or decoding scanning mode flags.",}, {513 ,"Error encoding or decoding the j direction increment.",}, {514 ,"Error encoding or decoding J,K,M pentagonal resolution parameters.",}, {515 ,"Error encoding or decoding representation type or mode.",}, {517 ,"Error encoding or decoding latitude or longitude of southern pole.",}, {518 ,"Error encoding or decoding angle of rotation.",}, {519 ,"Error encoding or decoding latitude or longitude of pole of stretching.",}, {520 ,"Error encoding or decoding stretching factor.",}, {521 ,"Error encoding or decoding vertical coordinate parameters.",}, {522 ,"Error encoding or decoding list of numbers of points.",}, {523 ,"Error encoding or decoding number of points along X or Y axis.",}, {524 ,"Error encoding or decoding X or Y axis grid lengths.",}, {525 ,"Error encoding or decoding the projection centre flag.",}, {526 ,"Error encoding or decoding the latitude or longitude of sub-satellite point.",}, {527 ,"Error encoding or decoding the diameter of the Earth in X or Y direction.",}, {528 ,"Error encoding or decoding the X or Y coordinate of the sub-satellite point.",}, {529 ,"Error encoding or decoding the orientation of the grid or camera angle.",}, {530 ,"Error encoding or decoding the X or Y coordinate of the origin of sector.",}, {531 ,"Error inserting/extracting Latin1 or Latin2 of secants points.",}, {532 ,"Error inserting/extracting latitude/longitude of southern pole.",}, {540 ,"Error encoding gaussian grid section 2.",}, {541 ,"Error decoding gaussian grid section 2.",}, {542 ,"Error encoding lat/long grid section 2.",}, {543 ,"Error decoding lat/long grid section .",}, {544 ,"Error encoding spectral section 2.",}, {545 ,"Error decoding spectral section 2.",}, {546 ,"Error encoding space view section 2, or error encoding ECMWF ocean section 2.",}, {547 ,"Error decoding space view section 2, or error decoding ECMWF ocean section 2.",}, {598 ,"Representation type not handled.",}, {599 ,"Error found when checking values for section 2 against expected valid GRIB values.",}, {601 ,"Error encoding or decoding length of section 3.",}, {602 ,"Error encoding or decoding the number of unused bits at the end of section 3.",}, {603 ,"Error encoding or decoding the bitmap reference table.",}, {604 ,"Error encoding or decoding the primary bitmap.",}, {605 ,"Unable to convert quasi-regular (reduced) gaussian grid with a bitmap.",}, {699 ,"Error found when checking values for section 3 against expected valid GRIB values.",}, {701 ,"Error encoding or decoding the length of section 4.",}, {703 ,"Second-order packing implies additional flags.",}, {704 ,"Function A or B invalid for second-order packed field (grid-point).",}, {705 ,"Only simple packing handled for grid data.",}, {706 ,"Error in decoding the section 4 flag field.",}, {707 ,"Error encoding or decoding the field scale factor.",}, {708 ,"Error encoding or decoding the field reference value.",}, {709 ,"Error encoding or decoding the number of bits per data value.",}, {710 ,"Output array too small.",}, {711 ,"Error encoding or decoding a REAL coefficient.",}, {712 ,"Error encoding or decoding data values.",}, {713 ,"Error encoding or decoding the 'flag and unused bits' field.",}, {714 ,"Decoding option is 'X' and number of values is illegal.",}, {715 ,"Decoding option is 'X' and scanning mode is not north to south and west to east.",}, {716 ,"Decoding option is 'X' and field is not a gaussian or latitude/longitude grid.",}, {717 ,"Decoding option is 'X' and a bitmap is encountered.",}, {720 ,"Error encoding or decoding the octet number at which packed data begins.",}, {721 ,"Error encoding or decoding the 'extended flag' field.",}, {722 ,"Error encoding or decoding the first or second dimension of a matrix.",}, {723 ,"Error encoding or decoding six fields, from first dimension coordinate value onwards.",}, {724 ,"Error encoding or decoding first or second dimension coefficients.",}, {725 ,"Error encoding a secondary bitmap.",}, {726 ,"Output array is not big enough to allow expansion using bitmap(s), or constant (0-bit) field without section 2 not supported.",}, {727 ,"Invalid predefined bitmap number given, or constant (0-bit) spectral field not supported.",}, {728 ,"Open of predefined bitmap failed, or constant (0-bit) field function G or B not supported.",}, {729 ,"Problem reading a predefined bitmap size, or output array is not big enough to allow expansion using bitmaps.",}, {730 ,"Problem allocating memory for a predefined bitmap.",}, {731 ,"Problem reading the number of non-missing points in a predefined bitmap.",}, {732 ,"Problem reading a predefined bitmap.",}, {733 ,"Problem closing a predefined bitmap file.",}, {734 ,"Error inserting padding zeroes at end of section 4.",}, {796 ,"Error extracting pointer to packed data.",}, {797 ,"Error extracting sub-set truncation or constant (0-bit) field function G or B not supported.",}, {798 ,"Function is X but no section 2 included.",}, {799 ,"Error found when checking values for section 4 against expected valid GRIB values.",}, {800 ,"Error encoding or decoding the 7777 group.",}, {801 ,"Error inserting or extracting the 7777 group.",}, {802 ,"Error encoding or decoding the length of the GRIB message.",}, {805 ,"End of message 7777 group not found.",}, {806 ,"Error in decoding a primary or secondary bitmap.",}, {807 ,"Inconsistent values specified for complex packing in KSEC2(6) and KSEC4(4).",}, {808 ,"Error converting quasi-regular (reduced) gaussian grid to regular.",}, {809 ,"Error padding unused part of GRIB to zero.",}, {810 ,"Error inserting dummy zero, or error padding GRIB to multiple of 120.",}, {811 ,"Unable to handle secondary bitmaps for J option.",}, }; extern double grib_power(long,long); double bucket_size(grib_handle* g) { long bits;long binary_scale_factor;long decimal_scale_factor; unsigned long zero = 0; unsigned long one = 1; double z,o,s,d,reference_value; int err=0; err=grib_get_double(g,"referenceValue",&reference_value); if (err!=GRIB_SUCCESS) return 0.0; grib_get_long(g,"binaryScaleFactor",&binary_scale_factor); if (err!=GRIB_SUCCESS) return 0.0; grib_get_long(g,"numberOfBitsContainingEachPackedValue",&bits); if (err!=GRIB_SUCCESS) return 0.0; grib_get_long(g,"decimalScaleFactor",&decimal_scale_factor); if (err!=GRIB_SUCCESS) return 0.0; s = grib_power(binary_scale_factor,2); d = grib_power(-decimal_scale_factor,10) ; z = (double) (((zero*s)+reference_value)*d); o = (double) (((one *s)+reference_value)*d); return o-z; } static double err(double A, double B) { double relativeError; if(fabs(A) <= maxAbsoluteError || fabs(B) <= maxAbsoluteError) relativeError = fabs(A-B); else if (fabs(B) > fabs(A)) relativeError = fabs((A - B) / B); else relativeError = fabs((A - B) / A); return relativeError; } #define NUMBER(a) (sizeof(a)/sizeof(a[0])) extern void grsrnd_(long*); extern void gribex_(long*, long*, long*, double*, long*, double*, long*, double*, long*, char *, long*, long*, const char*, long*,long ); static int cgribex( long miss, long ksec0[],long ksec1[], long ksec2[],double rsec2[], long ksec3[],double rsec3[], long ksec4[],double rsec4[], long sec4len, char *buffer, long *buflen, char *op) { long ret = 1; long in_bufsize = (*buflen) / sizeof(long); long out_bufsize = in_bufsize; long round = 0; rsec3[1] = miss, ksec3[1] = miss, grsrnd_(&round); /* no rounding !!! */ gribex_(ksec0, ksec1, ksec2, rsec2, ksec3, rsec3, ksec4, rsec4, &sec4len, buffer, &in_bufsize, &out_bufsize, op, &ret,1 ); *buflen = out_bufsize * sizeof(long); return ret; } static double proc_cpu() { struct rusage rup; if(getrusage(RUSAGE_SELF,&rup) != -1) { return (rup.ru_utime.tv_sec + rup.ru_utime.tv_usec / 1000000.0 + rup.ru_stime.tv_sec + rup.ru_stime.tv_usec / 1000000.0); } return clock()/(double)CLOCKS_PER_SEC; } int compare_values(double* v1,double* v2,int count,const char *what) { int i,j=0; double maxError = 0; for(i = 0; i < count ;i++) { double e = err(v1[i],v2[i]); if(e > maxError) { j = i; maxError = e; } } fprintf(stdout,"value no %d out of %d %s v1=%g v2=%g abs=%.30e rel=%g\n",j,count,what,v1[j], v2[j],v1[j]-v2[j],maxError); if(maxError > maxRelativeError) { fprintf(stderr,"value no %d out of %d %s v1=%g v2=%g abs=%g rel=%g\n",j,count,what,v1[j], v2[j],v1[j]-v2[j],maxError); return 1; } return 0; } void get_max_relative_error(double* v1,double* v2,int count,double* maxError,int* id) { int i; *id=0; *maxError = 0; for(i = 0; i < count ;i++) { double e = err(v1[i],v2[i]); if(e > *maxError) { *id = i; *maxError = e; } } } int compare_with_gribex(grib_handle* g) { long ed = 0; int err = 0; long missing = 0; long is_int; long miss; int i; size_t count; long sec4len; long len; double now,gtime,atime; size_t size; const void *buffer; long ksec0[ISECTION_0]; long ksec1[ISECTION_1]; long ksec2[ISECTION_2]; long ksec3[ISECTION_3]; long ksec4[ISECTION_4]; double rsec2[RSECTION_2]; double rsec3[RSECTION_3]; double *rsec4; double *values; if((err = grib_get_long(g,"editionNumber",&ed))) { fprintf(stderr,"Cannot get editionNumber %d\n",err); return err; } if(ed > 1) return 0; if((err = grib_get_long(g,"missingValue",&miss))) { fprintf(stderr,"Cannot get missingValue %d\n",err); return err; } if((err = grib_get_long(g,"integerPointValues",&is_int))) { fprintf(stderr,"Cannot get integerPointValues %d\n",err); return err; } if((err = grib_get_size(g,"values",&count))) { fprintf(stderr,"Cannot get number of values %d\n",err); return err; } printf("Number of values %ld\n",(long)count); /* printf("number of values is %d\n",count); */ rsec4 = (double*)malloc(sizeof(double)*(count+100000)); values = (double*)malloc(sizeof(double)*(count+100000)); sec4len = count +100000; if((err = grib_get_message(g,&buffer,&size))) { fprintf(stderr,"Cannot get message %d\n",err); goto foo; } len = size; now = proc_cpu(); err = cgribex( miss, ksec0,ksec1,ksec2,rsec2,ksec3,rsec3, ksec4,rsec4,sec4len, (char*)buffer,&len,"D"); gtime = proc_cpu() - now; if(err > 0) { char *e = "unknown"; for(i = 0; i < NUMBER(gribex_errors); i++) if(gribex_errors[i].code == err) e = gribex_errors[i].message; fprintf(stderr,"GRIBEX error %d: %s\n",err,e); fprintf(stderr,"Count api=%ld gribex=%ld\n", (long)count,(long)ksec4[0]); goto foo; } if(ksec4[0] == 0){ fprintf(stderr,"No GRIB processed - no values in KSEC4\n"); goto foo; } now = proc_cpu(); if((err = grib_get_double_array(g,"values",values,&count))) { fprintf(stderr,"Cannot get values %d (%s) count=%ld\n",err,grib_get_error_message(err),(long)count); return err; } atime = proc_cpu() - now; if( ksec4[0] < 0) { missing = 1; ksec4[0] = -ksec4[0]; } if(count != ksec4[0]) { fprintf(stderr,"Count mismatch api=%ld gribex=%ld\n", (long)count,ksec4[0]); err = 1; goto foo; } #if 0 printf("Values: %d gribex: %g api:%g\n",count,gtime,atime); if(atime > gtime*1.1) { fprintf(stderr,"Too slow: gribex: %g api:%g\n",gtime,atime); exit(1); } #endif if(missing) { for(i = 0; i < count ;i++) { if(values[i] != miss) { fprintf(stderr,"Missing field: value %g not set to %ld\n",values[i],miss); err = 1; goto foo; } } } else if(is_int) { long long* lval = (long long*)rsec4; for(i = 0; i < count ;i++) rsec4[i] = lval[i]; if(compare_values(rsec4,values,count,"gribex")) { err = 1; goto foo; } } else { if(compare_values(rsec4,values,count,"gribex")) { err = 1; goto foo; } } foo: free(rsec4); free(values); return err; } int compare_two_gribex(grib_handle* g) { long ed = 0; int err = 0; long is_int; long miss; int i; size_t count; long sec4len; long len; size_t size; const void *buffer; long ksec0[ISECTION_0]; long ksec1[ISECTION_1]; long ksec2[ISECTION_2]; long ksec3[ISECTION_3]; long ksec4[ISECTION_4]; double rsec2[RSECTION_2]; double rsec3[RSECTION_3]; double *rsec4_1 =0 ; double *rsec4_2 =0 ; double *rsec4_3 =0 ; char *b; if((err = grib_get_long(g,"editionNumber",&ed))) { fprintf(stderr,"Cannot get editionNumber %d\n",err); return err; } if(ed > 1) return 0; if((err = grib_get_long(g,"missingValue",&miss))) { fprintf(stderr,"Cannot get missingValue %d\n",err); return err; } if((err = grib_get_long(g,"integerPointValues",&is_int))) { fprintf(stderr,"Cannot get integerPointValues %d\n",err); return err; } if(is_int) return 0; if((err = grib_get_size(g,"values",&count))) { fprintf(stderr,"Cannot get number of values %d\n",err); return err; } printf("Number of values %ld\n",(long)count); /* printf("number of values is %d\n",count); */ rsec4_1 = (double*)malloc(sizeof(double)*(count+100000)); rsec4_2 = (double*)malloc(sizeof(double)*(count+100000)); rsec4_3 = (double*)malloc(sizeof(double)*(count+100000)); b = (char*)malloc(sizeof(char)*(count+100000)); sec4len = count +100000; if((err = grib_get_message(g,&buffer,&size))) { fprintf(stderr,"Cannot get message %d\n",err); goto foo; } len = size; err = cgribex( miss, ksec0,ksec1,ksec2,rsec2,ksec3,rsec3, ksec4,rsec4_1,sec4len, (char*)buffer,&len,"D"); if(err > 0) { char *e = "unknown"; for(i = 0; i < NUMBER(gribex_errors); i++) if(gribex_errors[i].code == err) e = gribex_errors[i].message; fprintf(stderr,"GRIBEX error %d: %s\n",err,e); fprintf(stderr,"Count api=%ld gribex=%ld\n", (long)count,ksec4[0]); goto foo; } memcpy(rsec4_3,rsec4_1,sizeof(double)*(count+100000)); if(ksec4[0] == 0){ fprintf(stderr,"No GRIB processed - no values in KSEC4\n"); goto foo; } if( ksec4[0] < 0) goto foo; if(count != ksec4[0]) { fprintf(stderr,"Count mismatch api=%ld gribex=%ld\n", (long)count,ksec4[0]); err = 1; goto foo; } len = count +100000; err = cgribex( miss, ksec0,ksec1,ksec2,rsec2,ksec3,rsec3, ksec4,rsec4_1,sec4len, (char*)b,&len,"C"); if(err > 0) { char *e = "unknown"; for(i = 0; i < NUMBER(gribex_errors); i++) if(gribex_errors[i].code == err) e = gribex_errors[i].message; fprintf(stderr,"GRIBEX error %d: %s\n",err,e); fprintf(stderr,"Count api=%ld gribex=%ld\n", (long)count,ksec4[0]); goto foo; } sec4len = count +100000; err = cgribex( miss, ksec0,ksec1,ksec2,rsec2,ksec3,rsec3, ksec4,rsec4_2,sec4len, (char*)buffer,&len,"D"); if(err > 0) { char *e = "unknown"; for(i = 0; i < NUMBER(gribex_errors); i++) if(gribex_errors[i].code == err) e = gribex_errors[i].message; fprintf(stderr,"GRIBEX error %d: %s\n",err,e); fprintf(stderr,"Count api=%ld gribex=%ld\n", (long)count,ksec4[0]); goto foo; } if(compare_values(rsec4_3,rsec4_2,count,"gribex2")) { err = 1; goto foo; } foo: free(rsec4_1); free(rsec4_2); free(rsec4_3); free(b); return err; } int compare_lists(grib_handle* g,grib_handle *h,const char* what) { int err; size_t count = 0; size_t ncount = 0; double *values1; double *values2; int e1 = grib_get_size(g,what,&count); int e2 = grib_get_size(h,what,&ncount); if(e1 == GRIB_NOT_FOUND && e2 == GRIB_NOT_FOUND) return 0; if(e1) { fprintf(stderr,"Cannot get number of '%s' %d\n",what,e1); return 1; } if(e2) { fprintf(stderr,"Cannot get number of '%s' %d\n",what,e2); return 1; } if(ncount != count) { fprintf(stderr,"'%s' count mismatch %ld -> %ld\n",what,(long)count,(long)ncount); return 1; } values1 = (double*)malloc(sizeof(double)*(count)); values2 = (double*)malloc(sizeof(double)*(count)); ncount = count; if((err = grib_get_double_array(g,what,values1,&count))) { fprintf(stderr,"Cannot get %s1 %d (%s) count=%ld\n",what,err,grib_get_error_message(err),(long)count); goto foo; } if(ncount != count) { fprintf(stderr,"%s count changes %ld -> %ld\n",what,(long)count,(long)ncount); err = 1; goto foo; } ncount = count; if((err = grib_get_double_array(h,what,values2,&count))) { fprintf(stderr,"Cannot get %s2 %d (%s) count=%ld\n",what,err,grib_get_error_message(err),(long)ncount); goto foo; } if(ncount != count) { fprintf(stderr,"%s count changes %ld -> %ld\n",what,(long)count,(long)ncount); err = 1; goto foo; } err = compare_values(values1,values2,count,what); foo: free(values2); free(values1); return err; } int compare_with_handle(grib_handle* g,grib_handle *h) { double s1,d1,r1,l1; double s2,d2,r2,l2; l1 = s1 = r1 = d1 = 1e-100; grib_get_double(g,"binaryScaleFactor",&s1); grib_get_double(g,"referenceValue",&r1); grib_get_double(g,"decimalScaleFactor",&d1); grib_get_double(g,"laplacianOperator",&l1); l2 = s2 = r2 = d2 = 1e-100; grib_get_double(h,"binaryScaleFactor",&s2); grib_get_double(h,"referenceValue",&r2); grib_get_double(h,"decimalScaleFactor",&d2); grib_get_double(g,"laplacianOperator",&l2); if(s1 != s2) { fprintf(stderr, "Scale factor 1 = %g " "Scale factor 2 = %g diff=%g" "\n" ,s1,s2,s1-s2); } if(d1 != d2) { fprintf(stderr, "Decimal Scale factor 1 = %g " "Decimal Scale factor 2 = %g diff=%g" "\n" ,d1,d2,d1-d2); } if(r1 != r2) { fprintf(stderr, "Ref. value 1 = %g " "Ref. value 2 = %g diff=%g" "\n" ,r1,r2,r1-r2); } if(l1 != l2) { fprintf(stderr, "Laplacian 1 = %g " "Laplacian 2 = %g diff=%g" "\n" ,l1,l2,l1-l2); } if(compare_lists(g,h,"values")) return 1; if(compare_lists(g,h,"pv")) return 1; if(compare_lists(g,h,"pl")) return 1; return 0; } int save_message(grib_handle* h,const char* path) { FILE* out = fopen(path,"w"); const void *p; size_t len; if(!out) { perror(path); return -1; } printf("save to %s\n",path); grib_get_message(h,&p,&len); if(fwrite(p,1,len,out) != len) { perror(path); fclose(out); return -1; } if(fclose(out)) { perror(path); return -1; } return 0; } const char* problem(const char* path,const char* name,int field) { static char buf[1024]; const char *p = path; const char *q = path; while(*p) { if(*p == '/') q=p+1; p++; } sprintf(buf,"/vol/dataserv/grib_samples/problems/%s_%s_%d.grib",q,name,field); return buf; }