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

Modernisation: C++ iterators

This commit is contained in:
Eugen Betke 2024-09-21 17:54:46 +02:00
parent 859399da2c
commit 1afda19f17
35 changed files with 1366 additions and 1949 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/CMakeLists.txt
View File

@ -331,20 +331,20 @@ list( APPEND eccodes_src_files
grib_nearest_class_mercator.cc
grib_nearest_class_polar_stereographic.cc
grib_nearest_class_space_view.cc
grib_iterator_class_polar_stereographic.cc
grib_iterator_class_lambert_azimuthal_equal_area.cc
grib_iterator_class_lambert_conformal.cc
grib_iterator_class_mercator.cc
grib_iterator.cc
iterator/grib_iterator_class_polar_stereographic.cc
iterator/grib_iterator_class_lambert_azimuthal_equal_area.cc
iterator/grib_iterator_class_lambert_conformal.cc
iterator/grib_iterator_class_mercator.cc
iterator/grib_iterator.cc
grib_iterator_class.cc
grib_iterator_class_gaussian.cc
grib_iterator_class_gaussian_reduced.cc
grib_iterator_class_latlon_reduced.cc
grib_iterator_class_gen.cc
grib_iterator_class_latlon.cc
grib_iterator_class_regular.cc
grib_iterator_class_space_view.cc
grib_iterator_class_healpix.cc
iterator/grib_iterator_class_gaussian.cc
iterator/grib_iterator_class_gaussian_reduced.cc
iterator/grib_iterator_class_latlon_reduced.cc
iterator/grib_iterator_class_gen.cc
iterator/grib_iterator_class_latlon.cc
iterator/grib_iterator_class_regular.cc
iterator/grib_iterator_class_space_view.cc
iterator/grib_iterator_class_healpix.cc
grib_expression.cc
codes_util.cc
grib_util.cc

12
src/eccodes_prototypes.h
View File

@ -811,12 +811,12 @@ grib_nearest* grib_nearest_factory(grib_handle* h, grib_arguments* args, int* er
/* grib_iterator.cc */
int grib_get_data(const grib_handle* h, double* lats, double* lons, double* values);
int grib_iterator_next(grib_iterator* i, double* lat, double* lon, double* value);
int grib_iterator_has_next(grib_iterator* i);
int grib_iterator_previous(grib_iterator* i, double* lat, double* lon, double* value);
int grib_iterator_reset(grib_iterator* i);
int grib_iterator_init(grib_iterator* i, grib_handle* h, grib_arguments* args);
int grib_iterator_delete(grib_iterator* i);
//int grib_iterator_next(grib_iterator* i, double* lat, double* lon, double* value);
//int grib_iterator_has_next(grib_iterator* i);
//int grib_iterator_previous(grib_iterator* i, double* lat, double* lon, double* value);
//int grib_iterator_reset(grib_iterator* i);
//int grib_iterator_init(grib_iterator* i, grib_handle* h, grib_arguments* args);
//int grib_iterator_delete(grib_iterator* i);
/* grib_iterator_class.cc */
grib_iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsigned long flags, int* error);

36
src/grib_api_internal.h
View File

@ -248,7 +248,7 @@ typedef struct grib_codetable grib_codetable;
typedef struct grib_smart_table grib_smart_table;
class grib_accessor;
typedef struct grib_iterator_class grib_iterator_class;
class grib_iterator;
typedef struct grib_nearest_class grib_nearest_class;
typedef struct grib_dumper grib_dumper;
typedef struct grib_dumper_class grib_dumper_class;
@ -267,14 +267,6 @@ typedef int (*nearest_find_proc)(grib_nearest* nearest, grib_handle* h,
double* distances, int* indexes, size_t* len);
typedef int (*nearest_destroy_proc)(grib_nearest* nearest);
typedef void (*iterator_init_class_proc)(grib_iterator_class*);
typedef int (*iterator_init_proc)(grib_iterator* i, grib_handle*, grib_arguments*);
typedef int (*iterator_next_proc)(grib_iterator* i, double* lat, double* lon, double* val);
typedef int (*iterator_previous_proc)(grib_iterator* i, double* lat, double* lon, double* val);
typedef int (*iterator_reset_proc)(grib_iterator* i);
typedef int (*iterator_destroy_proc)(grib_iterator* i);
typedef long (*iterator_has_next_proc)(grib_iterator* i);
typedef int (*grib_pack_proc)(grib_handle* h, const double* in, size_t inlen, void* out, size_t* outlen);
typedef int (*grib_unpack_proc)(grib_handle* h, const void* in, size_t inlen, double* out, size_t* outlen);
@ -459,20 +451,6 @@ struct grib_section
size_t padding;
};
struct grib_iterator_class
{
grib_iterator_class** super;
const char* name;
size_t size;
int inited;
iterator_init_class_proc init_class;
iterator_init_proc init;
iterator_destroy_proc destroy;
iterator_next_proc next;
iterator_previous_proc previous;
iterator_reset_proc reset;
iterator_has_next_proc has_next;
};
struct grib_nearest_class
{
@ -529,17 +507,6 @@ struct grib_dumper_class
dumper_footer_proc footer;
};
struct grib_iterator
{
grib_arguments* args; /** args of iterator */
grib_handle* h;
long e; /** current element */
size_t nv; /** number of values */
double* data; /** data values */
grib_iterator_class* cclass;
unsigned long flags;
};
struct grib_nearest
{
grib_arguments* args; /** args of iterator */
@ -1287,6 +1254,7 @@ typedef struct j2k_encode_helper
}
#include "accessor/grib_accessor.h"
#include "accessor/grib_accessors_list.h"
#include "iterator/grib_iterator.h"
#endif
#endif

168
src/grib_iterator.cc
View File

@ -1,168 +0,0 @@
/*
* (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.
*/
/***************************************************************************
* Jean Baptiste Filippi - 01.11.2005 *
***************************************************************************/
#include "grib_api_internal.h"
#if GRIB_PTHREADS
static pthread_once_t once = PTHREAD_ONCE_INIT;
static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
static void init_mutex()
{
pthread_mutexattr_t attr;
pthread_mutexattr_init(&attr);
pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
pthread_mutex_init(&mutex, &attr);
pthread_mutexattr_destroy(&attr);
}
#elif GRIB_OMP_THREADS
static int once = 0;
static omp_nest_lock_t mutex;
static void init_mutex()
{
GRIB_OMP_CRITICAL(lock_iterator_c)
{
if (once == 0) {
omp_init_nest_lock(&mutex);
once = 1;
}
}
}
#endif
int grib_get_data(const grib_handle* h, double* lats, double* lons, double* values)
{
int err = 0;
grib_iterator* iter = NULL;
double *lat, *lon, *val;
iter = grib_iterator_new(h, 0, &err);
if (!iter || err != GRIB_SUCCESS)
return err;
lat = lats;
lon = lons;
val = values;
while (grib_iterator_next(iter, lat++, lon++, val++)) {}
grib_iterator_delete(iter);
return err;
}
int grib_iterator_next(grib_iterator* i, double* lat, double* lon, double* value)
{
grib_iterator_class* c = i->cclass;
while (c) {
grib_iterator_class* s = c->super ? *(c->super) : NULL;
if (c->next)
return c->next(i, lat, lon, value);
c = s;
}
grib_context_log(i->h->context, GRIB_LOG_FATAL, "%s: No next() function in iterator '%s'", __func__, i->cclass->name);
return 0;
}
int grib_iterator_has_next(grib_iterator* i)
{
grib_iterator_class* c = i->cclass;
while (c) {
grib_iterator_class* s = c->super ? *(c->super) : NULL;
if (c->has_next)
return c->has_next(i);
c = s;
}
grib_context_log(i->h->context, GRIB_LOG_FATAL, "%s: No has_next() function in iterator '%s'", __func__, i->cclass->name);
return 0;
}
int grib_iterator_previous(grib_iterator* i, double* lat, double* lon, double* value)
{
grib_iterator_class* c = i->cclass;
while (c) {
grib_iterator_class* s = c->super ? *(c->super) : NULL;
if (c->previous)
return c->previous(i, lat, lon, value);
c = s;
}
grib_context_log(i->h->context, GRIB_LOG_FATAL, "%s: No previous() function in iterator '%s'", __func__, i->cclass->name);
return 0;
}
int grib_iterator_reset(grib_iterator* i)
{
grib_iterator_class* c = i->cclass;
while (c) {
grib_iterator_class* s = c->super ? *(c->super) : NULL;
if (c->reset)
return c->reset(i);
c = s;
}
grib_context_log(i->h->context, GRIB_LOG_FATAL, "%s: No reset() function in iterator '%s'", __func__, i->cclass->name);
return 0;
}
/* For this one, ALL init are called */
static int init_iterator(grib_iterator_class* c, grib_iterator* i, grib_handle* h, grib_arguments* args)
{
if (c) {
int ret = GRIB_SUCCESS;
grib_iterator_class* s = c->super ? *(c->super) : NULL;
if (!c->inited) {
if (c->init_class)
c->init_class(c);
c->inited = 1;
}
if (s)
ret = init_iterator(s, i, h, args);
if (ret != GRIB_SUCCESS)
return ret;
if (c->init)
return c->init(i, h, args);
}
return GRIB_INTERNAL_ERROR;
}
int grib_iterator_init(grib_iterator* i, grib_handle* h, grib_arguments* args)
{
int r = 0;
GRIB_MUTEX_INIT_ONCE(&once, &init_mutex);
GRIB_MUTEX_LOCK(&mutex);
r = init_iterator(i->cclass, i, h, args);
GRIB_MUTEX_UNLOCK(&mutex);
return r;
}
/* For this one, ALL destroy are called */
int grib_iterator_delete(grib_iterator* i)
{
if (i) {
grib_iterator_class* c = i->cclass;
while (c) {
grib_iterator_class* s = c->super ? *(c->super) : NULL;
if (c->destroy)
c->destroy(i);
c = s;
}
/* This should go in a top class */
grib_context_free(i->h->context, i);
}
else {
return GRIB_INVALID_ARGUMENT;
}
return 0;
}

11
src/grib_iterator_class.cc
View File

