eccodes/src/grib_iterator_class_latlon.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"

/*
   This is used by make_class.pl

   START_CLASS_DEF
   CLASS      = iterator
   SUPER      = grib_iterator_class_regular
   IMPLEMENTS = init;next
   END_CLASS_DEF

 */

/* START_CLASS_IMP */

/*

Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
Instead edit values between START_CLASS_DEF and END_CLASS_DEF
or edit "iterator.class" and rerun ./make_class.pl

*/


static void init_class              (grib_iterator_class*);

static int init               (grib_iterator* i,grib_handle*,grib_arguments*);
static int next               (grib_iterator* i, double *lat, double *lon, double *val);


typedef struct grib_iterator_latlon{
  grib_iterator it;
    /* Members defined in gen */
    int carg;
    const char* missingValue;
    /* Members defined in regular */
    double   *las;
    double   *los;
    long      Ni;
    long      Nj;
    long iScansNegatively;
    long isRotated;
    double angleOfRotation;
    double southPoleLat;
    double southPoleLon;
    long jPointsAreConsecutive;
    long disableUnrotate;
    /* Members defined in latlon */
} grib_iterator_latlon;

extern grib_iterator_class* grib_iterator_class_regular;

static grib_iterator_class _grib_iterator_class_latlon = {
    &grib_iterator_class_regular,                    /* super                     */
    "latlon",                    /* name                      */
    sizeof(grib_iterator_latlon),/* size of instance          */
    0,                           /* inited */
    &init_class,                 /* init_class */
    &init,                     /* constructor               */
    0,                  /* destructor                */
    &next,                     /* Next Value                */
    0,                 /*  Previous Value           */
    0,                    /* Reset the counter         */
    0,                 /* has next values           */
};

grib_iterator_class* grib_iterator_class_latlon = &_grib_iterator_class_latlon;


static void init_class(grib_iterator_class* c)
{
    c->previous    =    (*(c->super))->previous;
    c->reset    =    (*(c->super))->reset;
    c->has_next    =    (*(c->super))->has_next;
}
/* END_CLASS_IMP */

static int next(grib_iterator* iter, double* lat, double* lon, double* val)
{
    /* GRIB-238: Support rotated lat/lon grids */
    double ret_lat, ret_lon, ret_val=0;
    grib_iterator_latlon* self = (grib_iterator_latlon*)iter;

    if ((long)iter->e >= (long)(iter->nv - 1))
        return 0;

    iter->e++;

    /* Assumptions:
     *   All rows scan in the same direction (alternativeRowScanning==0)
     */
    if (!self->jPointsAreConsecutive) {
        /* Adjacent points in i (x) direction are consecutive */
        ret_lat = self->las[(long)floor(iter->e / self->Ni)];
        ret_lon = self->los[(long)iter->e % self->Ni];
        if (iter->data)
            ret_val = iter->data[iter->e];
    }
    else {
        /* Adjacent points in j (y) direction is consecutive */
        ret_lon = self->los[(long)iter->e / self->Nj];
        ret_lat = self->las[(long)floor(iter->e % self->Nj)];
        if (iter->data)
            ret_val = iter->data[iter->e];
    }

    /* See ECC-808: Some users want to disable the unrotate */
    if (self->isRotated && !self->disableUnrotate) {
        double new_lat = 0, new_lon = 0;
        unrotate(ret_lat, ret_lon,
                 self->angleOfRotation, self->southPoleLat, self->southPoleLon,
                 &new_lat, &new_lon);
        ret_lat = new_lat;
        ret_lon = new_lon;
    }

    *lat = ret_lat;
    *lon = ret_lon;
    if (val && iter->data) {
        *val = ret_val;
    }
    return 1;
}

static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
{
    grib_iterator_latlon* self = (grib_iterator_latlon*)iter;
    int err                    = 0;
    double jdir;
    double lat1=0, lat2=0, north=0, south=0;
    long jScansPositively;
    long lai;

    const char* s_lat1            = grib_arguments_get_name(h, args, self->carg++);
    const char* s_jdir            = grib_arguments_get_name(h, args, self->carg++);
    const char* s_jScansPos       = grib_arguments_get_name(h, args, self->carg++);
    const char* s_jPtsConsec      = grib_arguments_get_name(h, args, self->carg++);
    const char* s_isRotatedGrid   = grib_arguments_get_name(h, args, self->carg++);
    const char* s_angleOfRotation = grib_arguments_get_name(h, args, self->carg++);
    const char* s_latSouthernPole = grib_arguments_get_name(h, args, self->carg++);
    const char* s_lonSouthernPole = grib_arguments_get_name(h, args, self->carg++);

    self->angleOfRotation  = 0;
    self->isRotated        = 0;
    self->southPoleLat     = 0;
    self->southPoleLon     = 0;
    self->disableUnrotate  = 0; /* unrotate enabled by default */

    if ((err = grib_get_long(h, s_isRotatedGrid, &self->isRotated)))
        return err;
    if (self->isRotated) {
        if ((err = grib_get_double_internal(h, s_angleOfRotation, &self->angleOfRotation)))
            return err;
        if ((err = grib_get_double_internal(h, s_latSouthernPole, &self->southPoleLat)))
            return err;
        if ((err = grib_get_double_internal(h, s_lonSouthernPole, &self->southPoleLon)))
            return err;
    }

    if ((err = grib_get_double_internal(h, s_lat1, &lat1)))
        return err;
    if ((err = grib_get_double_internal(h, "latitudeLastInDegrees", &lat2)))
        return err;
    if ((err = grib_get_double_internal(h, s_jdir, &jdir))) //can be GRIB_MISSING_DOUBLE
        return err;
    if ((err = grib_get_long_internal(h, s_jScansPos, &jScansPositively)))
        return err;
    if ((err = grib_get_long_internal(h, s_jPtsConsec, &self->jPointsAreConsecutive)))
        return err;
    if ((err = grib_get_long(h, "iteratorDisableUnrotate", &self->disableUnrotate)))
        return err;

    /* ECC-984: If jDirectionIncrement is missing, then we cannot use it (See jDirectionIncrementGiven) */
    /* So try to compute the increment */
    if ( (grib_is_missing(h, s_jdir, &err) && err == GRIB_SUCCESS) || (jdir == GRIB_MISSING_DOUBLE) ) {
        const long Nj = self->Nj;
        Assert(Nj > 1);
        if (lat1 > lat2) {
            jdir = (lat1 - lat2) / (Nj - 1);
        }
        else {
            jdir = (lat2 - lat1) / (Nj - 1);
        }
        grib_context_log(h->context, GRIB_LOG_DEBUG,
                        "Cannot use jDirectionIncrement. Using value of %.6f obtained from La1, La2 and Nj", jdir);
    }

    if (jScansPositively) {
        north = lat2;
        south = lat1;
        jdir = -jdir;
    } else {
        north = lat1;
        south = lat2;
    }
    if (south > north) {
        grib_context_log(h->context, GRIB_LOG_ERROR,
                         "Lat/Lon Geoiterator: First and last latitudes are inconsistent with scanning order: lat1=%g, lat2=%g jScansPositively=%ld",
                         lat1, lat2, jScansPositively);
        return GRIB_WRONG_GRID;
    }

    for (lai = 0; lai < self->Nj; lai++) {
        self->las[lai] = lat1;
        lat1 -= jdir;
    }
    /* ECC-1406: Due to rounding, errors can accumulate.
     * So we ensure the last latitude is latitudeOfLastGridPointInDegrees
    */
    self->las[self->Nj-1] = lat2;

    iter->e = -1;
    return err;
}