eccodes/tests/grib_util_set_spec.cc

/*
 * (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 "eccodes.h"

static int get_packing_type_code(const char* packingType)
{
    int result = GRIB_UTIL_PACKING_TYPE_GRID_SIMPLE;
    if (packingType == NULL)
        return GRIB_UTIL_PACKING_TYPE_SAME_AS_INPUT;

    if (STR_EQUAL(packingType, "grid_jpeg"))
        return GRIB_UTIL_PACKING_TYPE_JPEG;
    else if (STR_EQUAL(packingType, "grid_ccsds"))
        return GRIB_UTIL_PACKING_TYPE_CCSDS;
    else if (STR_EQUAL(packingType, "grid_simple"))
        return GRIB_UTIL_PACKING_TYPE_GRID_SIMPLE;
    else if (STR_EQUAL(packingType, "grid_second_order"))
        return GRIB_UTIL_PACKING_TYPE_GRID_SECOND_ORDER;
    else if (STR_EQUAL(packingType, "grid_ieee"))
        return GRIB_UTIL_PACKING_TYPE_IEEE;
    else if (STR_EQUAL(packingType, "grid_complex"))
        return GRIB_UTIL_PACKING_TYPE_GRID_COMPLEX;

    Assert(!"Invalid packingType");
    return result;
}

static void test_reduced_gg(int remove_local_def, int edition, const char* packingType,
                            const char* input_filename, const char* output_filename)
{
    /* based on copy_spec_from_ksec */
    int err     = 0;
    size_t slen = 32, inlen = 0, outlen = 0;
    size_t i = 0, size = 0;
    int set_spec_flags = 0;
    double* values     = NULL;
    FILE* in           = NULL;
    FILE* out          = NULL;
    const void* buffer = NULL;
    char gridType[128] = {0,};

    grib_handle* handle     = 0;
    grib_handle* finalh     = 0;
    grib_util_grid_spec spec = {0,};
    grib_util_packing_spec packing_spec = {0,};

    Assert(input_filename);
    in = fopen(input_filename, "rb");
    Assert(in);
    handle = grib_handle_new_from_file(0, in, &err);
    Assert(handle);

    CODES_CHECK(grib_get_string(handle, "gridType", gridType, &slen), 0);
    if (!STR_EQUAL(gridType, "reduced_gg")) {
        grib_handle_delete(handle);
        return;
    }
    Assert(output_filename);
    out = fopen(output_filename, "wb");
    Assert(out);

    CODES_CHECK(grib_get_size(handle, "values", &inlen), 0);
    values = (double*)malloc(sizeof(double) * inlen);
    CODES_CHECK(grib_get_double_array(handle, "values", values, &inlen), 0);
    // make sure values are not constant
    values[0] = 4.4;
    values[1] = 5.5;
    for (i = 0; i < inlen; ++i) {
        values[i] *= 1.10;
    }

    spec.grid_type                          = GRIB_UTIL_GRID_SPEC_REDUCED_GG;
    spec.N                                  = 32; /* hardcoded for now */
    spec.Nj                                 = 2 * spec.N;
    outlen                                  = inlen;
    spec.iDirectionIncrementInDegrees       = 1.5;
    spec.jDirectionIncrementInDegrees       = 1.5;
    spec.latitudeOfFirstGridPointInDegrees  = 87.8637991;
    spec.longitudeOfFirstGridPointInDegrees = 0.0;
    spec.latitudeOfLastGridPointInDegrees   = -spec.latitudeOfFirstGridPointInDegrees;
    spec.longitudeOfLastGridPointInDegrees  = 357.187500;
    spec.bitmapPresent                      = 0;

    /*packing_spec.packing_type = GRIB_UTIL_PACKING_TYPE_GRID_SECOND_ORDER;*/
    packing_spec.packing_type = get_packing_type_code(packingType);
    packing_spec.bitsPerValue = 24;
    packing_spec.accuracy     = GRIB_UTIL_ACCURACY_USE_PROVIDED_BITS_PER_VALUES;
    if (packingType)
        packing_spec.packing  = GRIB_UTIL_PACKING_USE_PROVIDED;
    else
        packing_spec.packing  = GRIB_UTIL_PACKING_SAME_AS_INPUT;

    /*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,
        &spec,
        &packing_spec,
        set_spec_flags,
        values,
        outlen,
        &err);
    Assert(finalh);
    Assert(err == 0);

    /* Try some invalid inputs and check it is handled */
    {
        grib_handle* h2      = 0;
        packing_spec.accuracy = 999;
        h2                    = grib_util_set_spec(handle, &spec, &packing_spec, set_spec_flags, values, outlen, &err);
        Assert(err == GRIB_INTERNAL_ERROR);
        Assert(!h2);
        if (h2) exit(1);
#ifdef INFINITY
        packing_spec.accuracy = GRIB_UTIL_ACCURACY_USE_PROVIDED_BITS_PER_VALUES;
        values[0]             = INFINITY;
        h2                    = grib_util_set_spec(handle, &spec, &packing_spec, set_spec_flags, values, outlen, &err);
        Assert(err == GRIB_ENCODING_ERROR);
        Assert(!h2);
        if (h2) exit(1);

        values[0]             = -INFINITY;
        h2                    = grib_util_set_spec(handle, &spec, &packing_spec, set_spec_flags, values, outlen, &err);
        Assert(err == GRIB_ENCODING_ERROR);
        Assert(!h2);
        if (h2) exit(1);
#endif
    }

    /* Write out the message to the output file */
    CODES_CHECK(grib_get_message(finalh, &buffer, &size), 0);
    CODES_CHECK(codes_check_message_header(buffer, size, PRODUCT_GRIB), 0);
    CODES_CHECK(codes_check_message_footer(buffer, size, PRODUCT_GRIB), 0);
    if (fwrite(buffer, 1, size, out) != size) {
        Assert(0);
    }
    grib_handle_delete(handle);
    grib_handle_delete(finalh);
    fclose(in);
    fclose(out);
    free(values);
    /*printf("ALL OK: %s\n", __func__);*/
}