@ -21,7 +21,7 @@
struct table_entry
{
const char* type;
grib_iterator_class** cclass;
grib_iterator** cclass;
};
static const struct table_entry table[] = {
@ -38,11 +38,10 @@ grib_iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsig
num_table_entries = sizeof(table) / sizeof(table[0]);
for (i = 0; i < num_table_entries; i++) {
if (strcmp(type, table[i].type) == 0) {
grib_iterator_class* c = *(table[i].cclass);
grib_iterator* it = (grib_iterator*)grib_context_malloc_clear(h->context, c->size);
it->cclass = c;
it->flags = flags;
*error = grib_iterator_init(it, h, args);
grib_iterator* c = *(table[i].cclass);
grib_iterator* it = (grib_iterator*)grib_context_malloc_clear(h->context, c->size_);
it->flags_ = flags;
*error = it->init(h, args);
if (*error == GRIB_SUCCESS)
return it;
grib_context_log(h->context, GRIB_LOG_ERROR, "Geoiterator factory: Error instantiating iterator %s (%s)",

24
src/grib_iterator_class.h
View File

@ -1,13 +1,13 @@
/* This file is automatically generated by ./make_class.pl, do not edit */
extern grib_iterator_class* grib_iterator_class_gaussian;
extern grib_iterator_class* grib_iterator_class_gaussian_reduced;
extern grib_iterator_class* grib_iterator_class_gen;
extern grib_iterator_class* grib_iterator_class_healpix;
extern grib_iterator_class* grib_iterator_class_lambert_azimuthal_equal_area;
extern grib_iterator_class* grib_iterator_class_lambert_conformal;
extern grib_iterator_class* grib_iterator_class_latlon;
extern grib_iterator_class* grib_iterator_class_latlon_reduced;
extern grib_iterator_class* grib_iterator_class_mercator;
extern grib_iterator_class* grib_iterator_class_polar_stereographic;
extern grib_iterator_class* grib_iterator_class_regular;
extern grib_iterator_class* grib_iterator_class_space_view;
extern grib_iterator* grib_iterator_gaussian;
extern grib_iterator* grib_iterator_gaussian_reduced;
extern grib_iterator* grib_iterator_gen;
extern grib_iterator* grib_iterator_healpix;
extern grib_iterator* grib_iterator_lambert_azimuthal_equal_area;
extern grib_iterator* grib_iterator_lambert_conformal;
extern grib_iterator* grib_iterator_latlon;
extern grib_iterator* grib_iterator_latlon_reduced;
extern grib_iterator* grib_iterator_mercator;
extern grib_iterator* grib_iterator_polar_stereographic;
extern grib_iterator* grib_iterator_regular;
extern grib_iterator* grib_iterator_space_view;

224
src/grib_iterator_class_latlon.cc
View File

@ -1,224 +0,0 @@
/*
* (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;
}

205
src/grib_iterator_class_latlon_reduced.cc
View File

@ -1,205 +0,0 @@
/*
* (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 <cmath>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_gen
IMPLEMENTS = destroy
IMPLEMENTS = init;next
MEMBERS = double *las
MEMBERS = double *los
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);
static int destroy (grib_iterator* i);
typedef struct grib_iterator_latlon_reduced{
grib_iterator it;
/* Members defined in gen */
int carg;
const char* missingValue;
/* Members defined in latlon_reduced */
double *las;
double *los;
} grib_iterator_latlon_reduced;
extern grib_iterator_class* grib_iterator_class_gen;
static grib_iterator_class _grib_iterator_class_latlon_reduced = {
&grib_iterator_class_gen, /* super */
"latlon_reduced", /* name */
sizeof(grib_iterator_latlon_reduced),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
&next, /* Next Value */
0, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_latlon_reduced = &_grib_iterator_class_latlon_reduced;
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_iterator_latlon_reduced* self = (grib_iterator_latlon_reduced*)iter;
if ((long)iter->e >= (long)(iter->nv - 1))
return 0;
iter->e++;
*lat = self->las[iter->e];
*lon = self->los[iter->e];
if (val && iter->data) {
*val = iter->data[iter->e];
}
return 1;
}
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
{
grib_iterator_latlon_reduced* self = (grib_iterator_latlon_reduced*)iter;
int ret = GRIB_SUCCESS;
double laf;
double lal;
long nlats;
double lof, tlof;
double lol, dimin;
long* pl;
size_t plsize = 0;
long k, j, ii;
long nlons, plmax;
double jdirinc = 0;
double idirinc = 0;
double dlon = 0;
int islocal = 0;
long nlons2 = 0; /* adjusted num of longitudes */
const char* latofirst = grib_arguments_get_name(h, args, self->carg++);
const char* longoffirst = grib_arguments_get_name(h, args, self->carg++);
const char* latoflast = grib_arguments_get_name(h, args, self->carg++);
const char* longoflast = grib_arguments_get_name(h, args, self->carg++);
const char* nlats_name = grib_arguments_get_name(h, args, self->carg++);
const char* jdirec = grib_arguments_get_name(h, args, self->carg++);
const char* plac = grib_arguments_get_name(h, args, self->carg++);
if ((ret = grib_get_double_internal(h, latofirst, &laf)))
return ret;
if ((ret = grib_get_double_internal(h, longoffirst, &lof)))
return ret;
if ((ret = grib_get_double_internal(h, latoflast, &lal)))
return ret;
if ((ret = grib_get_double_internal(h, longoflast, &lol)))
return ret;
if ((ret = grib_get_long_internal(h, nlats_name, &nlats)))
return ret;
if ((ret = grib_get_double_internal(h, jdirec, &jdirinc)))
return ret;
plsize = nlats;
pl = (long*)grib_context_malloc(h->context, plsize * sizeof(long));
grib_get_long_array_internal(h, plac, pl, &plsize);
self->las = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
self->los = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
plmax = pl[0];
for (j = 0; j < nlats; j++)
if (plmax < pl[j])
plmax = pl[j];
dimin = 360.0 / plmax;
if (360 - fabs(lol - lof) < 2 * dimin) {
dlon = 360;
islocal = 0;
}
else if (lol < lof) {
/* handle something like 150 to -120 to treat as 150 to 240 */
/* so that dlon is 90 (not -270) */
dlon = lol + 360.0 - lof;
islocal = 1;
}
else {
dlon = lol - lof;
islocal = 1;
}
if (laf > lal)
jdirinc = -jdirinc;
k = 0;
for (j = 0; j < nlats; j++) {
nlons = pl[j];
tlof = lof;
nlons2 = nlons - islocal;
/*Sometimes there are no points on a latitude! Protect against div by zero*/
if (nlons2 < 1)
nlons2 = 1;
idirinc = dlon / nlons2;
for (ii = 0; ii < nlons; ii++) {
self->las[k] = laf;
self->los[k] = tlof;
tlof += idirinc;
k++;
}
laf += jdirinc;
}
iter->e = -1;
grib_context_free(h->context, pl);
return ret;
}
static int destroy(grib_iterator* iter)
{
grib_iterator_latlon_reduced* self = (grib_iterator_latlon_reduced*)iter;
const grib_context* c = iter->h->context;
grib_context_free(c, self->las);
grib_context_free(c, self->los);
return GRIB_SUCCESS;
}

246
src/grib_iterator_class_regular.cc
View File

@ -1,246 +0,0 @@
/*
* (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 <cmath>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_gen
IMPLEMENTS = previous;next
IMPLEMENTS = init;destroy
MEMBERS = double *las
MEMBERS = double *los
MEMBERS = long Ni
MEMBERS = long Nj
MEMBERS = long iScansNegatively
MEMBERS = long isRotated
MEMBERS = double angleOfRotation
MEMBERS = double southPoleLat
MEMBERS = double southPoleLon
MEMBERS = long jPointsAreConsecutive
MEMBERS = long disableUnrotate
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);
static int previous (grib_iterator* ei, double *lat, double *lon, double *val);
static int destroy (grib_iterator* i);
typedef struct grib_iterator_regular{
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;
} grib_iterator_regular;
extern grib_iterator_class* grib_iterator_class_gen;
static grib_iterator_class _grib_iterator_class_regular = {
&grib_iterator_class_gen, /* super */
"regular", /* name */
sizeof(grib_iterator_regular),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
&next, /* Next Value */
&previous, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_regular = &_grib_iterator_class_regular;
static void init_class(grib_iterator_class* c)
{
c->reset = (*(c->super))->reset;
c->has_next = (*(c->super))->has_next;
}
/* END_CLASS_IMP */
#define ITER "Regular grid Geoiterator"
static int next(grib_iterator* iter, double* lat, double* lon, double* val)
{
grib_iterator_regular* self = (grib_iterator_regular*)iter;
if ((long)iter->e >= (long)(iter->nv - 1))
return 0;
iter->e++;
*lat = self->las[(long)floor(iter->e / self->Ni)];
*lon = self->los[(long)iter->e % self->Ni];
if (val && iter->data) {
*val = iter->data[iter->e];
}
return 1;
}
static int previous(grib_iterator* iter, double* lat, double* lon, double* val)
{
grib_iterator_regular* self = (grib_iterator_regular*)iter;
if (iter->e < 0)
return 0;
*lat = self->las[(long)floor(iter->e / self->Ni)];
*lon = self->los[iter->e % self->Ni];
if (val && iter->data) {
*val = iter->data[iter->e];
}
iter->e--;
return 1;
}
static int destroy(grib_iterator* iter)
{
grib_iterator_regular* self = (grib_iterator_regular*)iter;
const grib_context* c = iter->h->context;
grib_context_free(c, self->las);
grib_context_free(c, self->los);
return GRIB_SUCCESS;
}
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
{
grib_iterator_regular* self = (grib_iterator_regular*)iter;
int ret = GRIB_SUCCESS;
long Ni; /* Number of points along a parallel = Nx */
long Nj; /* Number of points along a meridian = Ny */
double idir, idir_coded, lon1, lon2;
long loi;
const char* s_lon1 = grib_arguments_get_name(h, args, self->carg++);
const char* s_idir = grib_arguments_get_name(h, args, self->carg++);
const char* s_Ni = grib_arguments_get_name(h, args, self->carg++);
const char* s_Nj = grib_arguments_get_name(h, args, self->carg++);
const char* s_iScansNeg = grib_arguments_get_name(h, args, self->carg++);
if ((ret = grib_get_double_internal(h, s_lon1, &lon1)))
return ret;
if ((ret = grib_get_double_internal(h, "longitudeOfLastGridPointInDegrees", &lon2)))
return ret;
if ((ret = grib_get_double_internal(h, s_idir, &idir))) // can be GRIB_MISSING_DOUBLE
return ret;
idir_coded = idir;
if ((ret = grib_get_long_internal(h, s_Ni, &Ni)))
return ret;
if (grib_is_missing(h, s_Ni, &ret) && ret == GRIB_SUCCESS) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Key %s cannot be 'missing' for a regular grid!", ITER, s_Ni);
return GRIB_WRONG_GRID;
}
if ((ret = grib_get_long_internal(h, s_Nj, &Nj)))
return ret;
if (grib_is_missing(h, s_Nj, &ret) && ret == GRIB_SUCCESS) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Key %s cannot be 'missing' for a regular grid!", ITER, s_Nj);
return GRIB_WRONG_GRID;
}
if (Ni*Nj != iter->nv) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Ni*Nj!=numberOfDataPoints (%ld*%ld!=%zu)", ITER, Ni, Nj, iter->nv);
return GRIB_WRONG_GRID;
}
if ((ret = grib_get_long_internal(h, s_iScansNeg, &self->iScansNegatively)))
return ret;
/* GRIB-801: Careful of case with a single point! Ni==1 */
if (Ni > 1) {
/* Note: If first and last longitudes are equal I assume you wanna go round the globe */
if (self->iScansNegatively) {
if (lon1 > lon2) {
idir = (lon1 - lon2) / (Ni - 1);
}
else {
idir = (lon1 + 360.0 - lon2) / (Ni - 1);
}
}
else {
if (lon2 > lon1) {
idir = (lon2 - lon1) / (Ni - 1);
}
else {
idir = (lon2 + 360.0 - lon1) / (Ni - 1);
}
}
}
if (self->iScansNegatively) {
idir = -idir;
}
else {
if (lon1 + (Ni - 2) * idir > 360)
lon1 -= 360;
/*See ECC-704, GRIB-396*/
/*else if ( (lon1+(Ni-1)*idir)-360 > epsilon ){
idir=360.0/(float)Ni;
}*/
}
self->Ni = Ni;
self->Nj = Nj;
self->las = (double*)grib_context_malloc(h->context, Nj * sizeof(double));
self->los = (double*)grib_context_malloc(h->context, Ni * sizeof(double));
if (idir != idir_coded) {
grib_context_log(h->context, GRIB_LOG_DEBUG, "%s: Using idir=%g (coded value=%g)", ITER, idir, idir_coded);
}
for (loi = 0; loi < Ni; loi++) {
self->los[loi] = lon1;
lon1 += idir;
}
// ECC-1406: Due to rounding, errors can accumulate.
// So we ensure the last longitude is longitudeOfLastGridPointInDegrees
// Also see ECC-1671, ECC-1708
if (lon2 > 0) {
lon2 = normalise_longitude_in_degrees(lon2);
}
self->los[Ni-1] = lon2;
return ret;
}

26
src/grib_iterator_factory.h
View File

@ -1,13 +1,13 @@
/* This file is automatically generated by ./make_class.pl, do not edit */
{ "gaussian", &grib_iterator_class_gaussian, },
{ "gaussian_reduced", &grib_iterator_class_gaussian_reduced, },
{ "gen", &grib_iterator_class_gen, },
{ "healpix", &grib_iterator_class_healpix, },
{ "lambert_azimuthal_equal_area", &grib_iterator_class_lambert_azimuthal_equal_area, },
{ "lambert_conformal", &grib_iterator_class_lambert_conformal, },
{ "latlon", &grib_iterator_class_latlon, },
{ "latlon_reduced", &grib_iterator_class_latlon_reduced, },
{ "mercator", &grib_iterator_class_mercator, },
{ "polar_stereographic", &grib_iterator_class_polar_stereographic, },
{ "regular", &grib_iterator_class_regular, },
{ "space_view", &grib_iterator_class_space_view, },
/* This file is automatically generated by ./make.pl, do not edit */
{ "gaussian", &grib_iterator_gaussian, },
{ "gaussian_reduced", &grib_iterator_gaussian_reduced, },
{ "gen", &grib_iterator_gen, },
{ "healpix", &grib_iterator_healpix, },
{ "lambert_azimuthal_equal_area", &grib_iterator_lambert_azimuthal_equal_area, },
{ "lambert_conformal", &grib_iterator_lambert_conformal, },
{ "latlon", &grib_iterator_latlon, },
{ "latlon_reduced", &grib_iterator_latlon_reduced, },
{ "mercator", &grib_iterator_mercator, },
{ "polar_stereographic", &grib_iterator_polar_stereographic, },
{ "regular", &grib_iterator_regular, },
{ "space_view", &grib_iterator_space_view, },

40
src/iterator/grib_iterator.h Normal file
View File

@ -0,0 +1,40 @@
/*
* (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.
*/
#pragma once
#include "grib_api_internal.h"
/*! Grib geoiterator, structure supporting a geographic iteration of values on a GRIB message.
\ingroup grib_iterator
*/
class grib_iterator
{
public:
virtual int init(grib_handle*, grib_arguments*) = 0;
virtual int next(double*, double*, double*) = 0;
virtual int previous(double*, double*, double*) = 0;
virtual int reset() = 0;
virtual int destroy() = 0;
virtual long has_next() = 0;
public:
grib_context* context_;
grib_arguments* args_; /** args of iterator */
grib_handle* h_;
long e_; /** current element */
size_t nv_; /** number of values */
double* data_; /** data values */
unsigned long flags_;
const char* class_name_;
size_t size_;
int inited_;
};

98
src/grib_iterator_class_gaussian.cc → src/iterator/grib_iterator_class_gaussian.cc
View File

