core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

ECC-625: part1

This commit is contained in:
Shahram Najm 2018-02-05 16:05:31 +00:00
parent c6874c4cf7
commit 6eb9ef3b45
2 changed files with 124 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

128
src/grib_util.c
View File

@ -32,6 +32,8 @@
#endif
#endif
typedef enum {eROUND_ANGLE_UP, eROUND_ANGLE_DOWN} RoundingPolicy;
static void set_total_length(unsigned char* buffer,long *section_length,long *section_offset,int edition,size_t totalLength)
{
long off;
@ -349,16 +351,16 @@ static void print_values(grib_context* c, const grib_util_grid_spec2* spec,
{
switch(values[i].type)
{
case GRIB_TYPE_LONG: printf("ECCODES DEBUG grib_util: => %s = %ld;\n"
case GRIB_TYPE_LONG: printf("ECCODES DEBUG grib_util: => %s = %ld;\n"
,values[i].name,(long)values[i].long_value); break;
case GRIB_TYPE_DOUBLE: printf("ECCODES DEBUG grib_util: => %s = %.16e;\n"
case GRIB_TYPE_DOUBLE: printf("ECCODES DEBUG grib_util: => %s = %.16e;\n"
,values[i].name,values[i].double_value); break;
case GRIB_TYPE_STRING: printf("ECCODES DEBUG grib_util: => %s = \"%s\";\n"
case GRIB_TYPE_STRING: printf("ECCODES DEBUG grib_util: => %s = \"%s\";\n"
,values[i].name,values[i].string_value); break;
}
}
printf("ECCODES DEBUG grib_util: data_values_count=%lu;\n", (unsigned long)data_values_count);
printf("ECCODES DEBUG grib_util: data_values_count=%lu;\n", (unsigned long)data_values_count);
for (i=0; i<data_values_count; i++) {
if (i==0) v = data_values[i];
if (data_values[i] != spec->missingValue) {
@ -378,7 +380,7 @@ static void print_values(grib_context* c, const grib_util_grid_spec2* spec,
if (v > maxVal) maxVal=v;
}
}
printf("ECCODES DEBUG grib_util: data_values are CONSTANT? %d\t(min=%.16e, max=%.16e)\n",
printf("ECCODES DEBUG grib_util: data_values are CONSTANT? %d\t(min=%.16e, max=%.16e)\n",
isConstant, minVal, maxVal);
#if 0
@ -406,6 +408,7 @@ static int DBL_EQUAL(double d1, double d2, double tolerance)
return fabs(d1-d2) < tolerance;
}
*/
#if 0
/* Returns a boolean: 1 if angle can be encoded, 0 otherwise */
static int grib1_angle_can_be_encoded(const double angle)
{
@ -417,6 +420,80 @@ static int grib1_angle_can_be_encoded(const double angle)
if (angle == rounded) return 1;
return 0; /* sub millidegree. Cannot be encoded in grib1 */
}
#endif
/* Returns a boolean: 1 if angle can be encoded, 0 otherwise */
static int angle_can_be_encoded(const double angle, const double angular_precision)
{
const double angle_expanded = angle * angular_precision;
Assert(angular_precision>0);
double rounded = (int)(angle_expanded+0.5)/angular_precision;
if (angle<0) {
rounded = (int)(angle_expanded-0.5)/angular_precision;
}
if (angle == rounded) return 1;
/*printf(" ......... angle cannot be encoded: %.10e\n", angle);*/
return 0; /* Cannot be encoded */
}
static double adjust_angle(const double angle, const RoundingPolicy policy, const double angular_precision)
{
double result = 0;
Assert(angular_precision > 0);
result = angle * angular_precision;
if (policy == eROUND_ANGLE_UP)
result = round(result+0.5);
else
result = round(result-0.5);
result = result / angular_precision;
return result;
}
/* Search key=value array for:
* latitudeOfFirstGridPointInDegrees
* longitudeOfFirstGridPointInDegrees
* latitudeOfLastGridPointInDegrees
* longitudeOfLastGridPointInDegrees
* and change their values to expand the bounding box ??
*/
static int expand_bounding_box(grib_handle* h, grib_values *values, const size_t count)
{
int ret = GRIB_SUCCESS;
size_t i=0;
double new_angle = 0;
RoundingPolicy roundingPolicy = eROUND_ANGLE_UP;
long angular_precision = 0; /* e.g. 1e3 for grib1 and 1e6 for grib2 */
if((ret = grib_get_long(h, "angularPrecision", &angular_precision)) != 0)
return ret;
for(i=0; i<count; i++) {
int is_angle = 0;
if (strcmp(values[i].name, "longitudeOfFirstGridPointInDegrees")==0) {
roundingPolicy = eROUND_ANGLE_DOWN;
is_angle = 1;
}
else if (strcmp(values[i].name, "longitudeOfLastGridPointInDegrees")==0) {
roundingPolicy = eROUND_ANGLE_UP;
is_angle = 1;
}
else if (strcmp(values[i].name, "latitudeOfFirstGridPointInDegrees")==0) {
roundingPolicy = eROUND_ANGLE_UP;
is_angle = 1;
}
else if (strcmp(values[i].name, "latitudeOfLastGridPointInDegrees")==0) {
roundingPolicy = eROUND_ANGLE_DOWN;
is_angle = 1;
}
if (is_angle && !angle_can_be_encoded(values[i].double_value, angular_precision)) {
new_angle = adjust_angle(values[i].double_value, roundingPolicy, angular_precision);
if (h->context->debug)
printf("ECCODES DEBUG grib_util expand_bounding_box %s: old=%.15e new=%.15e\n", values[i].name, values[i].double_value, new_angle);
values[i].double_value = new_angle;
}
}
return ret;
}
/* Returns a boolean: 1 if angle is too small, 0 otherwise */
/*static int angle_too_small(const double angle, const double angular_precision)
@ -701,6 +778,7 @@ grib_handle* grib_util_set_spec2(grib_handle* h,
size_t slen=17;
int grib1_high_resolution_fix = 0; /* boolean: See GRIB-863 */
int global_grid = 0; /* boolean */
int expandBoundingBox = 0;
static grib_util_packing_spec default_packing_spec = {0, };
Assert(h);
@ -927,8 +1005,9 @@ grib_handle* grib_util_set_spec2(grib_handle* h,
SET_LONG_VALUE ("ijDirectionIncrementGiven", 1);
if (editionNumber == 1) {
/* GRIB-863: GRIB1 cannot represent increments less than a millidegree */
if (!grib1_angle_can_be_encoded(spec->iDirectionIncrementInDegrees) ||
!grib1_angle_can_be_encoded(spec->jDirectionIncrementInDegrees))
const double grib1_angular_precision = 1000;
if (!angle_can_be_encoded(spec->iDirectionIncrementInDegrees, grib1_angular_precision) ||
!angle_can_be_encoded(spec->jDirectionIncrementInDegrees, grib1_angular_precision))
{
grib1_high_resolution_fix = 1;
/* Set flag to compute the increments */
@ -1172,13 +1251,24 @@ grib_handle* grib_util_set_spec2(grib_handle* h,
if(packing_spec->extra_settings_count) {
for(i = 0; i < packing_spec->extra_settings_count; i++) {
Assert(count < 1024);
values[count++] = packing_spec->extra_settings[i];
if (strcmp(packing_spec->extra_settings[i].name, "global")==0 &&
packing_spec->extra_settings[i].long_value == 1)
if (strcmp(packing_spec->extra_settings[i].name, "expandBoundingBox")==0)
{
/* GRIB-922: Request is for a global grid. Setting this key will
* calculate the lat/lon values. So the spec's lat/lon can be ignored */
global_grid = 1;
if (packing_spec->extra_settings[i].long_value == 1) {
/* ECC-625: Request is for expansion of bounding box (sub-area).
* This is also called the "snap-out" policy */
expandBoundingBox = 1;
}
}
else
{
values[count++] = packing_spec->extra_settings[i];
if (strcmp(packing_spec->extra_settings[i].name, "global")==0 &&
packing_spec->extra_settings[i].long_value == 1)
{
/* GRIB-922: Request is for a global grid. Setting this key will
* calculate the lat/lon values. So the spec's lat/lon can be ignored */
global_grid = 1;
}
}
}
}
@ -1216,6 +1306,17 @@ grib_handle* grib_util_set_spec2(grib_handle* h,
if (h->context->debug==-1) {
printf("ECCODES DEBUG grib_util: global_grid = %d\n", global_grid);
printf("ECCODES DEBUG grib_util: expandBoundingBox = %d\n", expandBoundingBox);
print_values(h->context,spec,data_values,data_values_count,values,count);
}
/* Apply adjustments to bounding box if needed */
if (expandBoundingBox) {
if ((*err=expand_bounding_box(outh, values, count)) != 0)
{
fprintf(stderr,"SET_GRID_DATA_DESCRIPTION: Cannot expand bounding box: %s\n",grib_get_error_message(*err));
goto cleanup;
}
print_values(h->context,spec,data_values,data_values_count,values,count);
}
@ -1449,6 +1550,7 @@ int grib_moments(grib_handle* h,double east,double north,double west,double sout
grib_context_free(c,lat);
grib_context_free(c,lon);
grib_context_free(c,values);
(void)mass;
return ret;
}

10
tests/grib_util_set_spec.c
View File

@ -72,7 +72,7 @@ static void test_reduced_gg(int remove_local_def, int edition, const char* packi
outlen = inlen;
spec.iDirectionIncrementInDegrees = 1.5;
spec.jDirectionIncrementInDegrees = 1.5;
spec.latitudeOfFirstGridPointInDegrees = 87.863799;
spec.latitudeOfFirstGridPointInDegrees = 87.8637991;
spec.longitudeOfFirstGridPointInDegrees = 0.0;
spec.latitudeOfLastGridPointInDegrees = -spec.latitudeOfFirstGridPointInDegrees;
spec.longitudeOfLastGridPointInDegrees = 357.187500;
@ -83,6 +83,13 @@ static void test_reduced_gg(int remove_local_def, int edition, const char* packi
packing_spec.bitsPerValue = 24;
packing_spec.accuracy=GRIB_UTIL_ACCURACY_USE_PROVIDED_BITS_PER_VALUES;
packing_spec.packing=GRIB_UTIL_PACKING_USE_PROVIDED;
/*Extra settings
packing_spec.extra_settings_count++;
packing_spec.extra_settings[0].type = GRIB_TYPE_LONG;
packing_spec.extra_settings[0].name = "expandBoundingBox";
packing_spec.extra_settings[0].long_value = 1;
*/
finalh = grib_util_set_spec(
handle,
@ -104,6 +111,7 @@ static void test_reduced_gg(int remove_local_def, int edition, const char* packi
grib_handle_delete(finalh);
fclose(in);
fclose(out);
/*printf("ALL OK: %s\n", __func__);*/
}
static void test_regular_ll(int remove_local_def, int edition, const char* packingType,