static void test_regular_ll(int remove_local_def, int edition, const char* packingType,
                            const char* input_filename, const char* output_filename)
{
    /* based on copy_spec_from_ksec */
    int err     = 0;
    size_t slen = 32, inlen = 0, outlen = 0;
    size_t size        = 0;
    int set_spec_flags = 0;
    double* values     = NULL;
    FILE* in           = NULL;
    FILE* out          = NULL;
    const void* buffer = NULL;
    char gridType[128] = {0,};
    long input_edition = 0;

    grib_handle* handle     = 0;
    grib_handle* finalh     = 0;
    grib_util_grid_spec spec = {0,};
    grib_util_packing_spec packing_spec = {0,};

    Assert(input_filename);
    in = fopen(input_filename, "rb");
    Assert(in);
    handle = grib_handle_new_from_file(0, in, &err);
    Assert(handle);

    CODES_CHECK(grib_get_long(handle, "edition", &input_edition), 0);

    CODES_CHECK(grib_get_string(handle, "gridType", gridType, &slen), 0);
    if (!STR_EQUAL(gridType, "regular_ll")) {
        grib_handle_delete(handle);
        return;
    }
    Assert(output_filename);
    out = fopen(output_filename, "wb");
    Assert(out);

    CODES_CHECK(grib_get_size(handle, "values", &inlen), 0);
    values = (double*)malloc(sizeof(double) * inlen);
    CODES_CHECK(grib_get_double_array(handle, "values", values, &inlen), 0);
    // make sure values are not constant
    values[0] = 4.4;
    values[1] = 5.5;

    spec.grid_type = GRIB_UTIL_GRID_SPEC_REGULAR_LL;
    spec.Nj        = 14;
    spec.Ni        = 17;
    outlen         = spec.Nj * spec.Ni;
    /* outlen = inlen; */
    spec.iDirectionIncrementInDegrees       = 1.5;
    spec.jDirectionIncrementInDegrees       = 1.5;
    spec.latitudeOfFirstGridPointInDegrees  = 60.0001;
    spec.longitudeOfFirstGridPointInDegrees = -9.0;
    spec.latitudeOfLastGridPointInDegrees   = 40.5;
    spec.longitudeOfLastGridPointInDegrees  = 15.0;
    spec.bitmapPresent                      = 0;

    /*packing_spec.packing_type = GRIB_UTIL_PACKING_TYPE_GRID_SIMPLE;*/
    packing_spec.packing_type = get_packing_type_code(packingType);
    packing_spec.bitsPerValue = 24;
    packing_spec.accuracy     = GRIB_UTIL_ACCURACY_USE_PROVIDED_BITS_PER_VALUES;
    if (packingType)
        packing_spec.packing  = GRIB_UTIL_PACKING_USE_PROVIDED;
    else
        packing_spec.packing  = GRIB_UTIL_PACKING_SAME_AS_INPUT;

    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;

    if (edition != 0) {
        packing_spec.editionNumber = edition;
    }
    if (remove_local_def) {
        packing_spec.deleteLocalDefinition = 1;
    }

    finalh = grib_util_set_spec(
        handle,
        &spec,
        &packing_spec,
        set_spec_flags,
        values,
        outlen,
        &err);
    Assert(finalh);
    Assert(err == 0);

    /* Check expand_bounding_box worked.
     * Specified latitudeOfFirstGridPointInDegrees cannot be encoded in GRIB1
     * so it is rounded up to nearest millidegree */
    if (input_edition == 1) {
        const double expected_lat1 = 60.001;
        double lat1                = 0;
        CODES_CHECK(grib_get_double(finalh, "latitudeOfFirstGridPointInDegrees", &lat1), 0);
        Assert(fabs(lat1 - expected_lat1) < 1e-10);
    }

    /* Write out the message to the output file */
    CODES_CHECK(grib_get_message(finalh, &buffer, &size), 0);
    CODES_CHECK(codes_check_message_header(buffer, size, PRODUCT_GRIB), 0);
    CODES_CHECK(codes_check_message_footer(buffer, size, PRODUCT_GRIB), 0);
    if (fwrite(buffer, 1, size, out) != size) {
        Assert(0);
    }
    grib_handle_delete(handle);
    grib_handle_delete(finalh);
    free(values);
    fclose(in);
    fclose(out);
}

