eccodes/tools/grib_merge.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.
 */

/*
 * Author: Enrico Fucile
 * Description:
 * In Observations team we need a tool to merge the GRIB messages coming from GTS.
 * They come with compatible grid (lat/lon) but they are split in different quadrants.
 * This tool should be able to merge different quadrants in one bigger grid,
 * automatically recognising when this can be done by checking that the fields have the same
 * parameter, step, ... and compatible grid.
 * The current solution is to have a grib_merge working on an input file and writing in output the
 * merged fields working out automatically what has to be merged.
 */

#include "grib_tools.h"

#define MAX_KEY_VALUES 100

grib_handle* hh = 0;
grib_values key_values[MAX_KEY_VALUES];
int key_values_size = MAX_KEY_VALUES;

/* This key was added to provide a unique ID for the product description used by the merge
 * tool to identify the fields to be merged */
static const char* md5Key = "md5Product";

const char* tool_description = "Merge two fields with identical parameters and different geographical area";
const char* tool_name        = "grib_merge";
const char* tool_online_doc  = NULL;
const char* tool_usage       = "[options] file file ... output_file";

grib_option grib_options[] = {
    /*  {id, args, help}, on, command_line, value */
    { "f", 0, 0, 0, 1, 0 },
    { "c", 0, 0, 1, 0, 0 },
    { "r", 0, 0, 0, 1, 0 },
    { "q", 0, 0, 1, 0, 0 },
    { "p:", 0, 0, 1, 1, 0 },
    { "P:", 0, 0, 0, 1, 0 },
    { "B:", 0, 0, 1, 1, "md5Product" },
    { "V", 0, 0, 0, 1, 0 },
    { "W:", 0, 0, 0, 1, 0 },
    { "M", 0, 0, 0, 1, 0 },
    { "U", 0, 0, 1, 0, 0 },
    { "H", 0, 0, 1, 0, 0 },
    { "T:", 0, 0, 1, 0, "G" },
    { "S", 0, 0, 1, 0, 0 },
    { "g", 0, 0, 0, 1, 0 },
    { "G", 0, 0, 0, 1, 0 },
    { "7", 0, 0, 0, 1, 0 },
    { "v", 0, 0, 0, 1, 0 },
    { "h", 0, 0, 0, 1, 0 },
};

int grib_options_count = sizeof(grib_options) / sizeof(grib_option);

static int grib_key_equal(const grib_handle* h1, const grib_handle* h2, const char* key, int type, int* err)
{
    double d1 = 0, d2 = 0;
    long l1 = 0, l2 = 0;
    char s1[500] = {0,};
    char s2[500] = {0,};
    size_t len1, len2;

    if (type != GRIB_TYPE_DOUBLE &&
        type != GRIB_TYPE_LONG &&
        type != GRIB_TYPE_STRING) {
        *err = grib_get_native_type(h1, key, &type);
    }
    switch (type) {
        case GRIB_TYPE_DOUBLE:
            *err = grib_get_double(h1, key, &d1);
            *err = grib_get_double(h2, key, &d2);
            if (d1 != d2)
                return 0;
            break;
        case GRIB_TYPE_LONG:
            *err = grib_get_long(h1, key, &l1);
            *err = grib_get_long(h2, key, &l2);
            if (l1 != l2)
                return 0;
            break;
        default:
            len1 = sizeof(s1) / sizeof(*s1);
            len2 = sizeof(s2) / sizeof(*s2);
            *err = grib_get_string(h1, key, s1, &len1);
            *err = grib_get_string(h2, key, s2, &len2);
            if (!STR_EQUAL(s1, s2))
                return 0;
            break;
    }
    return 1;
}

int main(int argc, char* argv[])
{
    return grib_tool(argc, argv);
}

int grib_tool_before_getopt(grib_runtime_options* options)
{
    return 0;
}

int grib_tool_init(grib_runtime_options* options)
{
    return 0;
}

int grib_tool_new_filename_action(grib_runtime_options* options, const char* file)
{
    return 0;
}

int grib_tool_new_file_action(grib_runtime_options* options, grib_tools_file* file)
{
    return 0;
}