@ -8,89 +8,17 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_iterator_class_gaussian.h"
#include <cmath>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_regular
IMPLEMENTS = init
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*);
typedef struct grib_iterator_gaussian{
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 gaussian */
} grib_iterator_gaussian;
extern grib_iterator_class* grib_iterator_class_regular;
static grib_iterator_class _grib_iterator_class_gaussian = {
&grib_iterator_class_regular, /* super */
"gaussian", /* name */
sizeof(grib_iterator_gaussian),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
0, /* destructor */
0, /* Next Value */
0, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_gaussian = &_grib_iterator_class_gaussian;
static void init_class(grib_iterator_class* c)
{
c->next = (*(c->super))->next;
c->previous = (*(c->super))->previous;
c->reset = (*(c->super))->reset;
c->has_next = (*(c->super))->has_next;
}
/* END_CLASS_IMP */
grib_iterator_gaussian_t _grib_iterator_gaussian{};
grib_iterator* grib_iterator_gaussian = &_grib_iterator_gaussian;
static void binary_search_gaussian_latitudes(const double xx[], const unsigned long n, double x, long* j);
static int init(grib_iterator* i, grib_handle* h, grib_arguments* args)
int grib_iterator_gaussian_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gaussian* self = (grib_iterator_gaussian*)i;
grib_iterator_regular_t::init(h, args);
double* lats;
double laf; /* latitude of first point in degrees */
@ -103,10 +31,10 @@ static int init(grib_iterator* i, grib_handle* h, grib_arguments* args)
long istart = 0;
int ret = GRIB_SUCCESS;
const char* latofirst = grib_arguments_get_name(h, args, self->carg++);
const char* latoflast = grib_arguments_get_name(h, args, self->carg++);
const char* numtrunc = grib_arguments_get_name(h, args, self->carg++);
const char* s_jScansPositively = grib_arguments_get_name(h, args, self->carg++);
const char* latofirst = grib_arguments_get_name(h, args, carg_++);
const char* latoflast = grib_arguments_get_name(h, args, carg_++);
const char* numtrunc = grib_arguments_get_name(h, args, carg_++);
const char* s_jScansPositively = grib_arguments_get_name(h, args, carg_++);
if ((ret = grib_get_double_internal(h, latofirst, &laf)))
return ret;
@ -144,15 +72,15 @@ static int init(grib_iterator* i, grib_handle* h, grib_arguments* args)
}
if (jScansPositively) {
for (lai = 0; lai < self->Nj; lai++) {
for (lai = 0; lai < Nj_; lai++) {
DEBUG_ASSERT(istart >= 0);
self->las[lai] = lats[istart--];
las_[lai] = lats[istart--];
if (istart<0) istart=size-1;
}
}
else {
for (lai = 0; lai < self->Nj; lai++) {
self->las[lai] = lats[istart++];
for (lai = 0; lai < Nj_; lai++) {
las_[lai] = lats[istart++];
if (istart > size - 1)
istart = 0;
}

21
src/iterator/grib_iterator_class_gaussian.h Normal file
View File

@ -0,0 +1,21 @@
/*
* (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_iterator_class_regular.h"
class grib_iterator_gaussian_t : public grib_iterator_regular_t
{
public:
grib_iterator_gaussian_t() : grib_iterator_regular_t()
{
class_name_ = "gaussian";
}
int init(grib_handle*, grib_arguments*);
};

212
src/grib_iterator_class_gaussian_reduced.cc → src/iterator/grib_iterator_class_gaussian_reduced.cc
View File

@ -8,111 +8,33 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_iterator_class_gaussian_reduced.h"
#include <cmath>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_gen
IMPLEMENTS = destroy
IMPLEMENTS = init;next
MEMBERS = double *las
MEMBERS = double *los
MEMBERS = long Nj
MEMBERS = long isRotated
MEMBERS = double angleOfRotation
MEMBERS = double southPoleLat
MEMBERS = double southPoleLon
MEMBERS = long disableUnrotate
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);
static int destroy (grib_iterator* i);
typedef struct grib_iterator_gaussian_reduced{
grib_iterator it;
/* Members defined in gen */
int carg;
const char* missingValue;
/* Members defined in gaussian_reduced */
double *las;
double *los;
long Nj;
long isRotated;
double angleOfRotation;
double southPoleLat;
double southPoleLon;
long disableUnrotate;
} grib_iterator_gaussian_reduced;
extern grib_iterator_class* grib_iterator_class_gen;
static grib_iterator_class _grib_iterator_class_gaussian_reduced = {
&grib_iterator_class_gen, /* super */
"gaussian_reduced", /* name */
sizeof(grib_iterator_gaussian_reduced),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
&next, /* Next Value */
0, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_gaussian_reduced = &_grib_iterator_class_gaussian_reduced;
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 */
grib_iterator_gaussian_reduced_t _grib_iterator_gaussian_reduced{};
grib_iterator* grib_iterator_gaussian_reduced = &_grib_iterator_gaussian_reduced;
#define ITER "Reduced Gaussian grid Geoiterator"
static int next(grib_iterator* iter, double* lat, double* lon, double* val)
int grib_iterator_gaussian_reduced_t::next(double* lat, double* lon, double* val)
{
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter;
double ret_lat=0, ret_lon=0;
if (iter->e >= (long)(iter->nv - 1))
if (e_ >= (long)(nv_ - 1))
return 0;
iter->e++;
e_++;
ret_lat = self->las[iter->e];
ret_lon = self->los[iter->e];
if (val && iter->data) {
*val = iter->data[iter->e];
ret_lat = las_[e_];
ret_lon = los_[e_];
if (val && data_) {
*val = data_[e_];
}
if (self->isRotated && !self->disableUnrotate) {
if (isRotated_ && !disableUnrotate_) {
double new_lat = 0, new_lon = 0;
unrotate(ret_lat, ret_lon,
self->angleOfRotation, self->southPoleLat, self->southPoleLon,
angleOfRotation_, southPoleLat_, southPoleLon_,
&new_lat, &new_lon);
ret_lat = new_lat;
ret_lon = new_lon;
@ -172,7 +94,7 @@ static void binary_search(const double xx[], const unsigned long n, double x, lo
}
/* Use legacy way to compute the iterator latitude/longitude values */
static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_handle* h,
int grib_iterator_gaussian_reduced_t::iterate_reduced_gaussian_subarea_legacy(grib_handle* h,
double lat_first, double lon_first,
double lat_last, double lon_last,
double* lats, long* pl, size_t plsize)
@ -183,7 +105,6 @@ static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_han
long row_count = 0;
double d = 0;
long ilon_first, ilon_last, i;
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter;
/*get_reduced_row_proc get_reduced_row = &grib_get_reduced_row;*/
get_reduced_row_proc get_reduced_row = &grib_get_reduced_row_legacy; /* legacy algorithm */
@ -198,7 +119,7 @@ static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_han
l++;
}
iter->e = 0;
e_ = 0;
for (j = 0; j < plsize; j++) {
long k = 0;
row_count = 0;
@ -207,16 +128,16 @@ static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_han
if (ilon_first > ilon_last)
ilon_first -= pl[j];
for (i = ilon_first; i <= ilon_last; i++) {
if (iter->e >= iter->nv) {
if (e_ >= nv_) {
size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last);
grib_context_log(h->context, GRIB_LOG_ERROR,
"%s (sub-area legacy). Num points=%zu, size(values)=%zu", ITER, np, iter->nv);
"%s (sub-area legacy). Num points=%zu, size(values)=%zu", ITER, np, nv_);
return GRIB_WRONG_GRID;
}
self->los[iter->e] = ((i)*360.0) / pl[j];
self->las[iter->e] = lats[j + l];
iter->e++;
los_[e_] = ((i)*360.0) / pl[j];
las_[e_] = lats[j + l];
e_++;
k++;
if (k >= row_count) {
/* Ensure we exit the loop and only process 'row_count' points */
@ -235,7 +156,7 @@ static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_han
// err = iterate_reduced_gaussian_subarea_algorithm2(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize);
// }
// return err;
static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
int grib_iterator_gaussian_reduced_t::iterate_reduced_gaussian_subarea(grib_handle* h,
double lat_first, double lon_first,
double lat_last, double lon_last,
double* lats, long* pl, size_t plsize, size_t numlats)
@ -245,7 +166,6 @@ static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
size_t j = 0;
long row_count = 0, i = 0;
double olon_first, olon_last;
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter;
get_reduced_row_proc get_reduced_row = &grib_get_reduced_row;
if (h->context->debug) {
@ -265,33 +185,33 @@ static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
// }
// }
iter->e = 0;
e_ = 0;
for (j = 0; j < plsize; j++) {
const double delta = 360.0 / pl[j];
row_count = 0;
grib_get_reduced_row_p(pl[j], lon_first, lon_last, &row_count, &olon_first, &olon_last);
for (i = 0; i < row_count; ++i) {
double lon2 = olon_first + i * delta;
if (iter->e >= iter->nv) {
if (e_ >= nv_) {
/* Only print error message on the second pass */
size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last);
grib_context_log(h->context, GRIB_LOG_ERROR,
"%s (sub-area). Num points=%zu, size(values)=%zu", ITER, np, iter->nv);
"%s (sub-area). Num points=%zu, size(values)=%zu", ITER, np, nv_);
return GRIB_WRONG_GRID;
}
self->los[iter->e] = lon2;
los_[e_] = lon2;
DEBUG_ASSERT(j + l < numlats);
self->las[iter->e] = lats[j + l];
iter->e++;
las_[e_] = lats[j + l];
e_++;
}
}
if (iter->e != iter->nv) {
if (e_ != nv_) {
/* Fewer counted points in the sub-area than the number of data values */
const size_t legacy_count = count_subarea_points(h, grib_get_reduced_row_legacy, pl, plsize, lon_first, lon_last);
if (iter->nv == legacy_count) {
if (nv_ == legacy_count) {
/* Legacy (produced by PRODGEN/LIBEMOS) */
return iterate_reduced_gaussian_subarea_legacy(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize);
return iterate_reduced_gaussian_subarea_legacy(h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize);
}
else {
/* TODO: A gap exists! Not all values can be mapped. Inconsistent grid or error in calculating num. points! */
@ -300,8 +220,10 @@ static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
return err;
}
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
int grib_iterator_gaussian_reduced_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gen_t::init(h, args);
int ret = GRIB_SUCCESS, j, is_global = 0;
double lat_first = 0, lon_first = 0, lat_last = 0, lon_last = 0;
double angular_precision = 1.0 / 1000000.0;
@ -314,28 +236,27 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
long row_count = 0;
long angleSubdivisions = 0;
const grib_context* c = h->context;
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter;
const char* slat_first = grib_arguments_get_name(h, args, self->carg++);
const char* slon_first = grib_arguments_get_name(h, args, self->carg++);
const char* slat_last = grib_arguments_get_name(h, args, self->carg++);
const char* slon_last = grib_arguments_get_name(h, args, self->carg++);
const char* sorder = grib_arguments_get_name(h, args, self->carg++);
const char* spl = grib_arguments_get_name(h, args, self->carg++);
const char* snj = grib_arguments_get_name(h, args, self->carg++);
const char* slat_first = grib_arguments_get_name(h, args, carg_++);
const char* slon_first = grib_arguments_get_name(h, args, carg_++);
const char* slat_last = grib_arguments_get_name(h, args, carg_++);
const char* slon_last = grib_arguments_get_name(h, args, carg_++);
const char* sorder = grib_arguments_get_name(h, args, carg_++);
const char* spl = grib_arguments_get_name(h, args, carg_++);
const char* snj = grib_arguments_get_name(h, args, carg_++);
self->angleOfRotation = 0;
self->isRotated = 0;
self->southPoleLat = 0;
self->southPoleLon = 0;
self->disableUnrotate = 0; /* unrotate enabled by default */
angleOfRotation_ = 0;
isRotated_ = 0;
southPoleLat_ = 0;
southPoleLon_ = 0;
disableUnrotate_ = 0; /* unrotate enabled by default */
ret = grib_get_long(h, "isRotatedGrid", &self->isRotated);
if (ret == GRIB_SUCCESS && self->isRotated) {
if ((ret = grib_get_double_internal(h, "angleOfRotation", &self->angleOfRotation)))
ret = grib_get_long(h, "isRotatedGrid", &isRotated_);
if (ret == GRIB_SUCCESS && isRotated_) {
if ((ret = grib_get_double_internal(h, "angleOfRotation", &angleOfRotation_)))
return ret;
if ((ret = grib_get_double_internal(h, "latitudeOfSouthernPoleInDegrees", &self->southPoleLat)))
if ((ret = grib_get_double_internal(h, "latitudeOfSouthernPoleInDegrees", &southPoleLat_)))
return ret;
if ((ret = grib_get_double_internal(h, "longitudeOfSouthernPoleInDegrees", &self->southPoleLon)))
if ((ret = grib_get_double_internal(h, "longitudeOfSouthernPoleInDegrees", &southPoleLon_)))
return ret;
}
@ -379,11 +300,11 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
grib_get_long_array_internal(h, spl, pl, &plsize);
self->las = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
if (!self->las)
las_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (!las_)
return GRIB_OUT_OF_MEMORY;
self->los = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
if (!self->los)
los_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (!los_)
return GRIB_OUT_OF_MEMORY;
while (lon_last < 0)
@ -402,11 +323,11 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
is_global = is_gaussian_global(lat_first, lat_last, lon_first, lon_last, max_pl, lats, angular_precision);
if (!is_global) {
/*sub area*/
ret = iterate_reduced_gaussian_subarea(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize, numlats);
ret = iterate_reduced_gaussian_subarea(h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize, numlats);
}
else {
/*global*/
iter->e = 0;
e_ = 0;
if (h->context->debug) {
const size_t np = count_global_points(pl, plsize);
fprintf(stderr, "ECCODES DEBUG grib_iterator_class_gaussian_reduced: global num points=%zu\n", np);
@ -415,37 +336,38 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
for (j = 0; j < plsize; j++) {
row_count = pl[j];
for (i = 0; i < row_count; i++) {
if (iter->e >= iter->nv) {
if (e_ >= nv_) {
/*grib_context_log(h->context,GRIB_LOG_ERROR, "Failed to initialise reduced Gaussian iterator (global)");*/
/*return GRIB_WRONG_GRID;*/
/*Try now as NON-global*/
ret = iterate_reduced_gaussian_subarea(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize, numlats);
ret = iterate_reduced_gaussian_subarea(h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize, numlats);
if (ret != GRIB_SUCCESS)
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Failed to initialise iterator (global)", ITER);
goto finalise;
}
self->los[iter->e] = (i * 360.0) / row_count;
self->las[iter->e] = lats[j];
iter->e++;
los_[e_] = (i * 360.0) / row_count;
las_[e_] = lats[j];
e_++;
}
}
}
finalise:
iter->e = -1;
e_ = -1;
grib_context_free(h->context, lats);
grib_context_free(h->context, pl);
return ret;
}
static int destroy(grib_iterator* iter)
int grib_iterator_gaussian_reduced_t::destroy()
{
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter;
const grib_context* c = iter->h->context;
const grib_context* c = h_->context;
grib_context_free(c, self->las);
grib_context_free(c, self->los);
grib_context_free(c, las_);
grib_context_free(c, los_);
grib_iterator_gen_t::destroy();
return GRIB_SUCCESS;
}

46
src/iterator/grib_iterator_class_gaussian_reduced.h Normal file
View File

@ -0,0 +1,46 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_gen.h"
class grib_iterator_gaussian_reduced_t : public grib_iterator_gen_t
{
public:
grib_iterator_gaussian_reduced_t() : grib_iterator_gen_t()
{
class_name_ = "gaussian_reduced";
}
int init(grib_handle*,grib_arguments*) override;
int next(double *, double *, double *) override;
int destroy() override;
private:
double *las_;
double *los_;
long Nj_;
long isRotated_;
double angleOfRotation_;
double southPoleLat_;
double southPoleLon_;
long disableUnrotate_;
int iterate_reduced_gaussian_subarea_legacy(grib_handle*,
double, double,
double, double,
double*, long*, size_t);
int iterate_reduced_gaussian_subarea(grib_handle*,
double, double,
double, double,
double*, long*, size_t, size_t);
};

129
src/grib_iterator_class_gen.cc → src/iterator/grib_iterator_class_gen.cc
View File

@ -8,71 +8,11 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_iterator_class_gen.h"
/*
This is used by make_class.pl
//grib_iterator_gen_t _grib_iterator_class_gen{};
//grib_iterator* grib_iterator_class_gen = &_grib_iterator_class_gen;
START_CLASS_DEF
CLASS = iterator
IMPLEMENTS = destroy
IMPLEMENTS = has_next
IMPLEMENTS = init
IMPLEMENTS = reset
MEMBERS = int carg
MEMBERS = const char* missingValue;
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 destroy (grib_iterator* i);
static int reset (grib_iterator* i);
static long has_next (grib_iterator* i);
typedef struct grib_iterator_gen{
grib_iterator it;
/* Members defined in gen */
int carg;
const char* missingValue;
} grib_iterator_gen;
static grib_iterator_class _grib_iterator_class_gen = {
0, /* super */
"gen", /* name */
sizeof(grib_iterator_gen),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
0, /* Next Value */
0, /* Previous Value */
&reset, /* Reset the counter */
&has_next, /* has next values */
};
grib_iterator_class* grib_iterator_class_gen = &_grib_iterator_class_gen;
static void init_class(grib_iterator_class* c)
{
}
/* END_CLASS_IMP */
/*
* Return pointer to data at (i,j) (Fortran convention)
@ -170,23 +110,24 @@ int transform_iterator_data(grib_context* context, double* data,
return GRIB_SUCCESS;
}
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
int grib_iterator_gen_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gen* self = (grib_iterator_gen*)iter;
grib_iterator::init(h, args);
size_t dli = 0;
int err = GRIB_SUCCESS;
const char* s_rawData = NULL;
const char* s_numPoints = NULL;
long numberOfPoints = 0;
self->carg = 1;
carg_ = 1;
s_numPoints = grib_arguments_get_name(h, args, self->carg++);
self->missingValue = grib_arguments_get_name(h, args, self->carg++);
s_rawData = grib_arguments_get_name(h, args, self->carg++);
s_numPoints = grib_arguments_get_name(h, args, carg_++);
missingValue_ = grib_arguments_get_name(h, args, carg_++);
s_rawData = grib_arguments_get_name(h, args, carg_++);
iter->data = NULL;
iter->h = h; /* We may not need to keep them */
iter->args = args;
data_ = NULL;
h = h; /* We may not need to keep them */
args = args;
if ((err = grib_get_size(h, s_rawData, &dli)) != GRIB_SUCCESS)
return err;
@ -195,9 +136,9 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
// See ECC-1792. If we don't want to decode the Data Section, we should not do a check
// to see if it is consistent with the Grid Section
if (iter->flags & GRIB_GEOITERATOR_NO_VALUES) {
if (flags_ & GRIB_GEOITERATOR_NO_VALUES) {
// Iterator's number of values taken from the Grid Section
iter->nv = numberOfPoints;
nv_ = numberOfPoints;
} else {
// Check for consistency between the Grid and Data Sections
if (numberOfPoints != dli) {
@ -205,47 +146,59 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
s_numPoints, s_rawData, numberOfPoints, dli);
return GRIB_WRONG_GRID;
}
iter->nv = dli;
nv_ = dli;
}
if (iter->nv == 0) {
if (nv_ == 0) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Geoiterator: size(%s) is %ld", s_rawData, dli);
return GRIB_WRONG_GRID;
}
if ((iter->flags & GRIB_GEOITERATOR_NO_VALUES) == 0) {
if ((flags_ & GRIB_GEOITERATOR_NO_VALUES) == 0) {
// ECC-1525
// When the NO_VALUES flag is unset, decode the values and store them in the iterator.
// By default (and legacy) flags==0, so we decode
iter->data = (double*)grib_context_malloc(h->context, (iter->nv) * sizeof(double));
data_ = (double*)grib_context_malloc(h->context, (nv_) * sizeof(double));
if ((err = grib_get_double_array_internal(h, s_rawData, iter->data, &(iter->nv)))) {
if ((err = grib_get_double_array_internal(h, s_rawData, data_, &(nv_)))) {
return err;
}
}
iter->e = -1;
e_ = -1;
return err;
}
static int reset(grib_iterator* iter)
int grib_iterator_gen_t::reset()
{
iter->e = -1;
e_ = -1;
return 0;
}
static int destroy(grib_iterator* iter)
int grib_iterator_gen_t::destroy()
{
const grib_context* c = iter->h->context;
grib_context_free(c, iter->data);
const grib_context* c = h_->context;
grib_context_free(c, data_);
grib_iterator::destroy();
return GRIB_SUCCESS;
}
static long has_next(grib_iterator* iter)
long grib_iterator_gen_t::has_next()
{
if (iter->flags == 0 && iter->data == NULL)
if (flags_ == 0 && data_ == NULL)
return 0;
if (iter->e >= (long)(iter->nv - 1))
if (e_ >= (long)(nv_ - 1))
return 0;
return 1;
}
int grib_iterator_gen_t::previous(double*, double*, double*) {
return GRIB_NOT_IMPLEMENTED;
}
int grib_iterator_gen_t::next(double*, double*, double*) {
return GRIB_NOT_IMPLEMENTED;
}

33
src/iterator/grib_iterator_class_gen.h Normal file
View File

@ -0,0 +1,33 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator.h"
class grib_iterator_gen_t : public grib_iterator
{
public:
grib_iterator_gen_t() : grib_iterator()
{
class_name_ = "abstract_long_vector";
}
int init(grib_handle*,grib_arguments*) override;
int next(double*, double*, double*) override;
int previous(double*, double*, double*) override;
int reset() override;
int destroy() override;
long has_next() override;
public:
int carg_;
const char* missingValue_;
};

149
src/grib_iterator_class_healpix.cc → src/iterator/grib_iterator_class_healpix.cc
View File

@ -8,7 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_iterator_class_healpix.h"
#include <algorithm>
#include <bitset>
@ -17,77 +17,8 @@
#include <tuple>
#include <vector>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_gen
IMPLEMENTS = destroy
IMPLEMENTS = init;next
MEMBERS = double *lats
MEMBERS = double *lons
MEMBERS = long Nsides
MEMBERS = bool nested
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);
static int destroy (grib_iterator* i);
typedef struct grib_iterator_healpix{
grib_iterator it;
/* Members defined in gen */
int carg;
const char* missingValue;
/* Members defined in healpix */
double *lats;
double *lons;
long Nsides;
bool nested;
} grib_iterator_healpix;
extern grib_iterator_class* grib_iterator_class_gen;
static grib_iterator_class _grib_iterator_class_healpix = {
&grib_iterator_class_gen, /* super */
"healpix", /* name */
sizeof(grib_iterator_healpix),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
&next, /* Next Value */
0, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_healpix = &_grib_iterator_class_healpix;
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 */
grib_iterator_healpix_t _grib_iterator_healpix;
grib_iterator* grib_iterator_healpix = &_grib_iterator_healpix;
#define ITER "HEALPix Geoiterator"
constexpr double RAD2DEG = 57.29577951308232087684; // 180 over pi
@ -211,7 +142,9 @@ static std::vector<double> HEALPix_longitudes(size_t N, size_t i)
return longitudes;
}
static int iterate_healpix(grib_iterator_healpix* self, long N)
} // anonymous namespace
int grib_iterator_healpix_t::iterate_healpix(long N)
{
size_t Ny = 4 * static_cast<size_t>(N) - 1;
auto Nd = static_cast<double>(N);
@ -241,7 +174,7 @@ static int iterate_healpix(grib_iterator_healpix* self, long N)
// Equator
latitudes[2 * N - 1] = 0.0;
if (self->nested) {
if (nested_) {
if (!is_power_of_2(N)) {
grib_context* c = grib_context_get_default();
grib_context_log(c, GRIB_LOG_ERROR, "%s: For nested ordering, Nside must be a power of 2", ITER);
@ -313,8 +246,8 @@ static int iterate_healpix(grib_iterator_healpix* self, long N)
// Compute the longitudes at a given latitude
for (double longitude : HEALPix_longitudes(N, i)) {
Assert( ring_to_nest.at(j) < Npix );
self->lons[ring_to_nest.at(j)] = longitude;
self->lats[ring_to_nest.at(j)] = latitudes[i];
lons_[ring_to_nest.at(j)] = longitude;
lats_[ring_to_nest.at(j)] = latitudes[i];
++j;
}
}
@ -322,8 +255,8 @@ static int iterate_healpix(grib_iterator_healpix* self, long N)
for (size_t i = 0, j = 0; i < Ny; i++) {
// Compute the longitudes at a given latitude
for (double longitude : HEALPix_longitudes(N, i)) {
self->lons[j] = longitude;
self->lats[j] = latitudes[i];
lons_[j] = longitude;
lats_[j] = latitudes[i];
++j;
}
}
@ -332,15 +265,14 @@ static int iterate_healpix(grib_iterator_healpix* self, long N)
return GRIB_SUCCESS;
}
} // anonymous namespace
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
int grib_iterator_healpix_t::init(grib_handle* h, grib_arguments* args)
{
int err = 0;
grib_iterator_healpix* self = (grib_iterator_healpix*)iter;
grib_iterator_gen_t::init(h, args);
const char* snside = grib_arguments_get_name(h, args, self->carg++);
const char* sorder = grib_arguments_get_name(h, args, self->carg++);
int err = 0;
const char* snside = grib_arguments_get_name(h, args, carg_++);
const char* sorder = grib_arguments_get_name(h, args, carg_++);
long N = 0;
if ((err = grib_get_long_internal(h, snside, &N)) != GRIB_SUCCESS) {
@ -357,12 +289,12 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
return err;
}
self->nested = STR_EQUAL(ordering, "nested");
if (!STR_EQUAL(ordering, "ring") && !self->nested) {
nested_ = STR_EQUAL(ordering, "nested");
if (!STR_EQUAL(ordering, "ring") && !nested_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Only orderingConvention=(ring|nested) are supported", ITER);
return GRIB_GEOCALCULUS_PROBLEM;
}
// if (self->nested && N == 1) {
// if (nested && N == 1) {
// grib_context_log(h->context, GRIB_LOG_ERROR, "%s: For orderingConvention=nested, N must be greater than 1", ITER);
// return GRIB_GEOCALCULUS_PROBLEM;
// }
@ -372,58 +304,55 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
return GRIB_WRONG_GRID;
}
if (iter->nv != 12 * N * N) {
if (nv_ != 12 * N * N) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=12x%ldx%ld)",
ITER, iter->nv, N, N);
ITER, nv_, N, N);
return GRIB_WRONG_GRID;
}
self->lats = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
if (self->lats == nullptr) {
lats_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (lats_ == nullptr) {
return GRIB_OUT_OF_MEMORY;
}
self->lons = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
if (self->lons == nullptr) {
lons_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (lons_ == nullptr) {
return GRIB_OUT_OF_MEMORY;
}
try {
err = iterate_healpix(self, N);
err = iterate_healpix(N);
}
catch (...) {
return GRIB_INTERNAL_ERROR;
}
iter->e = -1;
e_ = -1;
return err;
}
static int next(grib_iterator* iter, double* lat, double* lon, double* val)
int grib_iterator_healpix_t::next(double* lat, double* lon, double* val)
{
auto* self = (grib_iterator_healpix*)iter;
if (iter->e >= static_cast<long>(iter->nv - 1)) {
if (e_ >= static_cast<long>(nv_ - 1)) {
return 0;
}
iter->e++;
e_++;
*lat = self->lats[iter->e];
*lon = self->lons[iter->e];
if (val != nullptr && iter->data != nullptr) {
*val = iter->data[iter->e];
*lat = lats_[e_];
*lon = lons_[e_];
if (val != nullptr && data_ != nullptr) {
*val = data_[e_];
}
return 1;
}
static int destroy(grib_iterator* iter)
int grib_iterator_healpix_t::destroy()
{
auto* self = (grib_iterator_healpix*)iter;
const auto* c = iter->h->context;
grib_context_free(context_, lats_);
grib_context_free(context_, lons_);
grib_context_free(c, self->lats);
grib_context_free(c, self->lons);
grib_iterator_gen_t::destroy();
return GRIB_SUCCESS;
}

35
src/iterator/grib_iterator_class_healpix.h Normal file
View File

@ -0,0 +1,35 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_gen.h"
class grib_iterator_healpix_t : public grib_iterator_gen_t
{
public:
grib_iterator_healpix_t() : grib_iterator_gen_t()
{
class_name_ = "healpix";
}
int init(grib_handle*, grib_arguments*) override;
int next(double *, double *, double *) override;
int destroy() override;
public:
double *lats_;
double *lons_;
long Nsides_;
bool nested_;
private:
int iterate_healpix(long N);
};

169
src/grib_iterator_class_lambert_azimuthal_equal_area.cc → src/iterator/grib_iterator_class_lambert_azimuthal_equal_area.cc
View File

@ -8,93 +8,26 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_iterator_class_lambert_azimuthal_equal_area.h"
#include <cmath>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_gen
IMPLEMENTS = destroy
IMPLEMENTS = init;next
MEMBERS = double *lats
MEMBERS = double *lons
MEMBERS = long Nj
END_CLASS_DEF
*/
grib_iterator_lambert_azimuthal_equal_area_t _grib_iterator_lambert_azimuthal_equal_area{};
grib_iterator* grib_iterator_lambert_azimuthal_equal_area = &_grib_iterator_lambert_azimuthal_equal_area;
/* 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);
static int destroy (grib_iterator* i);
typedef struct grib_iterator_lambert_azimuthal_equal_area{
grib_iterator it;
/* Members defined in gen */
int carg;
const char* missingValue;
/* Members defined in lambert_azimuthal_equal_area */
double *lats;
double *lons;
long Nj;
} grib_iterator_lambert_azimuthal_equal_area;
extern grib_iterator_class* grib_iterator_class_gen;
static grib_iterator_class _grib_iterator_class_lambert_azimuthal_equal_area = {
&grib_iterator_class_gen, /* super */
"lambert_azimuthal_equal_area", /* name */
sizeof(grib_iterator_lambert_azimuthal_equal_area),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
&next, /* Next Value */
0, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_lambert_azimuthal_equal_area = &_grib_iterator_class_lambert_azimuthal_equal_area;
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 */
#define ITER "Lambert azimuthal equal area Geoiterator"
static int next(grib_iterator* iter, double* lat, double* lon, double* val)
int grib_iterator_lambert_azimuthal_equal_area_t::next(double* lat, double* lon, double* val)
{
grib_iterator_lambert_azimuthal_equal_area* self = (grib_iterator_lambert_azimuthal_equal_area*)iter;
if ((long)iter->e >= (long)(iter->nv - 1))
if ((long)e_ >= (long)(nv_ - 1))
return 0;
iter->e++;
e_++;
*lat = self->lats[iter->e];
*lon = self->lons[iter->e];
if (val && iter->data) {
*val = iter->data[iter->e];
*lat = lats_[e_];
*lon = lons_[e_];
if (val && data_) {
*val = data_[e_];
}
return 1;
}
@ -161,7 +94,7 @@ static double pj_qsfn(double sinphi, double e, double one_es)
#define EPS10 1.e-10
static int init_oblate(grib_handle* h,
grib_iterator_lambert_azimuthal_equal_area* self,
grib_iterator_lambert_azimuthal_equal_area_t* self,
size_t nv, long nx, long ny,
double Dx, double Dy, double earthMinorAxisInMetres, double earthMajorAxisInMetres,
double latFirstInRadians, double lonFirstInRadians,
@ -231,18 +164,18 @@ static int init_oblate(grib_handle* h,
x0 = Q__xmf * b * cosb * sinlam;
/* Allocate latitude and longitude arrays */
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
self->lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
self->lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
lats = self->lats;
lons = self->lons;
lats = self->lats_;
lons = self->lons_;
/* Populate the lat and lon arrays */
{
@ -292,7 +225,7 @@ static int init_oblate(grib_handle* h,
}
static int init_sphere(grib_handle* h,
grib_iterator_lambert_azimuthal_equal_area* self,
grib_iterator_lambert_azimuthal_equal_area_t* self,
size_t nv, long nx, long ny,
double Dx, double Dy, double radius,
double latFirstInRadians, double lonFirstInRadians,
@ -317,18 +250,18 @@ static int init_sphere(grib_handle* h,
Dx = iScansNegatively == 0 ? Dx / 1000 : -Dx / 1000;
Dy = jScansPositively == 1 ? Dy / 1000 : -Dy / 1000;
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
self->lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
self->lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
lats = self->lats;
lons = self->lons;
lats = self->lats_;
lons = self->lons_;
/* compute xFirst,yFirst in metres */
sinphi = sin(latFirstInRadians);
@ -398,8 +331,10 @@ static int init_sphere(grib_handle* h,
return GRIB_SUCCESS;
}
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
int grib_iterator_lambert_azimuthal_equal_area_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gen_t::init(h, args);
int err = 0;
int is_oblate = 0;
double lonFirstInDegrees, latFirstInDegrees, lonFirstInRadians, latFirstInRadians, radius = 0;
@ -411,22 +346,19 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
long jScansPositively, jPointsAreConsecutive;
double earthMajorAxisInMetres = 0, earthMinorAxisInMetres = 0;
grib_iterator_lambert_azimuthal_equal_area* self =
(grib_iterator_lambert_azimuthal_equal_area*)iter;
const char* sradius = grib_arguments_get_name(h, args, self->carg++);
const char* snx = grib_arguments_get_name(h, args, self->carg++);
const char* sny = grib_arguments_get_name(h, args, self->carg++);
const char* slatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* slonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sstandardParallel = grib_arguments_get_name(h, args, self->carg++);
const char* scentralLongitude = grib_arguments_get_name(h, args, self->carg++);
const char* sDx = grib_arguments_get_name(h, args, self->carg++);
const char* sDy = grib_arguments_get_name(h, args, self->carg++);
const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++);
const char* sjScansPositively = grib_arguments_get_name(h, args, self->carg++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
const char* salternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
const char* sradius = grib_arguments_get_name(h, args, carg_++);
const char* snx = grib_arguments_get_name(h, args, carg_++);
const char* sny = grib_arguments_get_name(h, args, carg_++);
const char* slatFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* slonFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sstandardParallel = grib_arguments_get_name(h, args, carg_++);
const char* scentralLongitude = grib_arguments_get_name(h, args, carg_++);
const char* sDx = grib_arguments_get_name(h, args, carg_++);
const char* sDy = grib_arguments_get_name(h, args, carg_++);
const char* siScansNegatively = grib_arguments_get_name(h, args, carg_++);
const char* sjScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* salternativeRowScanning = grib_arguments_get_name(h, args, carg_++);
const double d2r = acos(0.0) / 90.0;
is_oblate = grib_is_earth_oblate(h);
@ -443,8 +375,8 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
if ((err = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS)
return err;
if (iter->nv != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, iter->nv, nx, ny);
if (nv_ != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, nv_, nx, ny);
return GRIB_WRONG_GRID;
}
if ((err = grib_get_double_internal(h, slatFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS)
@ -474,14 +406,14 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
standardParallelInRadians = standardParallelInDegrees * d2r;
if (is_oblate) {
err = init_oblate(h, self, iter->nv, nx, ny,
err = init_oblate(h, this, nv_, nx, ny,
Dx, Dy, earthMinorAxisInMetres, earthMajorAxisInMetres,
latFirstInRadians, lonFirstInRadians,
centralLongitudeInRadians, standardParallelInRadians,
iScansNegatively, jScansPositively, jPointsAreConsecutive);
}
else {
err = init_sphere(h, self, iter->nv, nx, ny,
err = init_sphere(h, this, nv_, nx, ny,
Dx, Dy, radius,
latFirstInRadians, lonFirstInRadians,
centralLongitudeInRadians, standardParallelInRadians,
@ -489,17 +421,18 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
}
if (err) return err;
iter->e = -1;
e_ = -1;
return GRIB_SUCCESS;
}
static int destroy(grib_iterator* iter)
int grib_iterator_lambert_azimuthal_equal_area_t::destroy()
{
grib_iterator_lambert_azimuthal_equal_area* self = (grib_iterator_lambert_azimuthal_equal_area*)iter;
const grib_context* c = iter->h->context;
const grib_context* c = h_->context;
grib_context_free(c, self->lats);
grib_context_free(c, self->lons);
grib_context_free(c, lats_);
grib_context_free(c, lons_);
grib_iterator_gen_t::destroy();
return GRIB_SUCCESS;
}

32
src/iterator/grib_iterator_class_lambert_azimuthal_equal_area.h Normal file
View File

@ -0,0 +1,32 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_gen.h"
class grib_iterator_lambert_azimuthal_equal_area_t : public grib_iterator_gen_t
{
public:
grib_iterator_lambert_azimuthal_equal_area_t() : grib_iterator_gen_t()
{
class_name_ = "lambert_azimuthal_equal_area";
}
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) override;
int previous(double*, double*, double*) override;
int destroy() override;
public:
double *lats_;
double *lons_;
long Nj_;
};

184
src/grib_iterator_class_lambert_conformal.cc → src/iterator/grib_iterator_class_lambert_conformal.cc
View File

@ -8,78 +8,11 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_iterator_class_lambert_conformal.h"
#include <cmath>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_gen
IMPLEMENTS = destroy
IMPLEMENTS = init;next
MEMBERS = double *lats
MEMBERS = double *lons
MEMBERS = long Nj
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);
static int destroy (grib_iterator* i);
typedef struct grib_iterator_lambert_conformal{
grib_iterator it;
/* Members defined in gen */
int carg;
const char* missingValue;
/* Members defined in lambert_conformal */
double *lats;
double *lons;
long Nj;
} grib_iterator_lambert_conformal;
extern grib_iterator_class* grib_iterator_class_gen;
static grib_iterator_class _grib_iterator_class_lambert_conformal = {
&grib_iterator_class_gen, /* super */
"lambert_conformal", /* name */
sizeof(grib_iterator_lambert_conformal),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
&next, /* Next Value */
0, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_lambert_conformal = &_grib_iterator_class_lambert_conformal;
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 */
grib_iterator_lambert_conformal_t _grib_iterator_lambert_conformal{};
grib_iterator* grib_iterator_lambert_conformal = &_grib_iterator_lambert_conformal;
#define ITER "Lambert conformal Geoiterator"
#define EPSILON 1.0e-10
@ -188,8 +121,7 @@ static void xy2lonlat(double radius, double n, double f, double rho0_bare, doubl
}
}
static int init_sphere(const grib_handle* h,
grib_iterator_lambert_conformal* self,
int grib_iterator_lambert_conformal_t::init_sphere(const grib_handle* h,
size_t nv, long nx, long ny,
double LoVInDegrees,
double Dx, double Dy, double radius,
@ -225,13 +157,13 @@ static int init_sphere(const grib_handle* h,
// Dy = jScansPositively == 1 ? Dy : -Dy;
// Allocate latitude and longitude arrays
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
@ -244,8 +176,8 @@ static int init_sphere(const grib_handle* h,
x = x0 + i * Dx;
xy2lonlat(radius, n, f, rho0_bare, LoVInRadians, x, y, &lonDeg, &latDeg);
lonDeg = normalise_longitude_in_degrees(lonDeg);
self->lons[index] = lonDeg;
self->lats[index] = latDeg;
lons_[index] = lonDeg;
lats_[index] = latDeg;
}
}
@ -278,8 +210,7 @@ static int init_sphere(const grib_handle* h,
}
// Oblate spheroid
static int init_oblate(const grib_handle* h,
grib_iterator_lambert_conformal* self,
int grib_iterator_lambert_conformal_t::init_oblate(const grib_handle* h,
size_t nv, long nx, long ny,
double LoVInDegrees,
double Dx, double Dy,
@ -349,13 +280,13 @@ static int init_oblate(const grib_handle* h,
y0 = -y0;
// Allocate latitude and longitude arrays
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
@ -388,8 +319,8 @@ static int init_oblate(const grib_handle* h,
if (err) {
grib_context_log(h->context, GRIB_LOG_ERROR,
"%s: Failed to compute the latitude angle, phi2, for the inverse", ITER);
grib_context_free(h->context, self->lats);
grib_context_free(h->context, self->lons);
grib_context_free(h->context, lats_);
grib_context_free(h->context, lons_);
return err;
}
} else {
@ -401,15 +332,17 @@ static int init_oblate(const grib_handle* h,
}
latDeg = latRad * RAD2DEG; // Convert to degrees
lonDeg = normalise_longitude_in_degrees(lonRad * RAD2DEG);
self->lons[index] = lonDeg;
self->lats[index] = latDeg;
lons_[index] = lonDeg;
lats_[index] = latDeg;
}
}
return GRIB_SUCCESS;
}
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
int grib_iterator_lambert_conformal_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gen_t::init(h, args);
int err = 0, is_oblate = 0;
long nx, ny, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning;
double LoVInDegrees, LaDInDegrees, Latin1InDegrees, Latin2InDegrees, latFirstInDegrees,
@ -418,24 +351,22 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
LaDInRadians;
double earthMajorAxisInMetres=0, earthMinorAxisInMetres=0;
grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)iter;
const char* sradius = grib_arguments_get_name(h, args, self->carg++);
const char* snx = grib_arguments_get_name(h, args, self->carg++);
const char* sny = grib_arguments_get_name(h, args, self->carg++);
const char* sLoVInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sLaDInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sLatin1InDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sLatin2InDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* slatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* slonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sradius = grib_arguments_get_name(h, args, carg_++);
const char* snx = grib_arguments_get_name(h, args, carg_++);
const char* sny = grib_arguments_get_name(h, args, carg_++);
const char* sLoVInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLaDInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLatin1InDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLatin2InDegrees = grib_arguments_get_name(h, args, carg_++);
const char* slatFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* slonFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
// Dx and Dy are in Metres
const char* sDx = grib_arguments_get_name(h, args, self->carg++);
const char* sDy = grib_arguments_get_name(h, args, self->carg++);
const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++);
const char* sjScansPositively = grib_arguments_get_name(h, args, self->carg++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
const char* salternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
const char* sDx = grib_arguments_get_name(h, args, carg_++);
const char* sDy = grib_arguments_get_name(h, args, carg_++);
const char* siScansNegatively = grib_arguments_get_name(h, args, carg_++);
const char* sjScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* salternativeRowScanning = grib_arguments_get_name(h, args, carg_++);
if ((err = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS) return err;
if ((err = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS) return err;
@ -449,8 +380,8 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
if ((err = grib_get_double_internal(h, sradius, &radius)) != GRIB_SUCCESS) return err;
}
if (iter->nv != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, iter->nv, nx, ny);
if (nv_ != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, nv_, nx, ny);
return GRIB_WRONG_GRID;
}
@ -497,14 +428,14 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
LoVInRadians = LoVInDegrees * DEG2RAD;
if (is_oblate) {
err = init_oblate(h, self, iter->nv, nx, ny,
err = init_oblate(h, nv_, nx, ny,
LoVInDegrees,
Dx, Dy, earthMinorAxisInMetres, earthMajorAxisInMetres,
latFirstInRadians, lonFirstInRadians,
LoVInRadians, Latin1InRadians, Latin2InRadians,
LaDInRadians);
} else {
err = init_sphere(h, self, iter->nv, nx, ny,
err = init_sphere(h, nv_, nx, ny,
LoVInDegrees,
Dx, Dy, radius,
latFirstInRadians, lonFirstInRadians,
@ -512,37 +443,36 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
}
if (err) return err;
iter->e = -1;
e_ = -1;
// Apply the scanning mode flags which may require data array to be transformed
err = transform_iterator_data(h->context, iter->data,
err = transform_iterator_data(h->context, data_,
iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning,
iter->nv, nx, ny);
nv_, nx, ny);
return err;
}
static int next(grib_iterator* iter, double* lat, double* lon, double* val)
int grib_iterator_lambert_conformal_t::next(double* lat, double* lon, double* val)
{
grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)iter;
if ((long)iter->e >= (long)(iter->nv - 1))
if ((long)e_ >= (long)(nv_ - 1))
return 0;
iter->e++;
e_++;
*lat = self->lats[iter->e];
*lon = self->lons[iter->e];
if (val && iter->data) {
*val = iter->data[iter->e];
*lat = lats_[e_];
*lon = lons_[e_];
if (val && data_) {
*val = data_[e_];
}
return 1;
}
static int destroy(grib_iterator* i)
int grib_iterator_lambert_conformal_t::destroy()
{
grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)i;
const grib_context* c = i->h->context;
const grib_context* c = h_->context;
grib_context_free(c, self->lats);
grib_context_free(c, self->lons);
grib_context_free(c, lats_);
grib_context_free(c, lons_);
grib_iterator_gen_t::destroy();
return GRIB_SUCCESS;
}

48
src/iterator/grib_iterator_class_lambert_conformal.h Normal file
View File

@ -0,0 +1,48 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_gen.h"
class grib_iterator_lambert_conformal_t : public grib_iterator_gen_t
{
public:
grib_iterator_lambert_conformal_t() : grib_iterator_gen_t()
{
class_name_ = "lambert_conformal";
}
int init(grib_handle*,grib_arguments*) override;
int next(double*, double*, double*) override;
int destroy() override;
private:
double* lats_;
double* lons_;
long Nj_;
int init_sphere(const grib_handle*,
size_t, long, long,
double,
double, double, double,
double, double,
double, double, double,
double);
int init_oblate(const grib_handle*,
size_t, long, long,
double,
double, double,
double, double,
double, double,
double, double, double,
double);
};

152
src/iterator/grib_iterator_class_latlon.cc Normal file
View File

@ -0,0 +1,152 @@
/*
* (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_iterator_class_latlon.h"
grib_iterator_latlon_t _grib_iterator_latlon{};
grib_iterator* grib_iterator_latlon = &_grib_iterator_latlon;
int grib_iterator_latlon_t::next(double* lat, double* lon, double* val)
{
/* GRIB-238: Support rotated lat/lon grids */
double ret_lat, ret_lon, ret_val=0;
if ((long)e_ >= (long)(nv_ - 1))
return 0;
e_++;
/* Assumptions:
* All rows scan in the same direction (alternativeRowScanning==0)
*/
if (!jPointsAreConsecutive_) {
/* Adjacent points in i (x) direction are consecutive */
ret_lat = las_[(long)floor(e_ / Ni_)];
ret_lon = los_[(long)e_ % Ni_];
if (data_)
ret_val = data_[e_];
}
else {
/* Adjacent points in j (y) direction is consecutive */
ret_lon = los_[(long)e_ / Nj_];
ret_lat = las_[(long)floor(e_ % Nj_)];
if (data_)
ret_val = data_[e_];
}
/* See ECC-808: Some users want to disable the unrotate */
if (isRotated_ && !disableUnrotate_) {
double new_lat = 0, new_lon = 0;
unrotate(ret_lat, ret_lon,
angleOfRotation_, southPoleLat_, southPoleLon_,
&new_lat, &new_lon);
ret_lat = new_lat;
ret_lon = new_lon;
}
*lat = ret_lat;
*lon = ret_lon;
if (val && data_) {
*val = ret_val;
}
return 1;
}
int grib_iterator_latlon_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_regular_t::init(h, args);
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, carg_++);
const char* s_jdir = grib_arguments_get_name(h, args, carg_++);
const char* s_jScansPos = grib_arguments_get_name(h, args, carg_++);
const char* s_jPtsConsec = grib_arguments_get_name(h, args, carg_++);
const char* s_isRotatedGrid = grib_arguments_get_name(h, args, carg_++);
const char* s_angleOfRotation = grib_arguments_get_name(h, args, carg_++);
const char* s_latSouthernPole = grib_arguments_get_name(h, args, carg_++);
const char* s_lonSouthernPole = grib_arguments_get_name(h, args, carg_++);
angleOfRotation_ = 0;
isRotated_ = 0;
southPoleLat_ = 0;
southPoleLon_ = 0;
disableUnrotate_ = 0; /* unrotate enabled by default */
if ((err = grib_get_long(h, s_isRotatedGrid, &isRotated_)))
return err;
if (isRotated_) {
if ((err = grib_get_double_internal(h, s_angleOfRotation, &angleOfRotation_)))
return err;
if ((err = grib_get_double_internal(h, s_latSouthernPole, &southPoleLat_)))
return err;
if ((err = grib_get_double_internal(h, s_lonSouthernPole, &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, &jPointsAreConsecutive_)))
return err;
if ((err = grib_get_long(h, "iteratorDisableUnrotate", &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 = 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 < Nj_; lai++) {
las_[lai] = lat1;
lat1 -= jdir;
}
/* ECC-1406: Due to rounding, errors can accumulate.
* So we ensure the last latitude is latitudeOfLastGridPointInDegrees
*/
las_[Nj_-1] = lat2;
e_ = -1;
return err;
}

25
src/iterator/grib_iterator_class_latlon.h Normal file
View File

@ -0,0 +1,25 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_regular.h"
class grib_iterator_latlon_t : public grib_iterator_regular_t
{
public:
grib_iterator_latlon_t() : grib_iterator_regular_t()
{
class_name_ = "latlon";
}
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) override;
};

139
src/iterator/grib_iterator_class_latlon_reduced.cc Normal file
View File

@ -0,0 +1,139 @@
/*
* (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_iterator_class_latlon_reduced.h"
#include <cmath>
grib_iterator_latlon_reduced_t _grib_iterator_latlon_reduced{};
grib_iterator* grib_iterator_latlon_reduced = &_grib_iterator_latlon_reduced;
int grib_iterator_latlon_reduced_t::next(double* lat, double* lon, double* val)
{
if ((long)e_ >= (long)(nv_ - 1))
return 0;
e_++;
*lat = las_[e_];
*lon = los_[e_];
if (val && data_) {
*val = data_[e_];
}
return 1;
}
int grib_iterator_latlon_reduced_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gen_t::init(h, args);
int ret = GRIB_SUCCESS;
double laf;
double lal;
long nlats;
double lof, tlof;
double lol, dimin;
long* pl;
size_t plsize = 0;
long k, j, ii;
long nlons, plmax;
double jdirinc = 0;
double idirinc = 0;
double dlon = 0;
int islocal = 0;
long nlons2 = 0; /* adjusted num of longitudes */
const char* latofirst = grib_arguments_get_name(h, args, carg_++);
const char* longoffirst = grib_arguments_get_name(h, args, carg_++);
const char* latoflast = grib_arguments_get_name(h, args, carg_++);
const char* longoflast = grib_arguments_get_name(h, args, carg_++);
const char* nlats_name = grib_arguments_get_name(h, args, carg_++);
const char* jdirec = grib_arguments_get_name(h, args, carg_++);
const char* plac = grib_arguments_get_name(h, args, carg_++);
if ((ret = grib_get_double_internal(h, latofirst, &laf)))
return ret;
if ((ret = grib_get_double_internal(h, longoffirst, &lof)))
return ret;
if ((ret = grib_get_double_internal(h, latoflast, &lal)))
return ret;
if ((ret = grib_get_double_internal(h, longoflast, &lol)))
return ret;
if ((ret = grib_get_long_internal(h, nlats_name, &nlats)))
return ret;
if ((ret = grib_get_double_internal(h, jdirec, &jdirinc)))
return ret;
plsize = nlats;
pl = (long*)grib_context_malloc(h->context, plsize * sizeof(long));
grib_get_long_array_internal(h, plac, pl, &plsize);
las_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
los_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
plmax = pl[0];
for (j = 0; j < nlats; j++)
if (plmax < pl[j])
plmax = pl[j];
dimin = 360.0 / plmax;
if (360 - fabs(lol - lof) < 2 * dimin) {
dlon = 360;
islocal = 0;
}
else if (lol < lof) {
/* handle something like 150 to -120 to treat as 150 to 240 */
/* so that dlon is 90 (not -270) */
dlon = lol + 360.0 - lof;
islocal = 1;
}
else {
dlon = lol - lof;
islocal = 1;
}
if (laf > lal)
jdirinc = -jdirinc;
k = 0;
for (j = 0; j < nlats; j++) {
nlons = pl[j];
tlof = lof;
nlons2 = nlons - islocal;
/*Sometimes there are no points on a latitude! Protect against div by zero*/
if (nlons2 < 1)
nlons2 = 1;
idirinc = dlon / nlons2;
for (ii = 0; ii < nlons; ii++) {
las_[k] = laf;
los_[k] = tlof;
tlof += idirinc;
k++;
}
laf += jdirinc;
}
e_ = -1;
grib_context_free(h->context, pl);
return ret;
}
int grib_iterator_latlon_reduced_t::destroy()
{
const grib_context* c = h_->context;
grib_context_free(c, las_);
grib_context_free(c, los_);
grib_iterator_gen_t::destroy();
return GRIB_SUCCESS;
}

30
src/iterator/grib_iterator_class_latlon_reduced.h Normal file
View File

@ -0,0 +1,30 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_gen.h"
class grib_iterator_latlon_reduced_t : public grib_iterator_gen_t
{
public:
grib_iterator_latlon_reduced_t() : grib_iterator_gen_t()
{
class_name_ = "latlon_reduced";
}
int init(grib_handle*,grib_arguments*) override;
int next(double *lat, double *lon, double *val) override;
int destroy() override;
private:
double* las_;
double* los_;
};

167
src/grib_iterator_class_mercator.cc → src/iterator/grib_iterator_class_mercator.cc
View File

@ -8,78 +8,11 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_iterator_class_mercator.h"
#include <cmath>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_gen
IMPLEMENTS = destroy
IMPLEMENTS = init;next
MEMBERS = double *lats
MEMBERS = double *lons
MEMBERS = long Nj
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);
static int destroy (grib_iterator* i);
typedef struct grib_iterator_mercator{
grib_iterator it;
/* Members defined in gen */
int carg;
const char* missingValue;
/* Members defined in mercator */
double *lats;
double *lons;
long Nj;
} grib_iterator_mercator;
extern grib_iterator_class* grib_iterator_class_gen;
static grib_iterator_class _grib_iterator_class_mercator = {
&grib_iterator_class_gen, /* super */
"mercator", /* name */
sizeof(grib_iterator_mercator),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
&next, /* Next Value */
0, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_mercator = &_grib_iterator_class_mercator;
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 */
grib_iterator_mercator_t _grib_iterator_mercator{};
grib_iterator* grib_iterator_mercator = &_grib_iterator_mercator;
#define ITER "Mercator Geoiterator"
#define EPSILON 1.0e-10
@ -148,8 +81,7 @@ static double compute_t(
return (tan(0.5 * (M_PI_2 - phi)) / con);
}
static int init_mercator(grib_handle* h,
grib_iterator_mercator* self,
int grib_iterator_mercator_t::init_mercator(grib_handle* h,
size_t nv, long nx, long ny,
double DiInMetres, double DjInMetres,
double earthMinorAxisInMetres, double earthMajorAxisInMetres,
@ -184,13 +116,13 @@ static int init_mercator(grib_handle* h,
y0 = -y0;
/* Allocate latitude and longitude arrays */
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
@ -211,8 +143,8 @@ static int init_mercator(grib_handle* h,
latRad = compute_phi(e, ts, &err);
if (err) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Failed to compute the latitude angle, phi2, for the inverse", ITER);
grib_context_free(h->context, self->lats);
grib_context_free(h->context, self->lons);
grib_context_free(h->context, lats_);
grib_context_free(h->context, lons_);
return err;
}
lonRad = adjust_lon_radians(orientationInRadians + _x / (earthMajorAxisInMetres * m1));
@ -221,15 +153,17 @@ static int init_mercator(grib_handle* h,
}
latDeg = latRad * RAD2DEG; /* Convert to degrees */
lonDeg = normalise_longitude_in_degrees(lonRad * RAD2DEG);
self->lons[index] = lonDeg;
self->lats[index] = latDeg;
lons_[index] = lonDeg;
lats_[index] = latDeg;
}
}
return GRIB_SUCCESS;
}
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
int grib_iterator_mercator_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gen_t::init(h, args);
int err = 0;
long ni, nj, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning;
double latFirstInDegrees, lonFirstInDegrees, LaDInDegrees;
@ -238,24 +172,22 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
LaDInRadians, orientationInRadians;
double earthMajorAxisInMetres = 0, earthMinorAxisInMetres = 0;
grib_iterator_mercator* self = (grib_iterator_mercator*)iter;
const char* sRadius = grib_arguments_get_name(h, args, self->carg++);
const char* sNi = grib_arguments_get_name(h, args, self->carg++);
const char* sNj = grib_arguments_get_name(h, args, self->carg++);
const char* sLatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sLonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sLaDInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sLatLastInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sLonLastInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sOrientationInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sRadius = grib_arguments_get_name(h, args, carg_++);
const char* sNi = grib_arguments_get_name(h, args, carg_++);
const char* sNj = grib_arguments_get_name(h, args, carg_++);
const char* sLatFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLonFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLaDInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLatLastInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLonLastInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sOrientationInDegrees = grib_arguments_get_name(h, args, carg_++);
/* Dx and Dy are in Metres */
const char* sDi = grib_arguments_get_name(h, args, self->carg++);
const char* sDj = grib_arguments_get_name(h, args, self->carg++);
const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++);
const char* sjScansPositively = grib_arguments_get_name(h, args, self->carg++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
const char* sAlternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
const char* sDi = grib_arguments_get_name(h, args, carg_++);
const char* sDj = grib_arguments_get_name(h, args, carg_++);
const char* siScansNegatively = grib_arguments_get_name(h, args, carg_++);
const char* sjScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* sAlternativeRowScanning = grib_arguments_get_name(h, args, carg_++);
if ((err = grib_get_long_internal(h, sNi, &ni)) != GRIB_SUCCESS) return err;
if ((err = grib_get_long_internal(h, sNj, &nj)) != GRIB_SUCCESS) return err;
@ -269,8 +201,8 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
earthMinorAxisInMetres = earthMajorAxisInMetres = radius;
}
if (iter->nv != ni * nj) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, iter->nv, ni, nj);
if (nv_ != ni * nj) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, nv_, ni, nj);
return GRIB_WRONG_GRID;
}
@ -307,44 +239,43 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
LaDInRadians = LaDInDegrees * DEG2RAD;
orientationInRadians = orientationInDegrees * DEG2RAD;
err = init_mercator(h, self, iter->nv, ni, nj,
err = init_mercator(h, nv_, ni, nj,
DiInMetres, DjInMetres, earthMinorAxisInMetres, earthMajorAxisInMetres,
latFirstInRadians, lonFirstInRadians,
latLastInRadians, lonLastInRadians,
LaDInRadians, orientationInRadians);
if (err) return err;
iter->e = -1;
e_ = -1;
/* Apply the scanning mode flags which may require data array to be transformed */
err = transform_iterator_data(h->context, iter->data,
err = transform_iterator_data(h->context, data_,
iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning,
iter->nv, ni, nj);
nv_, ni, nj);
return err;
}
static int next(grib_iterator* iter, double* lat, double* lon, double* val)
int grib_iterator_mercator_t::next(double* lat, double* lon, double* val)
{
grib_iterator_mercator* self = (grib_iterator_mercator*)iter;
if ((long)iter->e >= (long)(iter->nv - 1))
if ((long)e_ >= (long)(nv_ - 1))
return 0;
iter->e++;
e_++;
*lat = self->lats[iter->e];
*lon = self->lons[iter->e];
if (val && iter->data) {
*val = iter->data[iter->e];
*lat = lats_[e_];
*lon = lons_[e_];
if (val && data_) {
*val = data_[e_];
}
return 1;
}
static int destroy(grib_iterator* iter)
int grib_iterator_mercator_t::destroy()
{
grib_iterator_mercator* self = (grib_iterator_mercator*)iter;
const grib_context* c = iter->h->context;
const grib_context* c = h_->context;
grib_context_free(c, self->lats);
grib_context_free(c, self->lons);
grib_context_free(c, lats_);
grib_context_free(c, lons_);
grib_iterator_gen_t::destroy();
return GRIB_SUCCESS;
}

38
src/iterator/grib_iterator_class_mercator.h Normal file
View File

@ -0,0 +1,38 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_gen.h"
class grib_iterator_mercator_t : public grib_iterator_gen_t
{
public:
grib_iterator_mercator_t() {
class_name_ = "mercator";
}
int init(grib_handle*,grib_arguments*) override;
int next(double *lat, double *lon, double *val) override;
int destroy() override;
private:
double *lats_;
double *lons_;
long Nj_;
int init_mercator(grib_handle*,
size_t, long, long,
double, double,
double, double,
double, double,
double, double,
double, double);
};

160
src/grib_iterator_class_polar_stereographic.cc → src/iterator/grib_iterator_class_polar_stereographic.cc
View File

@ -8,93 +8,24 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_iterator_class_polar_stereographic.h"
#include <cmath>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_gen
IMPLEMENTS = destroy
IMPLEMENTS = init;next
MEMBERS = double *lats
MEMBERS = double *lons
MEMBERS = long Nj
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);
static int destroy (grib_iterator* i);
typedef struct grib_iterator_polar_stereographic{
grib_iterator it;
/* Members defined in gen */
int carg;
const char* missingValue;
/* Members defined in polar_stereographic */
double *lats;
double *lons;
long Nj;
} grib_iterator_polar_stereographic;
extern grib_iterator_class* grib_iterator_class_gen;
static grib_iterator_class _grib_iterator_class_polar_stereographic = {
&grib_iterator_class_gen, /* super */
"polar_stereographic", /* name */
sizeof(grib_iterator_polar_stereographic),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
&next, /* Next Value */
0, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_polar_stereographic = &_grib_iterator_class_polar_stereographic;
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 */
grib_iterator_polar_stereographic_t _grib_iterator_polar_stereographic{};
grib_iterator* grib_iterator_polar_stereographic = &_grib_iterator_polar_stereographic;
#define ITER "Polar stereographic Geoiterator"
static int next(grib_iterator* iter, double* lat, double* lon, double* val)
int grib_iterator_polar_stereographic_t::next(double* lat, double* lon, double* val)
{
grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)iter;
if ((long)iter->e >= (long)(iter->nv - 1))
if ((long)e_ >= (long)(nv_ - 1))
return 0;
iter->e++;
e_++;
*lat = self->lats[iter->e];
*lon = self->lons[iter->e];
if (val && iter->data) {
*val = iter->data[iter->e];
*lat = lats_[e_];
*lon = lons_[e_];
if (val && data_) {
*val = data_[e_];
}
return 1;
}
@ -117,8 +48,10 @@ typedef struct proj_data_t
#define PI_OVER_2 1.5707963267948966 /* half pi */
#define EPSILON 1.0e-10
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
int grib_iterator_polar_stereographic_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gen_t::init(h, args);
int ret = 0;
double *lats, *lons; /* arrays for latitudes and longitudes */
double lonFirstInDegrees, latFirstInDegrees, radius;
@ -135,22 +68,20 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
proj_data_t fwd_proj_data = {0,};
proj_data_t inv_proj_data = {0,};
grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)iter;
const char* s_radius = grib_arguments_get_name(h, args, self->carg++);
const char* s_nx = grib_arguments_get_name(h, args, self->carg++);
const char* s_ny = grib_arguments_get_name(h, args, self->carg++);
const char* s_latFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* s_lonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* s_southPoleOnPlane = grib_arguments_get_name(h, args, self->carg++);
const char* s_centralLongitude = grib_arguments_get_name(h, args, self->carg++);
const char* s_centralLatitude = grib_arguments_get_name(h, args, self->carg++);
const char* s_Dx = grib_arguments_get_name(h, args, self->carg++);
const char* s_Dy = grib_arguments_get_name(h, args, self->carg++);
const char* s_iScansNegatively = grib_arguments_get_name(h, args, self->carg++);
const char* s_jScansPositively = grib_arguments_get_name(h, args, self->carg++);
const char* s_jPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
const char* s_alternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
const char* s_radius = grib_arguments_get_name(h, args, carg_++);
const char* s_nx = grib_arguments_get_name(h, args, carg_++);
const char* s_ny = grib_arguments_get_name(h, args, carg_++);
const char* s_latFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* s_lonFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* s_southPoleOnPlane = grib_arguments_get_name(h, args, carg_++);
const char* s_centralLongitude = grib_arguments_get_name(h, args, carg_++);
const char* s_centralLatitude = grib_arguments_get_name(h, args, carg_++);
const char* s_Dx = grib_arguments_get_name(h, args, carg_++);
const char* s_Dy = grib_arguments_get_name(h, args, carg_++);
const char* s_iScansNegatively = grib_arguments_get_name(h, args, carg_++);
const char* s_jScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* s_jPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* s_alternativeRowScanning = grib_arguments_get_name(h, args, carg_++);
if (grib_is_earth_oblate(h)) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Only supported for spherical earth.", ITER);
@ -164,8 +95,8 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
if ((ret = grib_get_long_internal(h, s_ny, &ny)) != GRIB_SUCCESS)
return ret;
if (iter->nv != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, iter->nv, nx, ny);
if (nv_ != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, nv_, nx, ny);
return GRIB_WRONG_GRID;
}
if ((ret = grib_get_double_internal(h, s_latFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS)
@ -243,18 +174,18 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
inv_proj_data.mcs = cos(con1);
inv_proj_data.tcs = tan(0.5 * (PI_OVER_2 - con1));
}
self->lats = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
if (!self->lats) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, iter->nv * sizeof(double));
lats_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv_ * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
self->lons = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
if (!self->lats) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, iter->nv * sizeof(double));
lons_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv_ * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
lats = self->lats;
lons = self->lons;
lats = lats_;
lons = lons_;
/* These will be processed later in transform_iterator_data() */
/* Dx = iScansNegatively == 0 ? Dx : -Dx; */
/* Dy = jScansPositively == 1 ? Dy : -Dy; */
@ -353,22 +284,23 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
// }
// }
iter->e = -1;
e_ = -1;
/* Apply the scanning mode flags which may require data array to be transformed */
ret = transform_iterator_data(h->context, iter->data,
ret = transform_iterator_data(h->context, data_,
iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning,
iter->nv, nx, ny);
nv_, nx, ny);
return ret;
}
static int destroy(grib_iterator* i)
int grib_iterator_polar_stereographic_t::destroy()
{
grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)i;
const grib_context* c = i->h->context;
const grib_context* c = h_->context;
grib_context_free(c, self->lats);
grib_context_free(c, self->lons);
grib_context_free(c, lats_);
grib_context_free(c, lons_);
grib_iterator_gen_t::destroy();
return GRIB_SUCCESS;
}

32
src/iterator/grib_iterator_class_polar_stereographic.h Normal file
View File

@ -0,0 +1,32 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_gen.h"
class grib_iterator_polar_stereographic_t : public grib_iterator_gen_t
{
public:
grib_iterator_polar_stereographic_t() : grib_iterator_gen_t()
{
class_name_ = "polar_stereographic";
}
int init(grib_handle*,grib_arguments*) override;
int next(double *, double *, double *) override;
int destroy() override;
private:
double* lats_;
double* lons_;
long Nj_;
};

160
src/iterator/grib_iterator_class_regular.cc Normal file
View File

@ -0,0 +1,160 @@
/*
* (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_iterator_class_regular.h"
#include <cmath>
grib_iterator_regular_t _grib_iterator_regular{};
grib_iterator_regular_t* grib_iterator_regular = &_grib_iterator_regular;
#define ITER "Regular grid Geoiterator"
int grib_iterator_regular_t::next(double* lat, double* lon, double* val)
{
if ((long)e_ >= (long)(nv_ - 1))
return 0;
e_++;
*lat = las_[(long)floor(e_ / Ni_)];
*lon = los_[(long)e_ % Ni_];
if (val && data_) {
*val = data_[e_];
}
return 1;
}
int grib_iterator_regular_t::previous(double* lat, double* lon, double* val)
{
if (e_ < 0)
return 0;
*lat = las_[(long)floor(e_ / Ni_)];
*lon = los_[e_ % Ni_];
if (val && data_) {
*val = data_[e_];
}
e_--;
return 1;
}
int grib_iterator_regular_t::destroy()
{
const grib_context* c = h_->context;
grib_context_free(c, las_);
grib_context_free(c, los_);
grib_iterator_gen_t::destroy();
return GRIB_SUCCESS;
}
int grib_iterator_regular_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gen_t::init(h, args);
int ret = GRIB_SUCCESS;
long Ni; /* Number of points along a parallel = Nx */
long Nj; /* Number of points along a meridian = Ny */
double idir, idir_coded, lon1, lon2;
long loi;
const char* s_lon1 = grib_arguments_get_name(h, args, carg_++);
const char* s_idir = grib_arguments_get_name(h, args, carg_++);
const char* s_Ni = grib_arguments_get_name(h, args, carg_++);
const char* s_Nj = grib_arguments_get_name(h, args, carg_++);
const char* s_iScansNeg = grib_arguments_get_name(h, args, carg_++);
if ((ret = grib_get_double_internal(h, s_lon1, &lon1)))
return ret;
if ((ret = grib_get_double_internal(h, "longitudeOfLastGridPointInDegrees", &lon2)))
return ret;
if ((ret = grib_get_double_internal(h, s_idir, &idir))) // can be GRIB_MISSING_DOUBLE
return ret;
idir_coded = idir;
if ((ret = grib_get_long_internal(h, s_Ni, &Ni)))
return ret;
if (grib_is_missing(h, s_Ni, &ret) && ret == GRIB_SUCCESS) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Key %s cannot be 'missing' for a regular grid!", ITER, s_Ni);
return GRIB_WRONG_GRID;
}
if ((ret = grib_get_long_internal(h, s_Nj, &Nj)))
return ret;
if (grib_is_missing(h, s_Nj, &ret) && ret == GRIB_SUCCESS) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Key %s cannot be 'missing' for a regular grid!", ITER, s_Nj);
return GRIB_WRONG_GRID;
}
if (Ni*Nj != nv_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Ni*Nj!=numberOfDataPoints (%ld*%ld!=%zu)", ITER, Ni, Nj, nv_);
return GRIB_WRONG_GRID;
}
if ((ret = grib_get_long_internal(h, s_iScansNeg, &iScansNegatively_)))
return ret;
/* GRIB-801: Careful of case with a single point! Ni==1 */
if (Ni > 1) {
/* Note: If first and last longitudes are equal I assume you wanna go round the globe */
if (iScansNegatively_) {
if (lon1 > lon2) {
idir = (lon1 - lon2) / (Ni - 1);
}
else {
idir = (lon1 + 360.0 - lon2) / (Ni - 1);
}
}
else {
if (lon2 > lon1) {
idir = (lon2 - lon1) / (Ni - 1);
}
else {
idir = (lon2 + 360.0 - lon1) / (Ni - 1);
}
}
}
if (iScansNegatively_) {
idir = -idir;
}
else {
if (lon1 + (Ni - 2) * idir > 360)
lon1 -= 360;
/*See ECC-704, GRIB-396*/
/*else if ( (lon1+(Ni-1)*idir)-360 > epsilon ){
idir=360.0/(float)Ni;
}*/
}
Ni_ = Ni;
Nj_ = Nj;
las_ = (double*)grib_context_malloc(h->context, Nj * sizeof(double));
los_ = (double*)grib_context_malloc(h->context, Ni * sizeof(double));
if (idir != idir_coded) {
grib_context_log(h->context, GRIB_LOG_DEBUG, "%s: Using idir=%g (coded value=%g)", ITER, idir, idir_coded);
}
for (loi = 0; loi < Ni; loi++) {
los_[loi] = lon1;
lon1 += idir;
}
// ECC-1406: Due to rounding, errors can accumulate.
// So we ensure the last longitude is longitudeOfLastGridPointInDegrees
// Also see ECC-1671, ECC-1708
if (lon2 > 0) {
lon2 = normalise_longitude_in_degrees(lon2);
}
los_[Ni-1] = lon2;
return ret;
}

40
src/iterator/grib_iterator_class_regular.h Normal file
View File

@ -0,0 +1,40 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_gen.h"
class grib_iterator_regular_t : public grib_iterator_gen_t
{
public:
grib_iterator_regular_t()
{
class_name_ = "regular";
}
int init(grib_handle*,grib_arguments*) override;
int next(double *lat, double *lon, double *val) override;
int previous(double *lat, double *lon, double *val) override;
int destroy() override;
protected:
long Ni_;
long Nj_;
double *las_;
double *los_;
long iScansNegatively_;
long isRotated_;
double angleOfRotation_;
double southPoleLat_;
double southPoleLon_;
long jPointsAreConsecutive_;
long disableUnrotate_;
};

166
src/grib_iterator_class_space_view.cc → src/iterator/grib_iterator_class_space_view.cc
View File

@ -8,93 +8,24 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_iterator_class_space_view.h"
#include <cmath>
/*
This is used by make_class.pl
START_CLASS_DEF
CLASS = iterator
SUPER = grib_iterator_class_gen
IMPLEMENTS = destroy
IMPLEMENTS = init;next
MEMBERS = double *lats
MEMBERS = double *lons
MEMBERS = long Nj
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);
static int destroy (grib_iterator* i);
typedef struct grib_iterator_space_view{
grib_iterator it;
/* Members defined in gen */
int carg;
const char* missingValue;
/* Members defined in space_view */
double *lats;
double *lons;
long Nj;
} grib_iterator_space_view;
extern grib_iterator_class* grib_iterator_class_gen;
static grib_iterator_class _grib_iterator_class_space_view = {
&grib_iterator_class_gen, /* super */
"space_view", /* name */
sizeof(grib_iterator_space_view),/* size of instance */
0, /* inited */
&init_class, /* init_class */
&init, /* constructor */
&destroy, /* destructor */
&next, /* Next Value */
0, /* Previous Value */
0, /* Reset the counter */
0, /* has next values */
};
grib_iterator_class* grib_iterator_class_space_view = &_grib_iterator_class_space_view;
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 */
grib_iterator_space_view_t _grib_iterator_space_view{};
grib_iterator* grib_iterator_space_view = &_grib_iterator_space_view;
#define ITER "Space view Geoiterator"
static int next(grib_iterator* iter, double* lat, double* lon, double* val)
int grib_iterator_space_view_t::next(double* lat, double* lon, double* val)
{
grib_iterator_space_view* self = (grib_iterator_space_view*)iter;
if ((long)iter->e >= (long)(iter->nv - 1))
if ((long)e_ >= (long)(nv_ - 1))
return 0;
iter->e++;
e_++;
*lat = self->lats[iter->e];
*lon = self->lons[iter->e];
if (val && iter->data) {
*val = iter->data[iter->e];
*lat = lats_[e_];
*lon = lons_[e_];
if (val && data_) {
*val = data_[e_];
}
return 1;
}
@ -142,8 +73,10 @@ static int next(grib_iterator* iter, double* lat, double* lon, double* val)
#define RAD2DEG 57.29577951308232087684 /* 180 over pi */
#define DEG2RAD 0.01745329251994329576 /* pi over 180 */
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
int grib_iterator_space_view_t::init(grib_handle* h, grib_arguments* args)
{
grib_iterator_gen_t::init(h, args);
/* REFERENCE:
* LRIT/HRIT Global Specification (CGMS 03, Issue 2.6, 12.08.1999)
*/
@ -163,31 +96,29 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
double factor_1, factor_2, tmp1, Sd, Sn, Sxy, S1, S2, S3;
int x0, y0, ix, iy;
double *s_x, *c_x; /* arrays storing sin and cos values */
size_t array_size = (iter->nv * sizeof(double));
size_t array_size = (nv_ * sizeof(double));
grib_iterator_space_view* self = (grib_iterator_space_view*)iter;
const char* sradius = grib_arguments_get_name(h, args, carg_++);
const char* sEarthIsOblate = grib_arguments_get_name(h, args, carg_++);
const char* sMajorAxisInMetres = grib_arguments_get_name(h, args, carg_++);
const char* sMinorAxisInMetres = grib_arguments_get_name(h, args, carg_++);
const char* snx = grib_arguments_get_name(h, args, carg_++);
const char* sny = grib_arguments_get_name(h, args, carg_++);
const char* sLatOfSubSatellitePointInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLonOfSubSatellitePointInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sDx = grib_arguments_get_name(h, args, carg_++);
const char* sDy = grib_arguments_get_name(h, args, carg_++);
const char* sXpInGridLengths = grib_arguments_get_name(h, args, carg_++);
const char* sYpInGridLengths = grib_arguments_get_name(h, args, carg_++);
const char* sOrientationInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sNrInRadiusOfEarthScaled = grib_arguments_get_name(h, args, carg_++);
const char* sXo = grib_arguments_get_name(h, args, carg_++);
const char* sYo = grib_arguments_get_name(h, args, carg_++);
const char* sradius = grib_arguments_get_name(h, args, self->carg++);
const char* sEarthIsOblate = grib_arguments_get_name(h, args, self->carg++);
const char* sMajorAxisInMetres = grib_arguments_get_name(h, args, self->carg++);
const char* sMinorAxisInMetres = grib_arguments_get_name(h, args, self->carg++);
const char* snx = grib_arguments_get_name(h, args, self->carg++);
const char* sny = grib_arguments_get_name(h, args, self->carg++);
const char* sLatOfSubSatellitePointInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sLonOfSubSatellitePointInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sDx = grib_arguments_get_name(h, args, self->carg++);
const char* sDy = grib_arguments_get_name(h, args, self->carg++);
const char* sXpInGridLengths = grib_arguments_get_name(h, args, self->carg++);
const char* sYpInGridLengths = grib_arguments_get_name(h, args, self->carg++);
const char* sOrientationInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sNrInRadiusOfEarthScaled = grib_arguments_get_name(h, args, self->carg++);
const char* sXo = grib_arguments_get_name(h, args, self->carg++);
const char* sYo = grib_arguments_get_name(h, args, self->carg++);
const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++);
const char* sjScansPositively = grib_arguments_get_name(h, args, self->carg++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
const char* sAlternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
const char* siScansNegatively = grib_arguments_get_name(h, args, carg_++);
const char* sjScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* sAlternativeRowScanning = grib_arguments_get_name(h, args, carg_++);
if ((ret = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS)
return ret;
@ -207,8 +138,8 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
return ret;
}
if (iter->nv != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, iter->nv, nx, ny);
if (nv_ != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Wrong number of points (%zu!=%ldx%ld)", ITER, nv_, nx, ny);
return GRIB_WRONG_GRID;
}
if ((ret = grib_get_double_internal(h, sLatOfSubSatellitePointInDegrees, &latOfSubSatellitePointInDegrees)) != GRIB_SUCCESS)
@ -288,18 +219,18 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
rx = angular_size / dx;
ry = (r_pol / r_eq) * angular_size / dy;
self->lats = (double*)grib_context_malloc(h->context, array_size);
if (!self->lats) {
lats_ = (double*)grib_context_malloc(h->context, array_size);
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, array_size);
return GRIB_OUT_OF_MEMORY;
}
self->lons = (double*)grib_context_malloc(h->context, array_size);
if (!self->lats) {
lons_ = (double*)grib_context_malloc(h->context, array_size);
if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, array_size);
return GRIB_OUT_OF_MEMORY;
}
lats = self->lats;
lons = self->lons;
lats = lats_;
lons = lons_;
if (!iScansNegatively) {
xp = xp - x0;
@ -373,17 +304,18 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
}
grib_context_free(h->context, s_x);
grib_context_free(h->context, c_x);
iter->e = -1;
e_ = -1;
return ret;
}
static int destroy(grib_iterator* iter)
int grib_iterator_space_view_t::destroy()
{
grib_iterator_space_view* self = (grib_iterator_space_view*)iter;
const grib_context* c = iter->h->context;
const grib_context* c = h_->context;
grib_context_free(c, self->lats);
grib_context_free(c, self->lons);
grib_context_free(c, lats_);
grib_context_free(c, lons_);
grib_iterator_gen_t::destroy();
return GRIB_SUCCESS;
}

32
src/iterator/grib_iterator_class_space_view.h Normal file
View File

@ -0,0 +1,32 @@
/*
* (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.
*/
#pragma once
#include "grib_iterator_class_gen.h"
class grib_iterator_space_view_t : public grib_iterator_gen_t
{
public:
grib_iterator_space_view_t()
{
class_name_ = "space_view";
}
int init(grib_handle*, grib_arguments*) override;
int next(double *, double *, double *) override;
int destroy() override;
private:
double* lats_;
double* lons_;
long Nj_;
};