#if 0
static void test_grid_complex_spatial_differencing(int remove_local_def, int edition, const char* packingType,
                     const char* input_filename, const char* output_filename)
{
    /* based on copy_spec_from_ksec */
    int err = 0;
    size_t slen = 34, inlen = 0, outlen = 0;
    size_t size=0;
    int set_spec_flags=0;
    double* values = NULL;
    FILE* in = NULL;
    FILE* out = NULL;
    const void* buffer = NULL;
    char gridType[128] = {0,};
    double theMax,theMin,theAverage;

    grib_handle *handle = 0;
    grib_handle *finalh = 0;
    grib_util_grid_spec spec={0,};
    grib_util_packing_spec packing_spec={0,};

    in = fopen(input_filename,"rb");     Assert(in);
    handle = grib_handle_new_from_file(0, in, &err);    Assert(handle);

    CODES_CHECK(grib_get_string(handle, "packingType", gridType, &slen),0);
    if (!STR_EQUAL(gridType, "grid_complex_spatial_differencing")) {
        grib_handle_delete(handle);
        return;
    }
    out = fopen(output_filename,"wb");   Assert(out);

    CODES_CHECK(grib_get_size(handle,"values",&inlen), 0);
    values = (double*)malloc(sizeof(double)*inlen);
    CODES_CHECK(grib_get_double_array(handle, "values", values,&inlen), 0);

    CODES_CHECK(grib_get_double(handle, "max",    &theMax),0);
    CODES_CHECK(grib_get_double(handle, "min",    &theMin),0);
    CODES_CHECK(grib_get_double(handle, "average",&theAverage),0);
    printf("inlen=%lu \t max=%g \t min=%g \t avg=%g\n", inlen, theMax, theMin, theAverage);

    spec.grid_type = GRIB_UTIL_GRID_SPEC_REGULAR_LL;
    spec.Nj = 721;
    spec.Ni = 1440;
    outlen = spec.Nj * spec.Ni;
    spec.iDirectionIncrementInDegrees = 0.25;
    spec.jDirectionIncrementInDegrees = 0.25;
    spec.latitudeOfFirstGridPointInDegrees  = 90.0;
    spec.longitudeOfFirstGridPointInDegrees = 0.0;
    spec.latitudeOfLastGridPointInDegrees   = -90.0;
    spec.longitudeOfLastGridPointInDegrees  = 359.75;
    spec.bitmapPresent = 1; /* there are missing values inside the data section! */
    spec.missingValue = 9999;

    packing_spec.packing_type = GRIB_UTIL_PACKING_TYPE_GRID_SIMPLE; /*Convert to Grid Simple Packing*/
    packing_spec.bitsPerValue = 11;
    packing_spec.accuracy=GRIB_UTIL_ACCURACY_USE_PROVIDED_BITS_PER_VALUES;
    /*packing_spec.packing=GRIB_UTIL_PACKING_USE_PROVIDED;*/

    if (edition != 0) {
        packing_spec.editionNumber = edition;
    }
    if (remove_local_def) {
        packing_spec.deleteLocalDefinition = 1;
    }

    finalh = grib_util_set_spec(
            handle,
            &spec,
            &packing_spec,
            set_spec_flags,
            values,
            outlen,
            &err);
    Assert(finalh);
    Assert(err == 0);

    /* Write out the message to the output file */
    CODES_CHECK(grib_get_message(finalh, &buffer, &size),0);
    if(fwrite(buffer,1,size,out) != size) {
        Assert(0);
    }
    grib_handle_delete(handle);
    grib_handle_delete(finalh);
    fclose(in);
    fclose(out);
}
#endif