static int idx(double lat, double lon, double latFirst, double lonFirst, double latLast, double lonLast,
               long Ni, double di, double dj)
{
    long ilon, ilat;
    if ((ilon = (lon - lonFirst) / di) < 0)
        return -1;
    if ((ilat = (latFirst - lat) / dj) < 0)
        return -1;
    if (lon > lonLast) {
        if (lonLast == 180) {
            lon -= 360;
            ilon = (lon - lonFirst) / di;
        }
        else
            return -1;
    }
    if (lat < latLast)
        return -1;
    if ((ilat = (latFirst - lat) / dj) < 0)
        return -1;
    return ilon + ilat * Ni;
}

static grib_handle* merge(grib_handle* h1, grib_handle* h2)
{
    char s1[100] = {0,};
    size_t len1;
    char s2[100] = {0,};
    size_t len2;
    size_t sn;
    long l1, l2;
    double di, dj, di1, dj1, di2, dj2;
    long bitmapPresent;
    long n = 0, n1, n2;
    double latFirst1, latFirst2, latFirst;
    double latLast1, latLast2, latLast;
    double lonFirst1, lonFirst2, lonFirst;
    double lonLast1, lonLast2, lonLast;
    double missingValue;
    double *v = 0, *v1 = 0, *v2 = 0, *lat1 = 0, *lat2 = 0, *lon1 = 0, *lon2 = 0, *lat = 0, *lon = 0;
    long i, j, iscan, jscan, Ni, Nj, idj, idi;
    long Ni1, Nj1, Ni2, Nj2;
    grib_handle* h = NULL;
    int err        = 0;

    /* Same products? */
    if (grib_key_equal(h1, h2, md5Key, GRIB_TYPE_STRING, &err) == 0 && err == 0) {
        return NULL;
    }

    /* Can we do it? */
    len2 = sizeof(s2) / sizeof(*s2);
    err  = grib_get_string(h2, "gridType", s2, &len2);
    if (err) return NULL;
    if (strcmp(s2, "regular_ll")) {
        grib_context_log(h1->context, GRIB_LOG_ERROR, "gridType=%s not supported", s2);
        return NULL;
    }

    len1 = sizeof(s1) / sizeof(*s1);
    err  = grib_get_string(h1, "gridType", s1, &len1);
    if (err) return NULL;
    if (strcmp(s1, "regular_ll")) {
        grib_context_log(h1->context, GRIB_LOG_ERROR, "gridType=%s not supported", s1);
        return NULL;
    }

    if (!grib_key_equal(h1, h2, "iDirectionIncrementInDegrees", GRIB_TYPE_DOUBLE, &err)) {
        grib_context_log(h1->context, GRIB_LOG_ERROR,
                         "unable to merge: different iDirectionIncrementInDegrees");
        return NULL;
    }
    if (!grib_key_equal(h1, h2, "jDirectionIncrementInDegrees", GRIB_TYPE_DOUBLE, &err)) {
        grib_context_log(h1->context, GRIB_LOG_ERROR,
                         "unable to merge: different jDirectionIncrementInDegrees");
        return NULL;
    }

    grib_get_long(h1, "latitudeOfFirstGridPoint", &l1);
    grib_get_long(h2, "latitudeOfFirstGridPoint", &l2);
    grib_get_long(h2, "jDirectionIncrement", &idj);
    if ((l1 - l2) % idj) {
        grib_context_log(h1->context, GRIB_LOG_ERROR, "unable to merge: incompatible grid");
        return NULL;
    }

    grib_get_long(h1, "longitudeOfFirstGridPoint", &l1);
    grib_get_long(h2, "longitudeOfFirstGridPoint", &l2);
    grib_get_long(h2, "iDirectionIncrement", &idi);
    if ((l1 - l2) % idi) {
        grib_context_log(h1->context, GRIB_LOG_ERROR, "unable to merge: incompatible grid");
        return NULL;
    }

    /* Yes we can! Do we have something to do? */
    if (grib_key_equal(h1, h2, "latitudeOfFirstGridPointInDegrees", GRIB_TYPE_DOUBLE, &err) &&
        grib_key_equal(h1, h2, "latitudeOfLastGridPointInDegrees", GRIB_TYPE_DOUBLE, &err) &&
        grib_key_equal(h1, h2, "longitudeOfFirstGridPointInDegrees", GRIB_TYPE_DOUBLE, &err) &&
        grib_key_equal(h1, h2, "longitudeOfLastGridPointInDegrees", GRIB_TYPE_DOUBLE, &err)) {
        /* No we don't */
        return NULL;
    }

    /* Yes we do! */
    /* Check scanning mode */
    grib_get_long(h1, "iScansNegatively", &iscan);
    if (iscan)
        grib_set_long(h1, "swapScanningLon", 1);
    grib_get_long(h2, "iScansNegatively", &iscan);
    if (iscan)
        grib_set_long(h2, "swapScanningLon", 1);

    grib_get_long(h1, "jScansPositively", &jscan);
    if (jscan)
        grib_set_long(h1, "swapScanningLat", 1);
    grib_get_long(h2, "jScansPositively", &jscan);
    if (jscan)
        grib_set_long(h2, "swapScanningLat", 1);

    grib_get_double(h1, "latitudeOfFirstGridPointInDegrees", &latFirst1);
    grib_get_double(h2, "latitudeOfFirstGridPointInDegrees", &latFirst2);
    grib_get_double(h1, "longitudeOfFirstGridPointInDegrees", &lonFirst1);
    grib_get_double(h2, "longitudeOfFirstGridPointInDegrees", &lonFirst2);

    grib_get_double(h1, "latitudeOfLastGridPointInDegrees", &latLast1);
    grib_get_double(h2, "latitudeOfLastGridPointInDegrees", &latLast2);
    grib_get_double(h1, "longitudeOfLastGridPointInDegrees", &lonLast1);
    grib_get_double(h2, "longitudeOfLastGridPointInDegrees", &lonLast2);

    /* ECC-949 */
    if (lonFirst1 == 180 && lonLast1 == 180) {
        lonFirst1 = -180;
    }
    if (lonFirst2 == 180 && lonLast2 == 180) {
        lonFirst2 = -180;
    }

    latFirst = latFirst1 > latFirst2 ? latFirst1 : latFirst2;
    lonFirst = lonFirst1 < lonFirst2 ? lonFirst1 : lonFirst2;

    latLast = latLast1 < latLast2 ? latLast1 : latLast2;
    lonLast = lonLast1 > lonLast2 ? lonLast1 : lonLast2;

    grib_get_double(h1, "iDirectionIncrementInDegrees", &di1);
    grib_get_double(h1, "jDirectionIncrementInDegrees", &dj1);
    grib_get_double(h2, "iDirectionIncrementInDegrees", &di2);
    grib_get_double(h2, "jDirectionIncrementInDegrees", &dj2);
    grib_get_double(h1, "iDirectionIncrementInDegrees", &di);
    grib_get_double(h1, "jDirectionIncrementInDegrees", &dj);
    grib_get_long(h1, "Ni", &Ni1);
    grib_get_long(h1, "Nj", &Nj1);
    grib_get_long(h2, "Ni", &Ni2);
    grib_get_long(h2, "Nj", &Nj2);

    if (lonFirst == 0 && lonLast == 360)
        lonLast -= di;
    if (lonFirst == -180 && lonLast == 180) {
        lonFirst = 0;
        lonLast  = 360 - di;
    }

    /* Create new grib for bigger area*/
    h = grib_handle_clone(h1);
    grib_set_double(h, "latitudeOfFirstGridPointInDegrees", latFirst);
    grib_set_double(h, "longitudeOfFirstGridPointInDegrees", lonFirst);
    grib_set_double(h, "latitudeOfLastGridPointInDegrees", latLast);
    grib_set_double(h, "longitudeOfLastGridPointInDegrees", lonLast);
    Ni = fabs(lonLast - lonFirst) / di + 1;
    Nj = fabs(latLast - latFirst) / dj + 1;

    grib_set_long(h, "Ni", Ni);
    grib_set_long(h, "Nj", Nj);

    grib_get_long(h, "bitmapPresent", &bitmapPresent);
    if (!bitmapPresent)
        grib_set_long(h, "bitmapPresent", 1);

    grib_get_long(h, "numberOfPoints", &n);
    grib_get_double(h, "missingValue", &missingValue);
    v = (double*)grib_context_malloc_clear(h->context, sizeof(double) * n);
    for (i = 0; i < n; i++)
        v[i] = 1.0;
    sn = n;
    grib_set_double_array(h, "values", v, sn);
    for (i = 0; i < n; i++)
        v[i] = missingValue;
    lat = (double*)grib_context_malloc_clear(h->context, sizeof(double) * n);
    lon = (double*)grib_context_malloc_clear(h->context, sizeof(double) * n);
    sn  = n;
    grib_get_double_array(h, "latitudes", lat, &sn);
    grib_get_double_array(h, "longitudes", lon, &sn);

    grib_get_long(h1, "numberOfPoints", &n1);
    v1   = (double*)grib_context_malloc_clear(h->context, sizeof(double) * n1);
    lat1 = (double*)grib_context_malloc_clear(h->context, sizeof(double) * n1);
    lon1 = (double*)grib_context_malloc_clear(h->context, sizeof(double) * n1);
    sn   = n1;
    grib_get_double_array(h1, "latitudes", lat1, &sn);
    grib_get_double_array(h1, "longitudes", lon1, &sn);
    grib_get_double_array(h1, "values", v1, &sn);

    grib_get_long(h2, "numberOfPoints", &n2);
    v2   = (double*)grib_context_malloc_clear(h->context, sizeof(double) * n2);
    lat2 = (double*)grib_context_malloc_clear(h->context, sizeof(double) * n2);
    lon2 = (double*)grib_context_malloc_clear(h->context, sizeof(double) * n2);
    sn   = n2;
    grib_get_double_array(h2, "latitudes", lat2, &sn);
    grib_get_double_array(h2, "longitudes", lon2, &sn);
    grib_get_double_array(h2, "values", v2, &sn);

    for (i = 0; i < n; i++) {
        if ((j = idx(lat[i], lon[i], latFirst1, lonFirst1, latLast1, lonLast1, Ni1, di1, dj1)) >= 0) {
            v[i] = v1[j];
        }
        else if ((j = idx(lat[i], lon[i], latFirst2, lonFirst2, latLast2, lonLast2, Ni2, di2, dj2)) >= 0) {
            v[i] = v2[j];
        }
    }

    grib_set_double_array(h, "values", v, n);

    grib_context_free(h->context, v);
    grib_context_free(h->context, v1);
    grib_context_free(h->context, v2);
    grib_context_free(h->context, lat);
    grib_context_free(h->context, lat1);
    grib_context_free(h->context, lat2);
    grib_context_free(h->context, lon);
    grib_context_free(h->context, lon1);
    grib_context_free(h->context, lon2);

    return h;
}

