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

2252 lines
82 KiB
C
Raw Blame History

/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include <float.h>
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;
switch (edition) {
case 1:
if (totalLength < 0x800000) {
off = 32;
grib_encode_unsigned_long(buffer, (unsigned long)totalLength, &off, 24);
}
else {
long s4len, t120;
totalLength -= 4;
t120 = (totalLength + 119) / 120;
s4len = t120 * 120 - totalLength;
totalLength = 0x800000 | t120;
off = 32;
grib_encode_unsigned_long(buffer, (unsigned long)totalLength, &off, 24);
off = section_offset[4] * 8;
grib_encode_unsigned_long(buffer, (unsigned long)s4len, &off, 24);
}
break;
case 2:
off = 64;
grib_encode_unsigned_long(buffer, (unsigned long)totalLength, &off, 64);
break;
}
}
static grib_handle* grib_sections_copy_internal(grib_handle* hfrom, grib_handle* hto, int sections[], int* err)
{
int i;
size_t totalLength = 0;
unsigned char* buffer;
unsigned char* p;
long edition = 0;
long section_length[MAX_NUM_SECTIONS] = {0,};
long section_offset[MAX_NUM_SECTIONS] = {0,};
long off = 0;
grib_handle* h;
char section_length_str[64] = "section0Length";
char section_offset_str[64] = "offsetSection0";
long length, offset;
*err = grib_get_long(hfrom, "edition", &edition);
if (*err)
return NULL;
for (i = 0; i <= hfrom->sections_count; i++) {
if (sections[i]) {
h = hfrom;
}
else {
h = hto;
}
sprintf(section_length_str, "section%dLength", i);
if (grib_get_long(h, section_length_str, &length))
continue;
section_length[i] = length;
sprintf(section_offset_str, "offsetSection%d", i);
if (grib_get_long(h, section_offset_str, &offset))
continue;
section_offset[i] = offset;
totalLength += section_length[i];
}
buffer = (unsigned char*)grib_context_malloc_clear(hfrom->context, totalLength * sizeof(char));
p = buffer;
off = 0;
for (i = 0; i <= hfrom->sections_count; i++) {
grib_handle* hand = NULL;
if (sections[i])
hand = hfrom;
else
hand = hto;
p = (unsigned char*)memcpy(p, hand->buffer->data + section_offset[i], section_length[i]);
section_offset[i] = off;
off += section_length[i];
p += section_length[i];
}
/* copy section 3 present flag*/
if (edition == 1) {
const void* buffer_to = NULL;
size_t size_to = 0;
grib_get_message(hto, &buffer_to, &size_to);
memcpy(buffer + 15, ((unsigned char*)buffer_to) + 15, 1);
}
set_total_length(buffer, section_length, section_offset, edition, totalLength);
h = grib_handle_new_from_message(hfrom->context, buffer, totalLength);
/*to allow free of buffer*/
h->buffer->property = GRIB_MY_BUFFER;
switch (edition) {
case 1:
if (sections[1] && sections[2])
break;
if (sections[1]) {
long PVPresent;
grib_get_long(hfrom, "PVPresent", &PVPresent);
if (PVPresent) {
double* pv;
long numberOfVerticalCoordinateValues;
size_t size = 0;
grib_get_long(hfrom, "numberOfVerticalCoordinateValues", &numberOfVerticalCoordinateValues);
size = numberOfVerticalCoordinateValues;
pv = (double*)grib_context_malloc_clear(hfrom->context, numberOfVerticalCoordinateValues * sizeof(double));
grib_get_double_array(hfrom, "pv", pv, &size);
grib_set_long(h, "PVPresent", 1);
grib_set_double_array(h, "pv", pv, size);
grib_context_free(hfrom->context, pv);
}
else {
grib_set_long(h, "PVPresent", 0);
}
}
if (sections[2]) {
long PVPresent;
grib_get_long(hto, "PVPresent", &PVPresent);
if (PVPresent) {
double* pv;
long numberOfVerticalCoordinateValues;
size_t size = 0;
grib_get_long(hto, "numberOfVerticalCoordinateValues", &numberOfVerticalCoordinateValues);
size = numberOfVerticalCoordinateValues;
pv = (double*)grib_context_malloc_clear(hto->context, numberOfVerticalCoordinateValues * sizeof(double));
grib_get_double_array(hto, "pv", pv, &size);
grib_set_long(h, "PVPresent", 1);
grib_set_double_array(h, "pv", pv, size);
grib_context_free(hto->context, pv);
}
else {
grib_set_long(h, "PVPresent", 0);
}
}
break;
case 2:
if (sections[1]) {
long discipline;
grib_get_long(hfrom, "discipline", &discipline);
grib_set_long(h, "discipline", discipline);
}
break;
}
return h;
}
grib_handle* grib_util_sections_copy(grib_handle* hfrom, grib_handle* hto, int what, int* err)
{
long edition_from = 0;
long edition_to = 0;
long localDefinitionNumber = -1;
int sections_to_copy[MAX_NUM_SECTIONS] = {0,};
*err = grib_get_long(hfrom, "edition", &edition_from);
if (*err)
return NULL;
*err = grib_get_long(hto, "edition", &edition_to);
if (*err)
return NULL;
if (edition_to != 1 && edition_to != 2) {
*err = GRIB_NOT_IMPLEMENTED;
return NULL;
}
if (edition_from != edition_to) {
*err = GRIB_DIFFERENT_EDITION;
return NULL;
}
if (what & GRIB_SECTION_GRID) {
switch (edition_from) {
case 1:
sections_to_copy[2] = 1;
break;
case 2:
sections_to_copy[3] = 1;
break;
}
}
if (what & GRIB_SECTION_DATA) {
switch (edition_from) {
case 1:
sections_to_copy[3] = 1;
sections_to_copy[4] = 1;
break;
case 2:
sections_to_copy[5] = 1;
sections_to_copy[6] = 1;
sections_to_copy[7] = 1;
break;
}
}
if (what & GRIB_SECTION_LOCAL) {
switch (edition_from) {
case 1:
sections_to_copy[1] = 1;
break;
case 2:
sections_to_copy[2] = 1;
break;
}
}
if (what & GRIB_SECTION_PRODUCT) {
switch (edition_from) {
case 1:
grib_get_long(hfrom, "localDefinitionNumber", &localDefinitionNumber);
if (localDefinitionNumber == 13) {
sections_to_copy[4] = 1;
}
sections_to_copy[1] = 1;
break;
case 2:
sections_to_copy[1] = 1;
sections_to_copy[4] = 1;
break;
}
}
if (what & GRIB_SECTION_BITMAP) {
switch (edition_from) {
case 1:
sections_to_copy[3] = 1;
break;
case 2:
sections_to_copy[6] = 1;
break;
}
}
return grib_sections_copy_internal(hfrom, hto, sections_to_copy, err);
}
static grib_trie* init_list(const char* name);
static grib_trie* param_id_list = NULL;
static grib_trie* mars_param_list = NULL;
/* TODO thread safe */
grib_string_list* grib_util_get_param_id(const char* mars_param)
{
if (!mars_param_list && (mars_param_list = init_list("mars_param.table")) == NULL)
return NULL;
return (grib_string_list*)grib_trie_get(mars_param_list, mars_param);
}
grib_string_list* grib_util_get_mars_param(const char* param_id)
{
if (!param_id_list && (param_id_list = init_list("param_id.table")) == NULL)
return NULL;
return (grib_string_list*)grib_trie_get(param_id_list, param_id);
}
static grib_trie* init_list(const char* name)
{
char* full_path = 0;
FILE* fh;
char s[101];
char param[101];
grib_string_list* list = 0;
grib_string_list* next = 0;
grib_trie* trie_list;
grib_context* c = grib_context_get_default();
full_path = grib_context_full_defs_path(c, name);
fh = codes_fopen(full_path, "r");
if (!fh) {
grib_context_log(c, GRIB_LOG_PERROR, "unable to read %s", full_path);
return NULL;
}
list = (grib_string_list*)grib_context_malloc_clear(c, sizeof(grib_string_list));
trie_list = grib_trie_new(c);
if (fscanf(fh, "%100s", param) == EOF) {
fclose(fh);
return NULL;
}
while (fscanf(fh, "%100s", s) != EOF) {
if (!strcmp(s, "|")) {
grib_trie_insert(trie_list, param, list);
if (fscanf(fh, "%100s", param) == EOF) {
fclose(fh);
return trie_list;
}
list = NULL;
}
else {
if (!list) {
list = (grib_string_list*)grib_context_malloc_clear(c, sizeof(grib_string_list));
list->value = grib_context_strdup(c, s);
}
else {
next = list;
while (next->next)
next = next->next;
next->next = (grib_string_list*)grib_context_malloc_clear(c, sizeof(grib_string_list));
next->next->value = grib_context_strdup(c, s);
}
}
}
fclose(fh);
return 0;
}
/* For debugging purposes */
static void print_values(grib_context* c, const grib_util_grid_spec2* spec,
const double* data_values, const size_t data_values_count, /* the data pay load */
const grib_values* keyval_pairs, const size_t count) /* keys and their values */
{
size_t i = 0;
int isConstant = 1;
double v = 0, minVal = DBL_MAX, maxVal = -DBL_MAX;
fprintf(stderr, "ECCODES DEBUG grib_util grib_set_values: setting %lu key/value pairs\n", (unsigned long)count);
for (i = 0; i < count; i++) {
switch (keyval_pairs[i].type) {
case GRIB_TYPE_LONG:
fprintf(stderr, "ECCODES DEBUG grib_util: => %s = %ld;\n", keyval_pairs[i].name, (long)keyval_pairs[i].long_value);
break;
case GRIB_TYPE_DOUBLE:
fprintf(stderr, "ECCODES DEBUG grib_util: => %s = %.16e;\n", keyval_pairs[i].name, keyval_pairs[i].double_value);
break;
case GRIB_TYPE_STRING:
fprintf(stderr, "ECCODES DEBUG grib_util: => %s = \"%s\";\n", keyval_pairs[i].name, keyval_pairs[i].string_value);
break;
}
}
fprintf(stderr, "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) {
if (v == spec->missingValue) {
v = data_values[i];
}
else if (v != data_values[i]) {
isConstant = 0;
break;
}
}
}
for (i = 0; i < data_values_count; i++) {
v = data_values[i];
if (v != spec->missingValue) {
if (v < minVal)
minVal = v;
if (v > maxVal)
maxVal = v;
}
}
fprintf(stderr, "ECCODES DEBUG grib_util: data_values are CONSTANT? %d\t(min=%.16e, max=%.16e)\n",
isConstant, minVal, maxVal);
if (c->gribex_mode_on)
fprintf(stderr, "ECCODES DEBUG grib_util: GRIBEX mode is turned on!\n");
#if 0
if (spec->bitmapPresent) {
int missing = 0;
size_t j = 0;
double min = 1e100;
for(j = 0; j < data_values_count ; j++)
{
double d = data_values[j] - spec->missingValue;
if(d < 0) d = -d;
if(d < min) {
min = d;
}
if(data_values[j] == spec->missingValue)
missing++;
}
}
#endif
}
/*
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)
{
const double angle_milliDegrees = angle * 1000;
double rounded = (int)(angle_milliDegrees+0.5)/1000.0;
if (angle<0) {
rounded = (int)(angle_milliDegrees-0.5)/1000.0;
}
if (angle == rounded) return 1;
return 0; /* sub millidegree. Cannot be encoded in grib1 */
}
/* 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 = (long)(angle_expanded+0.5)/angular_precision;
if (angle<0) {
rounded = (long)(angle_expanded-0.5)/angular_precision;
}
if (angle == rounded) return 1;
/*printf(" ......... angle cannot be encoded: %.10e\n", angle);*/
return 0; /* Cannot be encoded */
}
#endif
/* Returns a boolean: 1 if angle can be encoded, 0 otherwise */
static int angle_can_be_encoded(grib_handle* h, const double angle)
{
int ret = 0;
int retval = 1;
grib_handle* h2 = NULL;
char sample_name[16] = {0,};
long angle_subdivisions = 0; /* e.g. 1e3 for grib1 and 1e6 for grib2 */
long edition = 0, coded = 0;
double expanded, diff;
if ((ret = grib_get_long(h, "edition", &edition)) != 0)
return ret;
if ((ret = grib_get_long(h, "angleSubdivisions", &angle_subdivisions)) != 0)
return ret;
Assert(angle_subdivisions > 0);
sprintf(sample_name, "GRIB%ld", edition);
h2 = grib_handle_new_from_samples(0, sample_name);
if ((ret = grib_set_double(h2, "latitudeOfFirstGridPointInDegrees", angle)) != 0)
return ret;
if ((ret = grib_get_long(h2, "latitudeOfFirstGridPoint", &coded)) != 0)
return ret;
grib_handle_delete(h2);
expanded = angle * angle_subdivisions;
diff = fabs(expanded - coded);
if (diff < 1.0 / angle_subdivisions)
retval = 1;
else
retval = 0;
return retval;
}
static double adjust_angle(const double angle, const RoundingPolicy policy, const double angle_subdivisions)
{
double result = 0;
Assert(angle_subdivisions > 0);
result = angle * angle_subdivisions;
if (policy == eROUND_ANGLE_UP)
result = round(result + 0.5);
else
result = round(result - 0.5);
result = result / angle_subdivisions;
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 angle_subdivisions = 0; /* e.g. 1e3 for grib1 and 1e6 for grib2 */
if ((ret = grib_get_long(h, "angleSubdivisions", &angle_subdivisions)) != 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(h, values[i].double_value)) {
new_angle = adjust_angle(values[i].double_value, roundingPolicy, angle_subdivisions);
if (h->context->debug) {
fprintf(stderr, "ECCODES DEBUG grib_util EXPAND_BOUNDING_BOX %s: old=%.15e new=%.15e (%s)\n",
values[i].name, values[i].double_value, new_angle,
(roundingPolicy == eROUND_ANGLE_UP ? "Up" : "Down"));
}
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)
{
const double a = fabs(angle);
if (a > 0 && a < angular_precision) return 1;
return 0;
}
static double normalise_angle(double angle)
{
while (angle<0) angle += 360;
while (angle>360) angle -= 360;
return angle;
}
static int check_values(const double* data_values, size_t data_values_count)
{
size_t i = 0;
for (i=0; i<data_values_count; i++) {
const double val = data_values[i];
if ( val >= DBL_MAX ||
val <= -DBL_MAX ||
isnan(val) )
{
fprintf(stderr,"GRIB_UTIL_SET_SPEC: Invalid data value: i=%lu, val=%g\n",i, val);
return GRIB_ENCODING_ERROR;
}
}
return GRIB_SUCCESS;
}*/
static int check_geometry(grib_handle* handle, const grib_util_grid_spec2* spec,
size_t data_values_count, int specified_as_global)
{
int err = 0;
if (spec->pl && spec->pl_size != 0 &&
(spec->grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_GG || spec->grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_ROTATED_GG)) {
if (specified_as_global) {
char msg[100] = {0,};
size_t sum = 0;
strcpy(msg, "Specified to be global (in spec)");
sum = sum_of_pl_array(spec->pl, spec->pl_size);
if (sum != data_values_count) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Invalid reduced gaussian grid: %s but data_values_count != sum_of_pl_array (%ld!=%ld)\n",
msg, (long)data_values_count, (long)sum);
return GRIB_WRONG_GRID;
}
}
}
return err;
}
#if 0
/* Check what is coded in the handle is what is requested by the spec. */
/* Return GRIB_SUCCESS if the geometry matches, otherwise the error code */
static int check_handle_against_spec(grib_handle* handle, const long edition,
const grib_util_grid_spec2* spec, int global_grid)
{
int err = 0;
int check_latitudes = 1;
int check_longitudes = 1;
long angleSubdivisions = 0;
double angular_precision = 1.0/1000.0; /* millidegree by default */
double tolerance = 0;
if (edition == 2) {
return GRIB_SUCCESS; /* For now only do checks on edition 1 */
}
if ((err = grib_get_long(handle, "angleSubdivisions", &angleSubdivisions))==GRIB_SUCCESS) {
angular_precision = 1.0/angleSubdivisions;
}
tolerance = angular_precision/2.0;
if (spec->grid_type == GRIB_UTIL_GRID_SPEC_POLAR_STEREOGRAPHIC ||
spec->grid_type == GRIB_UTIL_GRID_SPEC_SH)
{
return GRIB_SUCCESS;
}
/* Cannot check latitudes of Gaussian grids because are always sub-millidegree */
/* and for GRIB1 will differ from the encoded values. We accept this discrepancy! */
if (spec->grid_type == GRIB_UTIL_GRID_SPEC_REGULAR_GG ||
spec->grid_type == GRIB_UTIL_GRID_SPEC_ROTATED_GG ||
spec->grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_GG ||
spec->grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_ROTATED_GG)
{
if (edition == 1) {
check_latitudes = 0;
}
}
/* GRIB-922 */
/* Specification was to make the resulting grid global so no point checking for */
/* input lat/lon values as they would be reset by setting the "global" key to 1 */
if (global_grid)
{
check_latitudes = check_longitudes = 0;
}
if (check_latitudes) {
double lat1, lat2;
const double lat1spec = normalise_angle(spec->latitudeOfFirstGridPointInDegrees);
const double lat2spec = normalise_angle(spec->latitudeOfLastGridPointInDegrees);
if ((err = grib_get_double(handle, "latitudeOfFirstGridPointInDegrees", &lat1))!=0) return err;
if ((err = grib_get_double(handle, "latitudeOfLastGridPointInDegrees", &lat2))!=0) return err;
lat1 = normalise_angle(lat1);
lat2 = normalise_angle(lat2);
if (angle_too_small(lat1spec, angular_precision)) {
fprintf(stderr, "Failed to encode latitude of first grid point %.10e: less than angular precision\n",lat1spec);
return GRIB_WRONG_GRID;
}
if (angle_too_small(lat2spec, angular_precision)) {
fprintf(stderr, "Failed to encode latitude of last grid point %.10e: less than angular precision\n", lat2spec);
return GRIB_WRONG_GRID;
}
if (!DBL_EQUAL(lat1spec, lat1, tolerance)) {
fprintf(stderr, "Failed to encode latitude of first grid point: spec=%.10e val=%.10e\n", lat1spec, lat1);
return GRIB_WRONG_GRID;
}
if (!DBL_EQUAL(lat2spec, lat2, tolerance)) {
fprintf(stderr, "Failed to encode latitude of last grid point: spec=%.10e val=%.10e\n", lat2spec, lat2);
return GRIB_WRONG_GRID;
}
}
if (check_longitudes) {
double lon1, lon2;
const double lon1spec = normalise_angle(spec->longitudeOfFirstGridPointInDegrees);
const double lon2spec = normalise_angle(spec->longitudeOfLastGridPointInDegrees);
if ((err = grib_get_double(handle, "longitudeOfFirstGridPointInDegrees", &lon1))!=0) return err;
if ((err = grib_get_double(handle, "longitudeOfLastGridPointInDegrees", &lon2))!=0) return err;
lon1 = normalise_angle(lon1);
lon2 = normalise_angle(lon2);
if (angle_too_small(lon1spec, angular_precision)) {
fprintf(stderr, "Failed to encode longitude of first grid point %.10e: less than angular precision\n", lon1spec);
return GRIB_WRONG_GRID;
}
if (angle_too_small(lon2spec, angular_precision)) {
fprintf(stderr, "Failed to encode longitude of last grid point %.10e: less than angular precision\n", lon2spec);
return GRIB_WRONG_GRID;
}
if (!DBL_EQUAL(lon1spec, lon1, tolerance)) {
fprintf(stderr, "Failed to encode longitude of first grid point: spec=%.10e val=%.10e\n", lon1spec, lon1);
return GRIB_WRONG_GRID;
}
if (!DBL_EQUAL(lon2spec, lon2, tolerance)){
fprintf(stderr, "Failed to encode longitude of last grid point: spec=%.10e val=%.10e\n", lon2spec, lon2);
return GRIB_WRONG_GRID;
}
}
if (spec->grid_type == GRIB_UTIL_GRID_SPEC_ROTATED_LL ||
spec->grid_type == GRIB_UTIL_GRID_SPEC_ROTATED_GG ||
spec->grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_ROTATED_GG)
{
double latp, lonp;
const double latspec = normalise_angle(spec->latitudeOfSouthernPoleInDegrees);
const double lonspec = normalise_angle(spec->longitudeOfSouthernPoleInDegrees);
if ((err = grib_get_double(handle, "latitudeOfSouthernPoleInDegrees", &latp))!=0) return err;
if ((err = grib_get_double(handle, "longitudeOfSouthernPoleInDegrees", &lonp))!=0) return err;
latp = normalise_angle(latp);
lonp = normalise_angle(lonp);
if (!DBL_EQUAL(latspec, latp, tolerance)) {
fprintf(stderr, "Failed to encode latitude of southern pole: spec=%.10e val=%.10e\n",latspec,latp);
return GRIB_WRONG_GRID;
}
if (!DBL_EQUAL(lonspec, lonp, tolerance)) {
fprintf(stderr, "Failed to encode longitude of southern pole: spec=%.10e val=%.10e\n",lonspec,lonp);
return GRIB_WRONG_GRID;
}
}
return GRIB_SUCCESS;
}
#endif
static const char* get_grid_type_name(const int spec_grid_type)
{
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_REGULAR_LL)
return "regular_ll";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_ROTATED_LL)
return "rotated_ll";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_REGULAR_GG)
return "regular_gg";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_ROTATED_GG)
return "rotated_gg";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_LL)
return "reduced_ll";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_POLAR_STEREOGRAPHIC)
return "polar_stereographic";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_GG)
return "reduced_gg";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_SH)
return "sh";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_ROTATED_GG)
return "reduced_rotated_gg";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_LAMBERT_AZIMUTHAL_EQUAL_AREA)
return "lambert_azimuthal_equal_area";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_LAMBERT_CONFORMAL)
return "lambert";
if (spec_grid_type == GRIB_UTIL_GRID_SPEC_UNSTRUCTURED)
return "unstructured_grid";
return NULL;
}
static int is_constant_field(const double missingValue, const double* data_values, size_t data_values_count)
{
int ii = 0;
int constant = 1;
double value = missingValue;
for (ii = 0; ii < data_values_count; ii++) {
if (data_values[ii] != missingValue) {
if (value == missingValue) {
value = data_values[ii];
}
else {
if (value != data_values[ii]) {
constant = 0;
break;
}
}
}
}
return constant;
}
grib_handle* grib_util_set_spec(grib_handle* h,
const grib_util_grid_spec* spec,
const grib_util_packing_spec* packing_spec,
int flags,
const double* data_values,
size_t data_values_count,
int* err)
{
/* Create a spec v2 from spec v1 */
grib_util_grid_spec2 spec2;
Assert(h);
spec2.grid_type = spec->grid_type;
spec2.Ni = spec->Ni;
spec2.Nj = spec->Nj;
spec2.iDirectionIncrementInDegrees = spec->iDirectionIncrementInDegrees;
spec2.jDirectionIncrementInDegrees = spec->jDirectionIncrementInDegrees;
spec2.longitudeOfFirstGridPointInDegrees = spec->longitudeOfFirstGridPointInDegrees;
spec2.longitudeOfLastGridPointInDegrees = spec->longitudeOfLastGridPointInDegrees;
spec2.latitudeOfFirstGridPointInDegrees = spec->latitudeOfFirstGridPointInDegrees;
spec2.latitudeOfLastGridPointInDegrees = spec->latitudeOfLastGridPointInDegrees;
spec2.uvRelativeToGrid = spec->uvRelativeToGrid;
spec2.latitudeOfSouthernPoleInDegrees = spec->latitudeOfSouthernPoleInDegrees;
spec2.longitudeOfSouthernPoleInDegrees = spec->longitudeOfSouthernPoleInDegrees;
spec2.iScansNegatively = spec->iScansNegatively;
spec2.jScansPositively = spec->jScansPositively;
spec2.N = spec->N;
spec2.bitmapPresent = spec->bitmapPresent;
spec2.missingValue = spec->missingValue;
spec2.pl = spec->pl;
spec2.pl_size = spec->pl_size;
spec2.truncation = spec->truncation;
spec2.orientationOfTheGridInDegrees = spec->orientationOfTheGridInDegrees;
spec2.DyInMetres = spec->DyInMetres;
spec2.DxInMetres = spec->DxInMetres;
/* New data members of spec2 take default values */
spec2.angleOfRotationInDegrees = 0;
spec2.grid_name = NULL;
return grib_util_set_spec2(
h,
&spec2,
packing_spec,
flags,
data_values,
data_values_count,
err);
}
grib_handle* grib_util_set_spec2(grib_handle* h,
const grib_util_grid_spec2* spec,
const grib_util_packing_spec* packing_spec,
int flags,
const double* data_values,
size_t data_values_count,
int* err)
{
#define SET_LONG_VALUE(n, v) \
do { \
Assert(count < 1024); \
values[count].name = n; \
values[count].type = GRIB_TYPE_LONG; \
values[count].long_value = v; \
count++; \
} while (0)
#define SET_DOUBLE_VALUE(n, v) \
do { \
Assert(count < 1024); \
values[count].name = n; \
values[count].type = GRIB_TYPE_DOUBLE; \
values[count].double_value = v; \
count++; \
} while (0)
#define SET_STRING_VALUE(n, v) \
do { \
Assert(count < 1024); \
values[count].name = n; \
values[count].type = GRIB_TYPE_STRING; \
values[count].string_value = v; \
count++; \
} while (0)
#define COPY_SPEC_LONG(x) \
do { \
Assert(count < 1024); \
values[count].name = #x; \
values[count].type = GRIB_TYPE_LONG; \
values[count].long_value = spec->x; \
count++; \
} while (0)
#define COPY_SPEC_DOUBLE(x) \
do { \
Assert(count < 1024); \
values[count].name = #x; \
values[count].type = GRIB_TYPE_DOUBLE; \
values[count].double_value = spec->x; \
count++; \
} while (0)
grib_values values[1024] = {{0,},};
size_t count = 0;
int i;
long editionNumber;
grib_handle* outh = NULL;
grib_handle* tmp = NULL;
const char* grid_type = NULL;
char name[1024];
char input_grid_type[100];
char input_packing_type[100];
long input_bits_per_value = 0;
long input_decimal_scale_factor = 0;
size_t len = 100;
size_t input_grid_type_len = 100;
double laplacianOperator;
int packingTypeIsSet = 0;
int setSecondOrder = 0;
int setJpegPacking = 0;
int setCcsdsPacking = 0;
int convertEditionEarlier = 0; /* For cases when we cannot set some keys without converting */
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);
if (!packing_spec) {
packing_spec = &default_packing_spec;
}
/* Get edition number from input handle */
if ((*err = grib_get_long(h, "edition", &editionNumber)) != 0) {
grib_context* c = grib_context_get_default();
if (c->write_on_fail)
grib_write_message(h, "error.grib", "w");
return NULL;
}
if (packing_spec->deleteLocalDefinition) {
SET_LONG_VALUE("deleteLocalDefinition", 1);
}
len = 100;
grib_get_string(h, "packingType", input_packing_type, &len);
grib_get_long(h, "bitsPerValue", &input_bits_per_value);
grib_get_long(h, "decimalScaleFactor", &input_decimal_scale_factor);
if (h->context->debug == -1) {
fprintf(stderr, "ECCODES DEBUG grib_util: input_packing_type = %s\n", input_packing_type);
fprintf(stderr, "ECCODES DEBUG grib_util: input_bits_per_value = %ld\n", input_bits_per_value);
fprintf(stderr, "ECCODES DEBUG grib_util: input_decimal_scale_factor = %ld\n", input_decimal_scale_factor);
}
/*if ( (*err=check_values(data_values, data_values_count))!=GRIB_SUCCESS ) {
fprintf(stderr,"GRIB_UTIL_SET_SPEC: Data values check failed! %s\n", grib_get_error_message(*err));
goto cleanup;
}*/
if (flags & GRIB_UTIL_SET_SPEC_FLAGS_ONLY_PACKING) {
if (packing_spec->packing == GRIB_UTIL_PACKING_USE_PROVIDED &&
strcmp(input_packing_type, "grid_simple_matrix")) {
switch (packing_spec->packing_type) {
case GRIB_UTIL_PACKING_TYPE_SPECTRAL_COMPLEX:
if (strcmp(input_packing_type, "spectral_complex") && !strcmp(input_packing_type, "spectral_simple"))
SET_STRING_VALUE("packingType", "spectral_complex");
break;
case GRIB_UTIL_PACKING_TYPE_SPECTRAL_SIMPLE:
if (strcmp(input_packing_type, "spectral_simple") && !strcmp(input_packing_type, "spectral_complex"))
SET_STRING_VALUE("packingType", "spectral_simple");
break;
case GRIB_UTIL_PACKING_TYPE_GRID_SIMPLE:
if (strcmp(input_packing_type, "grid_simple") && !strcmp(input_packing_type, "grid_complex"))
SET_STRING_VALUE("packingType", "grid_simple");
break;
case GRIB_UTIL_PACKING_TYPE_GRID_SIMPLE_MATRIX:
SET_STRING_VALUE("packingType", "grid_simple_matrix");
break;
case GRIB_UTIL_PACKING_TYPE_GRID_COMPLEX:
if (strcmp(input_packing_type, "grid_complex") && !strcmp(input_packing_type, "grid_simple"))
SET_STRING_VALUE("packingType", "grid_complex");
break;
case GRIB_UTIL_PACKING_TYPE_JPEG:
if (strcmp(input_packing_type, "grid_jpeg") && !strcmp(input_packing_type, "grid_simple"))
SET_STRING_VALUE("packingType", "grid_jpeg");
break;
case GRIB_UTIL_PACKING_TYPE_CCSDS:
if (strcmp(input_packing_type, "grid_ccsds") && !strcmp(input_packing_type, "grid_simple"))
SET_STRING_VALUE("packingType", "grid_ccsds");
break;
case GRIB_UTIL_PACKING_TYPE_GRID_SECOND_ORDER:
/* we delay the set of grid_second_order because we don't want
to do it on a field with bitsPerValue=0 */
setSecondOrder = 1;
break;
default:
fprintf(stderr, "invalid packing_spec->packing_type = %ld\n", (long)packing_spec->packing_type);
*err = GRIB_INTERNAL_ERROR;
goto cleanup;
break;
}
}
switch (packing_spec->accuracy) {
case GRIB_UTIL_ACCURACY_SAME_BITS_PER_VALUES_AS_INPUT:
break;
case GRIB_UTIL_ACCURACY_USE_PROVIDED_BITS_PER_VALUES:
if (input_bits_per_value != packing_spec->bitsPerValue)
SET_LONG_VALUE("bitsPerValue", packing_spec->bitsPerValue);
break;
case GRIB_UTIL_ACCURACY_SAME_DECIMAL_SCALE_FACTOR_AS_INPUT:
break;
case GRIB_UTIL_ACCURACY_USE_PROVIDED_DECIMAL_SCALE_FACTOR:
if (input_decimal_scale_factor != packing_spec->decimalScaleFactor)
SET_LONG_VALUE("decimalScaleFactor", packing_spec->decimalScaleFactor);
break;
default:
fprintf(stderr, "invalid packing_spec->accuracy = %ld\n", (long)packing_spec->accuracy);
*err = GRIB_INTERNAL_ERROR;
goto cleanup;
break;
}
/*nothing to be changed*/
if (count == 0) {
*err = -1;
return h;
}
if (h->context->debug == -1) {
print_values(h->context, spec, data_values, data_values_count, values, count);
}
if ((*err = grib_set_values(h, values, count)) != 0) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Cannot set values %s\n", grib_get_error_message(*err));
for (i = 0; i < count; i++)
if (values[i].error)
fprintf(stderr, " %s %s\n", values[i].name, grib_get_error_message(values[i].error));
goto cleanup;
}
if (h->context->debug == -1) {
int j = 0;
fprintf(stderr, "ECCODES DEBUG grib_util: grib_set_double_array\n");
for (j = 0; j < 20; j++)
fprintf(stderr, "ECCODES DEBUG grib_util %g\n", data_values[j]);
fprintf(stderr, "ECCODES DEBUG grib_util: data_values_count=%d \n", (int)data_values_count);
}
if ((*err = grib_set_double_array(h, "values", data_values, data_values_count)) != 0) {
goto cleanup;
}
if (h->context->debug == -1)
fprintf(stderr, "ECCODES DEBUG grib_util: done grib_set_double_array \n");
/* convert to second_order if not constant field */
if (setSecondOrder) {
int constant = 0;
double missingValue = 0;
grib_get_double(h, "missingValue", &missingValue);
constant = is_constant_field(missingValue, data_values, data_values_count);
if (!constant) {
size_t packTypeLen;
if (editionNumber == 1) {
long numberOfGroups;
grib_handle* htmp = grib_handle_clone(h);
packTypeLen = 17;
grib_set_string(htmp, "packingType", "grid_second_order", &packTypeLen);
grib_get_long(htmp, "numberOfGroups", &numberOfGroups);
/* GRIBEX is not able to decode overflown numberOfGroups with SPD */
if (numberOfGroups > 65534 && h->context->no_spd) {
packTypeLen = 24;
grib_set_string(h, "packingType", "grid_second_order_no_SPD", &packTypeLen);
grib_handle_delete(htmp);
}
else {
grib_handle_delete(h);
h = htmp;
}
}
else {
packTypeLen = 17;
grib_set_string(h, "packingType", "grid_second_order", &packTypeLen);
grib_set_double_array(h, "values", data_values, data_values_count);
}
}
else {
if (h->context->gribex_mode_on) {
h->context->gribex_mode_on = 0;
grib_set_double_array(h, "values", data_values, data_values_count);
h->context->gribex_mode_on = 1;
}
}
}
return h;
}
grid_type = get_grid_type_name(spec->grid_type);
if (grid_type == NULL) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Unknown grid type: %d\n", spec->grid_type);
*err = GRIB_NOT_IMPLEMENTED;
return NULL;
}
SET_STRING_VALUE("gridType", grid_type);
/* The "pl" is given from the template, but "section_copy" will take care of setting the right headers */
{
switch (spec->grid_type) {
case GRIB_UTIL_GRID_SPEC_REDUCED_GG:
case GRIB_UTIL_GRID_SPEC_REDUCED_ROTATED_GG:
/* Choose a sample with the right Gaussian number and edition */
sprintf(name, "%s_pl_%ld_grib%ld", grid_type, spec->N, editionNumber);
if (spec->pl && spec->pl_size) {
/* GRIB-834: pl is given so can use any of the reduced_gg_pl samples */
sprintf(name, "%s_pl_grib%ld", grid_type, editionNumber);
}
break;
case GRIB_UTIL_GRID_SPEC_LAMBERT_AZIMUTHAL_EQUAL_AREA:
case GRIB_UTIL_GRID_SPEC_UNSTRUCTURED:
if (editionNumber == 1) { /* This grid type is not available in edition 1 */
if (h->context->debug == -1)
fprintf(stderr,
"ECCODES DEBUG grib_util: '%s' specified "
"but input is GRIB1. Output must be a higher edition!\n",
grid_type);
convertEditionEarlier = 1;
}
sprintf(name, "GRIB%ld", editionNumber);
break;
case GRIB_UTIL_GRID_SPEC_LAMBERT_CONFORMAL:
sprintf(name, "GRIB%ld", editionNumber);
break;
default:
sprintf(name, "%s_pl_grib%ld", grid_type, editionNumber);
}
if (spec->pl && spec->grid_name) {
/* Cannot have BOTH pl and grid name specified */
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Cannot set BOTH pl and grid_name.\n");
goto cleanup;
}
if (spec->grid_name) {
sprintf(name, "%s_grib%ld", spec->grid_name, editionNumber);
}
}
/* TODO: recycle tmp handle */
tmp = grib_handle_new_from_samples(NULL, name);
if (!tmp) {
*err = GRIB_INVALID_FILE;
return NULL;
}
/* Set grid */
switch (spec->grid_type) {
case GRIB_UTIL_GRID_SPEC_REGULAR_LL:
case GRIB_UTIL_GRID_SPEC_ROTATED_LL:
COPY_SPEC_LONG(bitmapPresent);
if (spec->missingValue)
COPY_SPEC_DOUBLE(missingValue);
SET_LONG_VALUE("ijDirectionIncrementGiven", 1);
if (editionNumber == 1) {
/* GRIB-863: GRIB1 cannot represent increments less than a millidegree */
if (!angle_can_be_encoded(h, spec->iDirectionIncrementInDegrees) ||
!angle_can_be_encoded(h, spec->jDirectionIncrementInDegrees)) {
grib1_high_resolution_fix = 1;
/* Set flag to compute the increments */
SET_LONG_VALUE("ijDirectionIncrementGiven", 0);
}
}
/* default iScansNegatively=0 jScansPositively=0 is ok */
COPY_SPEC_LONG(iScansNegatively);
COPY_SPEC_LONG(jScansPositively);
COPY_SPEC_LONG(Ni);
COPY_SPEC_LONG(Nj);
COPY_SPEC_DOUBLE(iDirectionIncrementInDegrees);
COPY_SPEC_DOUBLE(jDirectionIncrementInDegrees);
COPY_SPEC_DOUBLE(longitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(longitudeOfLastGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfLastGridPointInDegrees);
break;
case GRIB_UTIL_GRID_SPEC_REGULAR_GG:
case GRIB_UTIL_GRID_SPEC_ROTATED_GG:
COPY_SPEC_LONG(bitmapPresent);
if (spec->missingValue)
COPY_SPEC_DOUBLE(missingValue);
SET_LONG_VALUE("ijDirectionIncrementGiven", 1);
/* TODO: add Assert */
COPY_SPEC_LONG(Ni);
COPY_SPEC_DOUBLE(iDirectionIncrementInDegrees);
COPY_SPEC_LONG(Nj);
COPY_SPEC_LONG(N);
/* TODO: Compute here ... */
COPY_SPEC_DOUBLE(longitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(longitudeOfLastGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfLastGridPointInDegrees);
break;
case GRIB_UTIL_GRID_SPEC_REDUCED_LL:
COPY_SPEC_LONG(bitmapPresent);
if (spec->missingValue)
COPY_SPEC_DOUBLE(missingValue);
SET_LONG_VALUE("ijDirectionIncrementGiven", 0);
COPY_SPEC_LONG(Nj);
COPY_SPEC_DOUBLE(longitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(longitudeOfLastGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfLastGridPointInDegrees);
break;
case GRIB_UTIL_GRID_SPEC_POLAR_STEREOGRAPHIC:
COPY_SPEC_LONG(bitmapPresent);
if (spec->missingValue)
COPY_SPEC_DOUBLE(missingValue);
COPY_SPEC_DOUBLE(longitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfFirstGridPointInDegrees);
COPY_SPEC_LONG(Ni);
COPY_SPEC_LONG(Nj);
/* default iScansNegatively=0 jScansPositively=0 is ok */
COPY_SPEC_LONG(iScansNegatively);
COPY_SPEC_LONG(jScansPositively);
COPY_SPEC_DOUBLE(orientationOfTheGridInDegrees);
COPY_SPEC_LONG(DxInMetres);
COPY_SPEC_LONG(DyInMetres);
break;
case GRIB_UTIL_GRID_SPEC_LAMBERT_AZIMUTHAL_EQUAL_AREA:
COPY_SPEC_LONG(bitmapPresent);
if (spec->missingValue)
COPY_SPEC_DOUBLE(missingValue);
COPY_SPEC_DOUBLE(longitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfFirstGridPointInDegrees);
COPY_SPEC_LONG(Ni); /* same as Nx */
COPY_SPEC_LONG(Nj); /* same as Ny */
/* TODO
* pass in extra keys e.g. Dx, Dy, standardParallel and centralLongitude
*/
/*
COPY_SPEC_LONG(DxInMetres);
COPY_SPEC_LONG(DyInMetres);
COPY_SPEC_LONG(xDirectionGridLengthInMillimetres);
COPY_SPEC_LONG(yDirectionGridLengthInMillimetres);
COPY_SPEC_LONG(standardParallelInMicrodegrees);
COPY_SPEC_LONG(centralLongitudeInMicrodegrees);
*/
break;
case GRIB_UTIL_GRID_SPEC_UNSTRUCTURED:
COPY_SPEC_LONG(bitmapPresent);
if (spec->missingValue)
COPY_SPEC_DOUBLE(missingValue);
/*
* TODO: Other keys
*/
break;
case GRIB_UTIL_GRID_SPEC_LAMBERT_CONFORMAL:
COPY_SPEC_LONG(bitmapPresent);
if (spec->missingValue)
COPY_SPEC_DOUBLE(missingValue);
COPY_SPEC_DOUBLE(longitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfFirstGridPointInDegrees);
COPY_SPEC_LONG(Ni); /* same as Nx */
COPY_SPEC_LONG(Nj); /* same as Ny */
/*
* Note: DxInMetres and DyInMetres
* should be 'double' and not integer. WMO GRIB2 uses millimetres!
* TODO:
* Add other keys like Latin1, LoV etc
*err = GRIB_NOT_IMPLEMENTED;
goto cleanup;
*/
break;
case GRIB_UTIL_GRID_SPEC_REDUCED_GG:
case GRIB_UTIL_GRID_SPEC_REDUCED_ROTATED_GG:
COPY_SPEC_LONG(bitmapPresent);
if (spec->missingValue)
COPY_SPEC_DOUBLE(missingValue);
SET_LONG_VALUE("ijDirectionIncrementGiven", 0);
COPY_SPEC_LONG(Nj);
COPY_SPEC_LONG(N);
COPY_SPEC_DOUBLE(longitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(longitudeOfLastGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfFirstGridPointInDegrees);
COPY_SPEC_DOUBLE(latitudeOfLastGridPointInDegrees);
break;
case GRIB_UTIL_GRID_SPEC_SH:
*err = grib_get_string(h, "gridType", input_grid_type, &input_grid_type_len);
SET_LONG_VALUE("J", spec->truncation);
SET_LONG_VALUE("K", spec->truncation);
SET_LONG_VALUE("M", spec->truncation);
if (packing_spec->packing_type == GRIB_UTIL_PACKING_TYPE_SPECTRAL_COMPLEX) {
const long JS = spec->truncation < 20 ? spec->truncation : 20;
SET_STRING_VALUE("packingType", "spectral_complex");
packingTypeIsSet = 1;
SET_LONG_VALUE("JS", JS);
SET_LONG_VALUE("KS", JS);
SET_LONG_VALUE("MS", JS);
if (packing_spec->packing == GRIB_UTIL_PACKING_USE_PROVIDED && editionNumber == 2) {
SET_LONG_VALUE("computeLaplacianOperator", 1);
}
else if ((!(*err) && strcmp(input_grid_type, "sh")) || packing_spec->computeLaplacianOperator) {
SET_LONG_VALUE("computeLaplacianOperator", 1);
if (packing_spec->truncateLaplacian)
SET_LONG_VALUE("truncateLaplacian", 1);
}
else {
SET_LONG_VALUE("computeLaplacianOperator", 0);
*err = grib_get_double(h, "laplacianOperator", &laplacianOperator);
if (packing_spec->truncateLaplacian)
SET_LONG_VALUE("truncateLaplacian", 1);
SET_DOUBLE_VALUE("laplacianOperator", packing_spec->laplacianOperator);
if (laplacianOperator) {
SET_DOUBLE_VALUE("laplacianOperator", laplacianOperator);
}
}
}
break;
}
/* Set rotation */
switch (spec->grid_type) {
case GRIB_UTIL_GRID_SPEC_ROTATED_LL:
case GRIB_UTIL_GRID_SPEC_ROTATED_GG:
case GRIB_UTIL_GRID_SPEC_REDUCED_ROTATED_GG:
COPY_SPEC_LONG(uvRelativeToGrid);
COPY_SPEC_DOUBLE(latitudeOfSouthernPoleInDegrees);
COPY_SPEC_DOUBLE(longitudeOfSouthernPoleInDegrees);
COPY_SPEC_DOUBLE(angleOfRotationInDegrees);
break;
}
/* process packing options */
if (!packingTypeIsSet &&
packing_spec->packing == GRIB_UTIL_PACKING_USE_PROVIDED &&
strcmp(input_packing_type, "grid_simple_matrix")) {
switch (packing_spec->packing_type) {
case GRIB_UTIL_PACKING_TYPE_SPECTRAL_COMPLEX:
if (strcmp(input_packing_type, "spectral_complex") && !strcmp(input_packing_type, "spectral_simple"))
SET_STRING_VALUE("packingType", "spectral_complex");
break;
case GRIB_UTIL_PACKING_TYPE_SPECTRAL_SIMPLE:
if (strcmp(input_packing_type, "spectral_simple") && !strcmp(input_packing_type, "spectral_complex"))
SET_STRING_VALUE("packingType", "spectral_simple");
break;
case GRIB_UTIL_PACKING_TYPE_GRID_SIMPLE:
if (strcmp(input_packing_type, "grid_simple") && !strcmp(input_packing_type, "grid_complex"))
SET_STRING_VALUE("packingType", "grid_simple");
break;
case GRIB_UTIL_PACKING_TYPE_GRID_COMPLEX:
if (strcmp(input_packing_type, "grid_complex") && !strcmp(input_packing_type, "grid_simple"))
SET_STRING_VALUE("packingType", "grid_complex");
break;
case GRIB_UTIL_PACKING_TYPE_JPEG:
/* Have to delay JPEG packing:
* Reason 1: It is not available in GRIB1 and so we have to wait until we change edition
* Reason 2: It has to be done AFTER we set the data values
*/
if (strcmp(input_packing_type, "grid_jpeg") && !strcmp(input_packing_type, "grid_simple"))
setJpegPacking = 1;
break;
case GRIB_UTIL_PACKING_TYPE_CCSDS:
/* Have to delay CCSDS packing:
* Reason 1: It is not available in GRIB1 and so we have to wait until we change edition
* Reason 2: It has to be done AFTER we set the data values
*/
if (strcmp(input_packing_type, "grid_ccsds") && !strcmp(input_packing_type, "grid_simple"))
setCcsdsPacking = 1;
break;
case GRIB_UTIL_PACKING_TYPE_GRID_SECOND_ORDER:
/* we delay the set of grid_second_order because we don't want
to do it on a field with bitsPerValue=0 */
setSecondOrder = 1;
break;
default:
fprintf(stderr, "invalid packing_spec->packing_type = %ld\n", (long)packing_spec->packing_type);
*err = GRIB_INTERNAL_ERROR;
goto cleanup;
break;
}
}
if (strcmp(input_packing_type, "grid_simple_matrix") == 0) {
long numberOfDirections, numberOfFrequencies;
int keep_matrix = h->context->keep_matrix;
if (packing_spec->packing_type == GRIB_UTIL_PACKING_TYPE_GRID_SIMPLE) {
keep_matrix = 0; /* ECC-911 */
}
if (keep_matrix) {
int ret;
SET_STRING_VALUE("packingType", "grid_simple_matrix");
ret = grib_get_long(h, "numberOfDirections", &numberOfDirections);
if (ret == GRIB_SUCCESS) {
grib_get_long(h, "numberOfDirections", &numberOfDirections);
SET_LONG_VALUE("NC1", numberOfDirections);
grib_get_long(h, "numberOfFrequencies", &numberOfFrequencies);
SET_LONG_VALUE("NC2", numberOfFrequencies);
SET_LONG_VALUE("physicalFlag1", 1);
SET_LONG_VALUE("physicalFlag2", 2);
SET_LONG_VALUE("NR", 1);
SET_LONG_VALUE("NC", 1);
}
}
else {
SET_STRING_VALUE("packingType", "grid_simple");
}
}
switch (packing_spec->accuracy) {
case GRIB_UTIL_ACCURACY_SAME_BITS_PER_VALUES_AS_INPUT: {
long bitsPerValue = 0;
Assert(grib_get_long(h, "bitsPerValue", &bitsPerValue) == 0);
SET_LONG_VALUE("bitsPerValue", bitsPerValue);
} break;
case GRIB_UTIL_ACCURACY_USE_PROVIDED_BITS_PER_VALUES:
SET_LONG_VALUE("bitsPerValue", packing_spec->bitsPerValue);
break;
case GRIB_UTIL_ACCURACY_SAME_DECIMAL_SCALE_FACTOR_AS_INPUT: {
long decimalScaleFactor = 0;
Assert(grib_get_long(h, "decimalScaleFactor", &decimalScaleFactor) == 0);
SET_LONG_VALUE("decimalScaleFactor", decimalScaleFactor);
} break;
case GRIB_UTIL_ACCURACY_USE_PROVIDED_DECIMAL_SCALE_FACTOR:
SET_LONG_VALUE("decimalScaleFactor", packing_spec->decimalScaleFactor);
break;
default:
fprintf(stderr, "invalid packing_spec->accuracy = %ld\n", (long)packing_spec->accuracy);
grib_handle_delete(tmp);
*err = GRIB_INTERNAL_ERROR;
goto cleanup;
break;
}
if (packing_spec->extra_settings_count) {
for (i = 0; i < packing_spec->extra_settings_count; i++) {
Assert(count < 1024);
if (strcmp(packing_spec->extra_settings[i].name, "expandBoundingBox") == 0) {
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;
}
}
}
}
/* grib_write_message(h,"input.grib","w"); */
/* grib_write_message(tmp,"geo.grib","w"); */
/* copy product and local sections from h to tmp handle and store in outh */
if ((outh = grib_util_sections_copy(h, tmp, GRIB_SECTION_PRODUCT | GRIB_SECTION_LOCAL, err)) == NULL) {
goto cleanup;
}
grib_handle_delete(tmp);
Assert(*err == 0);
/* Set "pl" array if provided (For reduced Gaussian grids) */
/* See GRIB-857 */
Assert(spec->pl_size >= 0);
if (spec->pl && spec->pl_size == 0) {
fprintf(stderr, "pl array not NULL but pl_size == 0!\n");
goto cleanup;
}
if (spec->pl_size > 0 && spec->pl == NULL) {
fprintf(stderr, "pl_size not zero but pl array == NULL!\n");
goto cleanup;
}
if (spec->pl_size != 0 && (spec->grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_GG || spec->grid_type == GRIB_UTIL_GRID_SPEC_REDUCED_ROTATED_GG)) {
*err = grib_set_long_array(outh, "pl", spec->pl, spec->pl_size);
if (*err) {
fprintf(stderr, "SET_GRID_DATA_DESCRIPTION: Cannot set pl %s\n", grib_get_error_message(*err));
goto cleanup;
}
if (global_grid) {
size_t sum = sum_of_pl_array(spec->pl, spec->pl_size);
if (data_values_count != sum) {
fprintf(stderr, "invalid reduced gaussian grid: specified as global, data_values_count=%ld but sum of pl array=%ld\n",
(long)data_values_count, (long)sum);
*err = GRIB_WRONG_GRID;
goto cleanup;
}
}
}
if (h->context->debug == -1) {
fprintf(stderr, "ECCODES DEBUG grib_util: global_grid = %d\n", global_grid);
fprintf(stderr, "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));
if (h->context->write_on_fail)
grib_write_message(outh, "error.grib", "w");
goto cleanup;
}
}
if (convertEditionEarlier && packing_spec->editionNumber > 1) {
*err = grib_set_long(outh, "edition", packing_spec->editionNumber);
if (*err) {
fprintf(stderr, "SET_GRID_DATA_DESCRIPTION: Cannot convert to edition %ld.\n", packing_spec->editionNumber);
goto cleanup;
}
}
if ((*err = grib_set_values(outh, values, count)) != 0) {
fprintf(stderr, "SET_GRID_DATA_DESCRIPTION: Cannot set key values: %s\n", grib_get_error_message(*err));
for (i = 0; i < count; i++)
if (values[i].error)
fprintf(stderr, " %s %s\n", values[i].name, grib_get_error_message(values[i].error));
goto cleanup;
}
if ((*err = grib_set_double_array(outh, "values", data_values, data_values_count)) != 0) {
FILE* ferror;
size_t ii, lcount;
grib_context* c = grib_context_get_default();
ferror = fopen("error.data", "w");
lcount = 0;
fprintf(ferror, "# data_values_count=%ld\n", (long)data_values_count);
fprintf(ferror, "set values={ ");
for (ii = 0; ii < data_values_count - 1; ii++) {
fprintf(ferror, "%g, ", data_values[ii]);
if (lcount > 10) {
fprintf(ferror, "\n");
lcount = 0;
}
lcount++;
}
fprintf(ferror, "%g }", data_values[data_values_count - 1]);
fclose(ferror);
if (c->write_on_fail)
grib_write_message(outh, "error.grib", "w");
goto cleanup;
}
/* grib_write_message(outh,"h.grib","w"); */
/* if the field is empty GRIBEX is packing as simple*/
/* if (!strcmp(input_packing_type,"grid_simple_matrix")) {
long numberOfValues;
grib_get_long(outh,"numberOfValues",&numberOfValues);
if (numberOfValues==0) {
slen=11;
grib_set_string(outh,"packingType","grid_simple",&slen);
}
} */
if (grib1_high_resolution_fix) {
/* GRIB-863: must set increments to MISSING */
/* increments are not coded in message but computed */
if ((*err = grib_set_missing(outh, "iDirectionIncrement")) != 0) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Cannot set Di to missing: %s\n", grib_get_error_message(*err));
goto cleanup;
}
if ((*err = grib_set_missing(outh, "jDirectionIncrement")) != 0) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Cannot set Dj to missing: %s\n", grib_get_error_message(*err));
goto cleanup;
}
}
/*grib_dump_content(outh, stdout,"debug", ~0, NULL);*/
/* convert to second_order if not constant field. (Also see ECC-326) */
if (setSecondOrder) {
int constant = 0;
double missingValue = 0;
grib_get_double(outh, "missingValue", &missingValue);
constant = is_constant_field(missingValue, data_values, data_values_count);
if (!constant) {
if (editionNumber == 1) {
long numberOfGroups = 0;
grib_handle* htmp = grib_handle_clone(outh);
slen = 17;
grib_set_string(htmp, "packingType", "grid_second_order", &slen);
grib_get_long(htmp, "numberOfGroups", &numberOfGroups);
/* GRIBEX is not able to decode overflown numberOfGroups with SPD */
if (numberOfGroups > 65534 && outh->context->no_spd) {
slen = 24;
grib_set_string(outh, "packingType", "grid_second_order_no_SPD", &slen);
grib_handle_delete(htmp);
}
else {
grib_handle_delete(outh);
outh = htmp;
}
}
else {
slen = 17;
grib_set_string(outh, "packingType", "grid_second_order", &slen);
*err = grib_set_double_array(outh, "values", data_values, data_values_count);
if (*err != GRIB_SUCCESS) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: setting data values failed! %s\n", grib_get_error_message(*err));
goto cleanup;
}
}
}
else {
if (outh->context->gribex_mode_on) {
outh->context->gribex_mode_on = 0;
grib_set_double_array(outh, "values", data_values, data_values_count);
outh->context->gribex_mode_on = 1;
}
}
}
if (packing_spec->editionNumber && packing_spec->editionNumber != editionNumber) {
*err = grib_set_long(outh, "edition", packing_spec->editionNumber);
if (*err != GRIB_SUCCESS) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Failed to change edition to %ld: %s\n",
packing_spec->editionNumber, grib_get_error_message(*err));
if (h->context->write_on_fail)
grib_write_message(outh, "error.grib", "w");
goto cleanup;
}
}
if (editionNumber > 1 || packing_spec->editionNumber > 1) {
/* ECC-353 */
/* JPEG packing is not available in GRIB edition 1 and has to be done AFTER we set data values */
if (setJpegPacking == 1) {
*err = grib_set_string(outh, "packingType", "grid_jpeg", &slen);
if (*err != GRIB_SUCCESS) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Failed to change packingType to JPEG: %s\n",
grib_get_error_message(*err));
goto cleanup;
}
}
if (setCcsdsPacking == 1) {
*err = grib_set_string(outh, "packingType", "grid_ccsds", &slen);
if (*err != GRIB_SUCCESS) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Failed to change packingType to CCSDS: %s\n",
grib_get_error_message(*err));
goto cleanup;
}
}
}
if (packing_spec->deleteLocalDefinition) {
grib_set_long(outh, "deleteLocalDefinition", 1);
}
/* ECC-445 */
if (expandBoundingBox) {
Assert(!global_grid); /* ECC-576: "global" should not be set */
}
if ((*err = check_geometry(outh, spec, data_values_count, global_grid)) != GRIB_SUCCESS) {
fprintf(stderr, "GRIB_UTIL_SET_SPEC: Geometry check failed! %s\n", grib_get_error_message(*err));
if (h->context->write_on_fail)
grib_write_message(outh, "error.grib", "w");
goto cleanup;
}
/* Disable check: need to re-examine GRIB-864 */
#if 0
if ( (*err = check_handle_against_spec(outh, editionNumber, spec, global_grid)) != GRIB_SUCCESS)
{
grib_context* c=grib_context_get_default();
fprintf(stderr,"GRIB_UTIL_SET_SPEC: Geometry check failed! %s\n", grib_get_error_message(*err));
if (editionNumber == 1) {
fprintf(stderr,"Note: in GRIB edition 1 latitude and longitude values cannot be represented with sub-millidegree precision.\n");
}
if (c->write_on_fail)
grib_write_message(outh,"error.grib","w");
goto cleanup;
}
#endif
if (h->context->debug == -1)
fprintf(stderr, "ECCODES DEBUG: grib_util_set_spec end\n");
return outh;
cleanup:
if (outh)
grib_handle_delete(outh);
return NULL;
}
int grib_moments(grib_handle* h, double east, double north, double west, double south, int order, double* moments, long* count)
{
grib_iterator* iter = NULL;
int ret = 0, i, j, l;
size_t n = 0, numberOfPoints = 0;
double *lat, *lon, *values;
double vlat, vlon, val;
double dx, dy, ddx, ddy;
double mass, centroidX, centroidY;
double missingValue;
grib_context* c = grib_context_get_default();
ret = grib_get_size(h, "values", &n);
if (ret)
return ret;
lat = (double*)grib_context_malloc_clear(c, sizeof(double) * n);
lon = (double*)grib_context_malloc_clear(c, sizeof(double) * n);
values = (double*)grib_context_malloc_clear(c, sizeof(double) * n);
iter = grib_iterator_new(h, 0, &ret);
numberOfPoints = 0;
while (grib_iterator_next(iter, &vlat, &vlon, &val)) {
if (vlon >= east && vlon <= west && vlat >= south && vlat <= north) {
lat[numberOfPoints] = vlat;
lon[numberOfPoints] = vlon;
values[numberOfPoints] = val;
numberOfPoints++;
}
}
grib_iterator_delete(iter);
ret = grib_get_double(h, "missingValue", &missingValue);
centroidX = 0;
centroidY = 0;
mass = 0;
*count = 0;
for (i = 0; i < numberOfPoints; i++) {
if (values[i] != missingValue) {
centroidX += lon[i] * values[i];
centroidY += lat[i] * values[i];
mass += values[i];
(*count)++;
}
}
centroidX /= mass;
centroidY /= mass;
mass /= *count;
for (j = 0; j < order * order; j++)
moments[j] = 0;
for (i = 0; i < numberOfPoints; i++) {
if (values[i] != missingValue) {
dx = (lon[i] - centroidX);
dy = (lat[i] - centroidY);
ddx = 1;
for (j = 0; j < order; j++) {
ddy = 1;
for (l = 0; l < order; l++) {
moments[j * order + l] += ddx * ddy * values[i];
ddy *= dy;
}
ddx *= dx;
}
}
}
for (j = 0; j < order; j++) {
for (l = 0; l < order; l++) {
if (j + l > 1) {
moments[j * order + l] = pow(fabs(moments[j * order + l]), 1.0 / (j + l)) / *count;
}
else {
moments[j * order + l] /= *count;
}
}
}
grib_context_free(c, lat);
grib_context_free(c, lon);
grib_context_free(c, values);
(void)mass;
return ret;
}
/* Helper function for 'parse_keyval_string' */
static void set_value(grib_values* value, char* str, int equal)
{
char *p = 0, *q = 0, *s = 0;
char buf[1000] = {0,};
grib_context* c = grib_context_get_default();
value->equal = equal;
q = str;
while (*q != '/' && *q != 0)
q++;
if (*q == '/') {
s = grib_context_strdup(c, q + 1);
value->next = (grib_values*)grib_context_malloc_clear(c, sizeof(grib_values));
value->next->type = value->type;
value->next->name = grib_context_strdup(c, value->name);
set_value(value->next, s, equal);
grib_context_free(c, s);
}
memcpy(buf, str, q - str);
switch (value->type) {
case GRIB_TYPE_DOUBLE:
value->double_value = strtod(buf, &p);
if (*p != 0)
value->has_value = 1;
else if (!strcmp(str, "missing") ||
!strcmp(str, "MISSING") ||
!strcmp(str, "Missing")) {
value->type = GRIB_TYPE_MISSING;
value->has_value = 1;
}
break;
case GRIB_TYPE_LONG:
value->long_value = strtol(buf, &p, 10);
if (*p != 0)
value->has_value = 1;
else if (!strcmp(buf, "missing") ||
!strcmp(buf, "MISSING") ||
!strcmp(buf, "Missing")) {
value->type = GRIB_TYPE_MISSING;
value->has_value = 1;
}
break;
case GRIB_TYPE_STRING:
if (!strcmp(buf, "missing") ||
!strcmp(buf, "MISSING") ||
!strcmp(buf, "Missing")) {
value->type = GRIB_TYPE_MISSING;
value->has_value = 1;
}
else {
value->string_value = grib_context_strdup(c, buf);
value->has_value = 1;
}
break;
case GRIB_TYPE_UNDEFINED:
value->long_value = strtol(buf, &p, 10);
if (*p == 0) {
value->type = GRIB_TYPE_LONG;
value->has_value = 1;
}
else {
value->double_value = strtod(buf, &p);
if (*p == 0) {
value->type = GRIB_TYPE_DOUBLE;
value->has_value = 1;
}
else if (!strcmp(buf, "missing") ||
!strcmp(buf, "MISSING") ||
!strcmp(buf, "Missing")) {
value->type = GRIB_TYPE_MISSING;
value->has_value = 1;
}
else {
value->string_value = grib_context_strdup(c, buf);
value->type = GRIB_TYPE_STRING;
value->has_value = 1;
}
}
break;
}
}
/*
'grib_tool' Optional tool name which is printed on error. Can be NULL
'arg' The string to be parsed e.g. key1=value1,key2!=value2 etc
'values_required' If true then each key must have a value after it
'default_type' The default type e.g. GRIB_TYPE_UNDEFINED or GRIB_TYPE_DOUBLE
'values' The array we populate and return
'count' The number of elements
*/
int parse_keyval_string(const char* grib_tool, char* arg, int values_required, int default_type, grib_values values[], int* count)
{
char* p;
int i = 0;
if (arg == NULL) {
*count = 0;
return GRIB_SUCCESS;
}
p = strtok(arg, ",");
while (p != NULL) {
values[i].name = (char*)calloc(1, strlen(p) + 1);
Assert(values[i].name);
strcpy((char*)values[i].name, p);
p = strtok(NULL, ",");
i++;
if (i > *count)
return GRIB_ARRAY_TOO_SMALL;
}
*count = i;
for (i = 0; i < *count; i++) {
int equal = 1;
char* value = NULL;
if (values_required) {
/* Can be either k=v or k!=v */
p = (char*)values[i].name;
while (*p != '=' && *p != '!' && *p != '\0')
p++;
if (*p == '=') {
*p = '\0';
p++;
value = p;
equal = 1;
}
else if (*p == '!' && *(++p) == '=') {
*p = '\0';
*(p - 1) = '\0';
p++;
value = p;
equal = 0;
}
else {
return GRIB_INVALID_ARGUMENT;
}
}
p = (char*)values[i].name;
while (*p != ':' && *p != '\0')
p++;
if (*p == ':') {
values[i].type = grib_type_to_int(*(p + 1));
if (*(p + 1) == 'n')
values[i].type = GRIB_NAMESPACE;
*p = '\0';
p++;
}
else {
values[i].type = default_type;
}
if (values_required) {
if (strlen(value) == 0) {
if (grib_tool)
printf("%s error: no value provided for key \"%s\"\n", grib_tool, values[i].name);
else
printf("Error: no value provided for key \"%s\"\n", values[i].name);
return GRIB_INVALID_ARGUMENT;
}
set_value(&values[i], value, equal);
}
}
return GRIB_SUCCESS;
}
/* Return 1 if the productDefinitionTemplateNumber (GRIB2) is related to EPS */
int grib2_is_PDTN_EPS(long productDefinitionTemplateNumber)
{
return (
productDefinitionTemplateNumber == 1 || productDefinitionTemplateNumber == 11 ||
productDefinitionTemplateNumber == 33 || productDefinitionTemplateNumber == 34 || /*simulated (synthetic) satellite data*/
productDefinitionTemplateNumber == 41 || productDefinitionTemplateNumber == 43 || /*atmospheric chemical constituents*/
productDefinitionTemplateNumber == 45 || productDefinitionTemplateNumber == 47 /*aerosols*/
);
}
/* Return 1 if the productDefinitionTemplateNumber (GRIB2) is for atmospheric chemical constituents */
int grib2_is_PDTN_Chemical(long productDefinitionTemplateNumber)
{
return (
productDefinitionTemplateNumber == 40 ||
productDefinitionTemplateNumber == 41 ||
productDefinitionTemplateNumber == 42 ||
productDefinitionTemplateNumber == 43);
}
/* Return 1 if the productDefinitionTemplateNumber (GRIB2) is for
* atmospheric chemical constituents based on a distribution function */
int grib2_is_PDTN_ChemicalDistFunc(long productDefinitionTemplateNumber)
{
return (
productDefinitionTemplateNumber == 57 ||
productDefinitionTemplateNumber == 58 ||
productDefinitionTemplateNumber == 67 ||
productDefinitionTemplateNumber == 68);
}
/* Return 1 if the productDefinitionTemplateNumber (GRIB2) is for aerosols */
int grib2_is_PDTN_Aerosol(long productDefinitionTemplateNumber)
{
return (
productDefinitionTemplateNumber == 44 || /* Note: PDT 44 is deprecated. Use 48 instead */
productDefinitionTemplateNumber == 48 ||
productDefinitionTemplateNumber == 49 ||
productDefinitionTemplateNumber == 45 ||
productDefinitionTemplateNumber == 46 ||
productDefinitionTemplateNumber == 47);
}
/* Return 1 if the productDefinitionTemplateNumber (GRIB2) is for optical properties of aerosol */
int grib2_is_PDTN_AerosolOptical(long productDefinitionTemplateNumber)
{
/* Note: PDT 48 can be used for both plain aerosols as well as optical properties of aerosol.
* For the former user must set the optical wavelength range to missing.
*/
return (
productDefinitionTemplateNumber == 48 ||
productDefinitionTemplateNumber == 49);
}
/* Given some information about the type of grib2 parameter, return the productDefinitionTemplateNumber to use.
* All arguments are booleans (0 or 1)
* is_eps: ensemble or deterministic
* is_instant: instantaneous or interval-based
* etc
*/
int grib2_select_PDTN(int is_eps, int is_instant,
int is_chemical,
int is_chemical_distfn,
int is_aerosol,
int is_aerosol_optical)
{
/* At most one has to be set. All could be 0 */
/* Unfortunately if PDTN=48 then both aerosol and aerosol_optical can be 1! */
const int sum = is_chemical + is_chemical_distfn + is_aerosol + is_aerosol_optical;
Assert(sum == 0 || sum == 1 || sum == 2);
if (is_chemical) {
if (is_eps) {
if (is_instant)
return 41;
else
return 43;
}
else {
if (is_instant)
return 40;
else
return 42;
}
}
if (is_chemical_distfn) {
if (is_eps) {
if (is_instant)
return 58;
else
return 68;
}
else {
if (is_instant)
return 57;
else
return 67;
}
}
if (is_aerosol_optical) {
if (is_eps) {
if (is_instant)
return 49;
/* WMO does not have a non-instantaneous case here! */
}
else {
if (is_instant)
return 48;
/* WMO does not have a non-instantaneous case here! */
}
}
if (is_aerosol) {
if (is_eps) {
if (is_instant)
return 45;
else
return 47;
}
else {
if (is_instant)
return 48; /*44 is deprecated*/
else
return 46;
}
}
/* Fallthru case: default */
if (is_eps) {
if (is_instant)
return 1;
else
return 11;
}
else {
if (is_instant)
return 0;
else
return 8;
}
}
int is_grib_index_file(const char* filename)
{
FILE* fh;
char buf[8] = {0,};
const char* str = "GRBIDX";
int ret = 0;
size_t size = 0;
fh = fopen(filename, "r");
if (!fh)
return 0;
size = fread(buf, 1, 1, fh);
if (size != 1) {
fclose(fh);
return 0;
}
size = fread(buf, 6, 1, fh);
if (size != 1) {
fclose(fh);
return 0;
}
ret = !strcmp(buf, str);
fclose(fh);
return ret;
}
size_t sum_of_pl_array(const long* pl, size_t plsize)
{
long i, count = 0;
for (i = 0; i < plsize; i++) {
count += pl[i];
}
return count;
}
int grib_is_earth_oblate(grib_handle* h)
{
long oblate = 0;
int err = grib_get_long(h, "earthIsOblate", &oblate);
if (!err && oblate == 1) {
return 1;
}
return 0;
}
int grib_util_grib_data_quality_check(grib_handle* h, double min_val, double max_val)
{
int err = 0;
double min_field_value_allowed = 0, max_field_value_allowed = 0;
long paramId = 0;
grib_context* ctx = h->context;
int is_error = 1;
char description[1024] = {0,};
char step[32] = "unknown";
size_t len = 32;
/*
* If grib_data_quality_checks == 1, limits failure results in an error
* If grib_data_quality_checks == 2, limits failure results in a warning
*/
Assert(ctx->grib_data_quality_checks == 1 || ctx->grib_data_quality_checks == 2);
is_error = (ctx->grib_data_quality_checks == 1);
/* The limit keys must exist if we are here */
err = grib_get_double(h, "param_value_min", &min_field_value_allowed);
if (err) {
grib_context_log(ctx, GRIB_LOG_ERROR, "grib_data_quality_check: Could not get param_value_min");
return err;
}
err = grib_get_double(h, "param_value_max", &max_field_value_allowed);
if (err) {
grib_context_log(ctx, GRIB_LOG_ERROR, "grib_data_quality_check: Could not get param_value_max");
return err;
}
if (ctx->debug) {
if (get_concept_condition_string(h, "param_value_max", NULL, description) == GRIB_SUCCESS) {
printf("ECCODES DEBUG grib_data_quality_check: Checking condition '%s' (min=%g, max=%g)\n",
description, min_field_value_allowed, max_field_value_allowed);
}
}
if (min_val < min_field_value_allowed) {
grib_get_string(h, "step", step, &len);
if (get_concept_condition_string(h, "param_value_min", NULL, description) == GRIB_SUCCESS) {
fprintf(stderr, "ECCODES %s : (%s, step=%s): minimum (%g) is less than the allowable limit (%g)\n",
(is_error ? "ERROR" : "WARNING"), description, step, min_val, min_field_value_allowed);
}
else {
if (grib_get_long(h, "paramId", &paramId) == GRIB_SUCCESS) {
fprintf(stderr, "ECCODES %s : (paramId=%ld, step=%s): minimum (%g) is less than the default allowable limit (%g)\n",
(is_error ? "ERROR" : "WARNING"), paramId, step, min_val, min_field_value_allowed);
}
}
if (is_error) {
return GRIB_OUT_OF_RANGE; /* Failure */
}
}
if (max_val > max_field_value_allowed) {
grib_get_string(h, "step", step, &len);
if (get_concept_condition_string(h, "param_value_max", NULL, description) == GRIB_SUCCESS) {
fprintf(stderr, "ECCODES %s : (%s, step=%s): maximum (%g) is more than the allowable limit (%g)\n",
(is_error ? "ERROR" : "WARNING"), description, step, max_val, max_field_value_allowed);
}
else {
if (grib_get_long(h, "paramId", &paramId) == GRIB_SUCCESS) {
fprintf(stderr, "ECCODES %s : (paramId=%ld, step=%s): maximum (%g) is more than the default allowable limit (%g)\n",
(is_error ? "ERROR" : "WARNING"), paramId, step, max_val, max_field_value_allowed);
}
}
if (is_error) {
return GRIB_OUT_OF_RANGE; /* Failure */
}
}
return GRIB_SUCCESS;
}
Reference in New Issue View Git Blame Copy Permalink