// static void usage(const char* prog)
// {
//     fprintf(stderr, "%s: [-p packingType] [-r] [-e edition] in.grib out.grib\n", prog);
//     fprintf(stderr, "-p  packingType: one of grid_jpeg, grid_ccsds, grid_second_order or grid_simple\n");
//     fprintf(stderr, "-r  remove local definition\n");
//     fprintf(stderr, "-e  edition: 1 or 2\n");
//     exit(1);
// }

int main(int argc, char* argv[])
{
    int i = 0, remove_local_def = 0;
    int edition        = 0;
    char* packingType  = NULL;
    char* infile_name  = NULL;
    char* outfile_name = NULL;

    if (argc == 1 || argc > 8) return 1;// see usage

    for (i = 1; i < argc; i++) {
        if (strcmp(argv[i], "-p") == 0) {
            packingType = argv[i + 1];
            ++i;
        }
        else if (strcmp(argv[i], "-e") == 0) {
            edition = atoi(argv[i + 1]);
            ++i;
        }
        else if (strcmp(argv[i], "-r") == 0) {
            remove_local_def = 1;
        }
        else {
            /* Expect 2 filenames */
            infile_name  = argv[i];
            outfile_name = argv[i + 1];
            break;
        }
    }
#if 0
    printf("DEBUG remove_local_def = %d\n", remove_local_def);
    printf("DEBUG edition          = %d\n", edition);
    printf("DEBUG packingType      = %s\n", packingType);
    printf("DEBUG infile_name      = %s\n", infile_name);
    printf("DEBUG outfile_name     = %s\n", outfile_name);
#endif
    test_regular_ll(remove_local_def, edition, packingType, infile_name, outfile_name);
    test_reduced_gg(remove_local_def, edition, packingType, infile_name, outfile_name);
    /*test_grid_complex_spatial_differencing(remove_local_def, edition, packingType, infile_name, outfile_name);*/

    return 0;
}