int grib_tool_new_handle_action(grib_runtime_options* options, grib_handle* h)
{
    int err         = 0;
    grib_handle* hm = 0;
    char md5[200]   = {0,};
    char fname[256] = {0,};
    size_t lmd5 = 32;

    if (!hh) {
        hh = grib_handle_clone(h);
        return 0;
    }
    err = grib_get_string(h, md5Key, md5, &lmd5);
    if (err) {
        fprintf(stderr, "grib_merge error getting key %s: %s\n",
                md5Key, grib_get_error_message(err));
        exit(err);
    }
    snprintf(fname, sizeof(fname), "_%s.orig.grib", md5);
    err = grib_write_message(h, fname, "a");

    if ((hm = merge(h, hh)) == NULL) {
        grib_tools_write_message(options, hh);
        lmd5 = sizeof(md5) / sizeof(*md5);
        err  = grib_get_string(hh, md5Key, md5, &lmd5);
        if (err) {
            fprintf(stderr, "grib_merge error getting key %s: %s\n",
                    md5Key, grib_get_error_message(err));
            exit(err);
        }
        snprintf(fname, sizeof(fname), "_%s.merge.grib", md5);
        grib_write_message(hh, fname, "a");
    }
    grib_handle_delete(hh);
    hh = hm != NULL ? hm : grib_handle_clone(h);

    return err;
}

int grib_tool_skip_handle(grib_runtime_options* options, grib_handle* h)
{
    grib_handle_delete(h);
    return 0;
}

int grib_tool_finalise_action(grib_runtime_options* options)
{
    grib_tools_write_message(options, hh);
    if (options->outfile->file) {
        fclose(options->outfile->file);
    }
    return 0;
}

int grib_no_handle_action(grib_runtime_options* options, int err)
{
    fprintf(dump_file, "\t\t\"ERROR: unreadable message\"\n");
    return 0;
}