Merge pull request #256 from ecmwf/feature/modernisation_iterators

ECC-1948: Conversion of geo-iterators to C++
This commit is contained in:
shahramn 2024-10-21 22:24:47 +01:00 committed by GitHub
commit bc034fc726
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
40 changed files with 2046 additions and 2359 deletions

View File

@ -434,7 +434,10 @@ endforeach()
### include directories ### include directories
include_directories( ${ECCODES_INCLUDE_DIRS} ${ECCODES_EXTRA_INCLUDE_DIRS} ) include_directories(
${ECCODES_INCLUDE_DIRS}
${ECCODES_EXTRA_INCLUDE_DIRS}
)
add_subdirectory( definitions ) # must be before memfs add_subdirectory( definitions ) # must be before memfs
add_subdirectory( memfs ) add_subdirectory( memfs )

View File

@ -8,6 +8,12 @@
# granted to it by virtue of its status as an intergovernmental organisation # granted to it by virtue of its status as an intergovernmental organisation
# nor does it submit to any jurisdiction. # nor does it submit to any jurisdiction.
# #
include_directories(
"${CMAKE_CURRENT_SOURCE_DIR}/accessor"
"${CMAKE_CURRENT_SOURCE_DIR}/geo_iterator"
)
list( APPEND eccodes_src_files list( APPEND eccodes_src_files
accessor/grib_accessor.cc accessor/grib_accessor.cc
accessor/grib_accessors_list.cc accessor/grib_accessors_list.cc
@ -331,20 +337,20 @@ list( APPEND eccodes_src_files
grib_nearest_class_mercator.cc grib_nearest_class_mercator.cc
grib_nearest_class_polar_stereographic.cc grib_nearest_class_polar_stereographic.cc
grib_nearest_class_space_view.cc grib_nearest_class_space_view.cc
grib_iterator_class_polar_stereographic.cc geo_iterator/grib_iterator_class_polar_stereographic.cc
grib_iterator_class_lambert_azimuthal_equal_area.cc geo_iterator/grib_iterator_class_lambert_azimuthal_equal_area.cc
grib_iterator_class_lambert_conformal.cc geo_iterator/grib_iterator_class_lambert_conformal.cc
grib_iterator_class_mercator.cc geo_iterator/grib_iterator_class_mercator.cc
grib_iterator.cc geo_iterator/grib_iterator.cc
grib_iterator_class.cc grib_iterator_class.cc
grib_iterator_class_gaussian.cc geo_iterator/grib_iterator_class_gaussian.cc
grib_iterator_class_gaussian_reduced.cc geo_iterator/grib_iterator_class_gaussian_reduced.cc
grib_iterator_class_latlon_reduced.cc geo_iterator/grib_iterator_class_latlon_reduced.cc
grib_iterator_class_gen.cc geo_iterator/grib_iterator_class_gen.cc
grib_iterator_class_latlon.cc geo_iterator/grib_iterator_class_latlon.cc
grib_iterator_class_regular.cc geo_iterator/grib_iterator_class_regular.cc
grib_iterator_class_space_view.cc geo_iterator/grib_iterator_class_space_view.cc
grib_iterator_class_healpix.cc geo_iterator/grib_iterator_class_healpix.cc
grib_expression.cc grib_expression.cc
codes_util.cc codes_util.cc
grib_util.cc grib_util.cc

View File

@ -24,35 +24,3 @@ void grib_accessor_iterator_t::dump(grib_dumper* dumper)
/* TODO: pass args */ /* TODO: pass args */
grib_dump_label(dumper, this, NULL); grib_dump_label(dumper, this, NULL);
} }
#if defined(HAVE_GEOGRAPHY)
grib_iterator* grib_iterator_new(const grib_handle* ch, unsigned long flags, int* error)
{
grib_handle* h = (grib_handle*)ch;
grib_accessor* a = NULL;
grib_accessor_iterator_t* ita = NULL;
grib_iterator* iter = NULL;
*error = GRIB_NOT_IMPLEMENTED;
a = grib_find_accessor(h, "ITERATOR");
ita = (grib_accessor_iterator_t*)a;
if (!a)
return NULL;
iter = grib_iterator_factory(h, ita->args_, flags, error);
if (iter)
*error = GRIB_SUCCESS;
return iter;
}
#else
grib_iterator* grib_iterator_new(const grib_handle* ch, unsigned long flags, int* error)
{
*error = GRIB_FUNCTIONALITY_NOT_ENABLED;
grib_context_log(ch->context, GRIB_LOG_ERROR,
"Geoiterator functionality not enabled. Please rebuild with -DENABLE_GEOGRAPHY=ON");
return NULL;
}
#endif

View File

@ -23,8 +23,5 @@ public:
private: private:
grib_arguments* args_; grib_arguments* args_;
friend eccodes::geo_iterator::Iterator* eccodes::geo_iterator::gribIteratorNew(const grib_handle*, unsigned long, int*);
friend grib_iterator* grib_iterator_new(const grib_handle* ch, unsigned long flags, int* error);
}; };
// grib_iterator* grib_iterator_new(const grib_handle* ch, unsigned long flags, int* error)

View File

@ -812,15 +812,15 @@ grib_nearest* grib_nearest_factory(grib_handle* h, grib_arguments* args, int* er
/* grib_iterator.cc */ /* grib_iterator.cc */
int grib_get_data(const grib_handle* h, double* lats, double* lons, double* values); 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_next(grib_iterator* i, double* lat, double* lon, double* value);
int grib_iterator_has_next(grib_iterator* i); //int grib_iterator_has_next(grib_iterator* i);
int grib_iterator_previous(grib_iterator* i, double* lat, double* lon, double* value); //int grib_iterator_previous(grib_iterator* i, double* lat, double* lon, double* value);
int grib_iterator_reset(grib_iterator* i); //int grib_iterator_reset(grib_iterator* i);
int grib_iterator_init(grib_iterator* i, grib_handle* h, grib_arguments* args); //int grib_iterator_init(grib_iterator* i, grib_handle* h, grib_arguments* args);
int grib_iterator_delete(grib_iterator* i); //int grib_iterator_delete(grib_iterator* i);
/* grib_iterator_class.cc */ /* grib_iterator_class.cc */
grib_iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsigned long flags, int* error); eccodes::geo_iterator::Iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsigned long flags, int* error);
/* grib_iterator_class_gen.cc */ /* grib_iterator_class_gen.cc */
int transform_iterator_data(grib_context* c, double* data, long iScansNegatively, long jScansPositively, long jPointsAreConsecutive, long alternativeRowScanning, size_t numPoints, long nx, long ny); int transform_iterator_data(grib_context* c, double* data, long iScansNegatively, long jScansPositively, long jPointsAreConsecutive, long alternativeRowScanning, size_t numPoints, long nx, long ny);

View File

@ -0,0 +1,133 @@
/*
* (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_iterator.h"
#include "accessor/grib_accessor_class_iterator.h"
namespace eccodes::geo_iterator {
int Iterator::init(grib_handle* h, grib_arguments* args)
{
h_ = h;
return GRIB_SUCCESS;
}
/* For this one, ALL destroy are called */
int Iterator::destroy()
{
delete context_;
delete this;
return GRIB_SUCCESS;
}
eccodes::geo_iterator::Iterator* gribIteratorNew(const grib_handle* ch, unsigned long flags, int* error)
{
grib_handle* h = (grib_handle*)ch;
grib_accessor* a = NULL;
grib_accessor_iterator_t* ita = NULL;
*error = GRIB_NOT_IMPLEMENTED;
a = grib_find_accessor(h, "ITERATOR");
ita = (grib_accessor_iterator_t*)a;
if (!a)
return NULL;
eccodes::geo_iterator::Iterator* iter = grib_iterator_factory(h, ita->args_, flags, error);
if (iter)
*error = GRIB_SUCCESS;
return iter;
}
int gribIteratorDelete(eccodes::geo_iterator::Iterator* i)
{
if (i)
i->destroy();
return GRIB_SUCCESS;
}
} // namespace eccodes::geo_iterator
/*
* C API Implementation
* codes_iterator_* wrappers are in eccodes.h and eccodes.cc
* grib_iterator_* declarations are in grib_api.h
*/
int grib_iterator_reset(grib_iterator* i)
{
return i->iterator->reset();
}
int grib_iterator_has_next(grib_iterator* i)
{
return i->iterator->has_next();
}
int grib_iterator_next(grib_iterator* i, double* lat, double* lon, double* value)
{
return i->iterator->next(lat, lon, value);
}
int grib_iterator_previous(grib_iterator* i, double* lat, double* lon, double* value)
{
return i->iterator->previous(lat, lon, value);
}
int grib_iterator_destroy(grib_context* c, grib_iterator* i)
{
return i->iterator->destroy();
}
#if defined(HAVE_GEOGRAPHY)
grib_iterator* grib_iterator_new(const grib_handle* ch, unsigned long flags, int* error)
{
grib_iterator* i = (grib_iterator*)grib_context_malloc_clear(ch->context, sizeof(grib_iterator));
i->iterator = eccodes::geo_iterator::gribIteratorNew(ch, flags, error);
if (!i->iterator) {
grib_context_free(ch->context, i);
return NULL;
}
return i;
}
int grib_iterator_delete(grib_iterator* i)
{
if (i) {
grib_context* c = grib_context_get_default();
gribIteratorDelete(i->iterator);
grib_context_free(c, i);
}
return GRIB_SUCCESS;
}
#else
grib_iterator* grib_iterator_new(const grib_handle* ch, unsigned long flags, int* error)
{
*error = GRIB_FUNCTIONALITY_NOT_ENABLED;
grib_context_log(ch->context, GRIB_LOG_ERROR,
"Geoiterator functionality not enabled. Please rebuild with -DENABLE_GEOGRAPHY=ON");
return NULL;
}
int grib_iterator_delete(grib_iterator* i)
{
grib_context_log(ch->context, GRIB_LOG_ERROR,
"Geoiterator functionality not enabled. Please rebuild with -DENABLE_GEOGRAPHY=ON");
return GRIB_FUNCTIONALITY_NOT_ENABLED;
}
#endif

View File

@ -0,0 +1,52 @@
/*
* (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
*/
namespace eccodes::geo_iterator {
class Iterator
{
public:
virtual ~Iterator() {}
virtual int init(grib_handle*, grib_arguments*) = 0;
virtual int next(double*, double*, double*) const = 0;
virtual int previous(double*, double*, double*) const = 0;
virtual int reset() = 0;
virtual int destroy() = 0;
virtual bool has_next() const = 0;
virtual Iterator* create() const = 0;
unsigned long flags_;
protected:
grib_context* context_;
grib_handle* h_;
double* data_; /** data values */
mutable long e_; /** current element */
size_t nv_; /** number of values */
const char* class_name_;
private:
//grib_arguments* args_; [>* args of iterator <]
size_t size_;
int inited_;
};
eccodes::geo_iterator::Iterator* gribIteratorNew(const grib_handle*, unsigned long, int*);
int gribIteratorDelete(eccodes::geo_iterator::Iterator*);
} // namespace eccodes::geo_iterator

View File

@ -8,105 +8,35 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. * 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>
/* eccodes::geo_iterator::Gaussian _grib_iterator_gaussian{};
This is used by make_class.pl eccodes::geo_iterator::Iterator* grib_iterator_gaussian = &_grib_iterator_gaussian;
START_CLASS_DEF namespace eccodes::geo_iterator {
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 */
static void binary_search_gaussian_latitudes(const double xx[], const unsigned long n, double x, long* j); 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 Gaussian::init(grib_handle* h, grib_arguments* args)
{ {
grib_iterator_gaussian* self = (grib_iterator_gaussian*)i; int ret = GRIB_SUCCESS;
if ((ret = Regular::init(h, args)) != GRIB_SUCCESS)
return ret;
double* lats; double* lats;
double laf; /* latitude of first point in degrees */ double laf; /* latitude of first point in degrees */
double lal; /* latitude of last point in degrees */ double lal; /* latitude of last point in degrees */
long trunc; /* number of parallels between a pole and the equator */ long trunc; /* number of parallels between a pole and the equator */
long lai = 0; long lai = 0;
long jScansPositively = 0; long jScansPositively = 0;
int size = 0; int size = 0;
double start; double start;
long istart = 0; long istart = 0;
int ret = GRIB_SUCCESS;
const char* latofirst = 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, self->carg++); const char* latoflast = grib_arguments_get_name(h, args, carg_++);
const char* numtrunc = grib_arguments_get_name(h, args, self->carg++); const char* numtrunc = grib_arguments_get_name(h, args, carg_++);
const char* s_jScansPositively = grib_arguments_get_name(h, args, self->carg++); const char* s_jScansPositively = grib_arguments_get_name(h, args, carg_++);
if ((ret = grib_get_double_internal(h, latofirst, &laf))) if ((ret = grib_get_double_internal(h, latofirst, &laf)))
return ret; return ret;
@ -137,22 +67,22 @@ static int init(grib_iterator* i, grib_handle* h, grib_arguments* args)
} }
*/ */
binary_search_gaussian_latitudes(lats, size-1, start, &istart); binary_search_gaussian_latitudes(lats, size - 1, start, &istart);
if (istart < 0 || istart >= size) { if (istart < 0 || istart >= size) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Failed to find index for latitude=%g", start); grib_context_log(h->context, GRIB_LOG_ERROR, "Failed to find index for latitude=%g", start);
return GRIB_GEOCALCULUS_PROBLEM; return GRIB_GEOCALCULUS_PROBLEM;
} }
if (jScansPositively) { if (jScansPositively) {
for (lai = 0; lai < self->Nj; lai++) { for (lai = 0; lai < Nj_; lai++) {
DEBUG_ASSERT(istart >= 0); DEBUG_ASSERT(istart >= 0);
self->las[lai] = lats[istart--]; lats_[lai] = lats[istart--];
if (istart<0) istart=size-1; if (istart < 0) istart = size - 1;
} }
} }
else { else {
for (lai = 0; lai < self->Nj; lai++) { for (lai = 0; lai < Nj_; lai++) {
self->las[lai] = lats[istart++]; lats_[lai] = lats[istart++];
if (istart > size - 1) if (istart > size - 1)
istart = 0; istart = 0;
} }
@ -167,7 +97,7 @@ static int init(grib_iterator* i, grib_handle* h, grib_arguments* args)
/* Note: the argument 'n' is NOT the size of the 'xx' array but its LAST index i.e. size of xx - 1 */ /* Note: the argument 'n' is NOT the size of the 'xx' array but its LAST index i.e. size of xx - 1 */
static void binary_search_gaussian_latitudes(const double array[], const unsigned long n, double x, long* j) static void binary_search_gaussian_latitudes(const double array[], const unsigned long n, double x, long* j)
{ {
unsigned long low = 0; unsigned long low = 0;
unsigned long high = n; unsigned long high = n;
unsigned long mid; unsigned long mid;
const int descending = (array[n] < array[0]); const int descending = (array[n] < array[0]);
@ -216,3 +146,5 @@ static void binary_search_gaussian_latitudes(const double array[], const unsigne
// } // }
// *j = jl; // *j = jl;
// } // }
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,28 @@
/*
* (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"
namespace eccodes::geo_iterator {
class Gaussian : public Regular
{
public:
Gaussian() :
Regular() { class_name_ = "gaussian"; }
Iterator* create() const override { return new Gaussian(); }
int init(grib_handle*, grib_arguments*) override;
};
} // namespace eccodes::geo_iterator

View File

@ -8,111 +8,34 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. * 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> #include <cmath>
/* eccodes::geo_iterator::GaussianReduced _grib_iterator_gaussian_reduced{};
This is used by make_class.pl eccodes::geo_iterator::Iterator* grib_iterator_gaussian_reduced = &_grib_iterator_gaussian_reduced;
START_CLASS_DEF namespace eccodes::geo_iterator {
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 */
#define ITER "Reduced Gaussian grid Geoiterator" #define ITER "Reduced Gaussian grid Geoiterator"
static int next(grib_iterator* iter, double* lat, double* lon, double* val) int GaussianReduced::next(double* lat, double* lon, double* val) const
{ {
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter; double ret_lat = 0, ret_lon = 0;
double ret_lat=0, ret_lon=0;
if (iter->e >= (long)(iter->nv - 1)) if (e_ >= (long)(nv_ - 1))
return 0; return 0;
iter->e++; e_++;
ret_lat = self->las[iter->e]; ret_lat = lats_[e_];
ret_lon = self->los[iter->e]; ret_lon = lons_[e_];
if (val && iter->data) { if (val && data_) {
*val = iter->data[iter->e]; *val = data_[e_];
} }
if (self->isRotated && !self->disableUnrotate) { if (isRotated_ && !disableUnrotate_) {
double new_lat = 0, new_lon = 0; double new_lat = 0, new_lon = 0;
unrotate(ret_lat, ret_lon, unrotate(ret_lat, ret_lon,
self->angleOfRotation, self->southPoleLat, self->southPoleLon, angleOfRotation_, southPoleLat_, southPoleLon_,
&new_lat, &new_lon); &new_lat, &new_lon);
ret_lat = new_lat; ret_lat = new_lat;
ret_lon = new_lon; ret_lon = new_lon;
@ -172,10 +95,10 @@ static void binary_search(const double xx[], const unsigned long n, double x, lo
} }
/* Use legacy way to compute the iterator latitude/longitude values */ /* Use legacy way to compute the iterator latitude/longitude values */
static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_handle* h, int GaussianReduced::iterate_reduced_gaussian_subarea_legacy(grib_handle* h,
double lat_first, double lon_first, double lat_first, double lon_first,
double lat_last, double lon_last, double lat_last, double lon_last,
double* lats, long* pl, size_t plsize) double* lats, long* pl, size_t plsize)
{ {
int err = 0; int err = 0;
int l = 0; int l = 0;
@ -183,9 +106,8 @@ static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_han
long row_count = 0; long row_count = 0;
double d = 0; double d = 0;
long ilon_first, ilon_last, i; 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;*/
get_reduced_row_proc get_reduced_row = &grib_get_reduced_row_legacy; /* legacy algorithm */ get_reduced_row_proc get_reduced_row = &grib_get_reduced_row_legacy; /* legacy algorithm */
if (h->context->debug) { if (h->context->debug) {
const size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last); const size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last);
@ -198,7 +120,7 @@ static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_han
l++; l++;
} }
iter->e = 0; e_ = 0;
for (j = 0; j < plsize; j++) { for (j = 0; j < plsize; j++) {
long k = 0; long k = 0;
row_count = 0; row_count = 0;
@ -207,16 +129,16 @@ static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_han
if (ilon_first > ilon_last) if (ilon_first > ilon_last)
ilon_first -= pl[j]; ilon_first -= pl[j];
for (i = ilon_first; i <= ilon_last; i++) { 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); size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last);
grib_context_log(h->context, GRIB_LOG_ERROR, 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; return GRIB_WRONG_GRID;
} }
self->los[iter->e] = ((i)*360.0) / pl[j]; lons_[e_] = ((i) * 360.0) / pl[j];
self->las[iter->e] = lats[j + l]; lats_[e_] = lats[j + l];
iter->e++; e_++;
k++; k++;
if (k >= row_count) { if (k >= row_count) {
/* Ensure we exit the loop and only process 'row_count' points */ /* Ensure we exit the loop and only process 'row_count' points */
@ -235,17 +157,16 @@ 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); // err = iterate_reduced_gaussian_subarea_algorithm2(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize);
// } // }
// return err; // return err;
static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h, int GaussianReduced::iterate_reduced_gaussian_subarea(grib_handle* h,
double lat_first, double lon_first, double lat_first, double lon_first,
double lat_last, double lon_last, double lat_last, double lon_last,
double* lats, long* pl, size_t plsize, size_t numlats) double* lats, long* pl, size_t plsize, size_t numlats)
{ {
int err = 0; int err = 0;
long l = 0; long l = 0;
size_t j = 0; size_t j = 0;
long row_count = 0, i = 0; long row_count = 0, i = 0;
double olon_first, olon_last; 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; get_reduced_row_proc get_reduced_row = &grib_get_reduced_row;
if (h->context->debug) { if (h->context->debug) {
@ -257,41 +178,41 @@ static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
binary_search(lats, numlats - 1, lat_first, &l); binary_search(lats, numlats - 1, lat_first, &l);
Assert(l < numlats); Assert(l < numlats);
// for(il=0; il<numlats; ++il) { // for(il=0; il<numlats; ++il) {
// const double diff = fabs(lat_first-lats[il]); // const double diff = fabs(lat_first-lats[il]);
// if (diff < min_d) { // if (diff < min_d) {
// min_d = diff; // min_d = diff;
// l = il; /* index of the latitude */ // l = il; /* index of the latitude */
// } // }
// } // }
iter->e = 0; e_ = 0;
for (j = 0; j < plsize; j++) { for (j = 0; j < plsize; j++) {
const double delta = 360.0 / pl[j]; const double delta = 360.0 / pl[j];
row_count = 0; row_count = 0;
grib_get_reduced_row_p(pl[j], lon_first, lon_last, &row_count, &olon_first, &olon_last); grib_get_reduced_row_p(pl[j], lon_first, lon_last, &row_count, &olon_first, &olon_last);
for (i = 0; i < row_count; ++i) { for (i = 0; i < row_count; ++i) {
double lon2 = olon_first + i * delta; double lon2 = olon_first + i * delta;
if (iter->e >= iter->nv) { if (e_ >= nv_) {
/* Only print error message on the second pass */ /* Only print error message on the second pass */
size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last); size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last);
grib_context_log(h->context, GRIB_LOG_ERROR, 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; return GRIB_WRONG_GRID;
} }
self->los[iter->e] = lon2; lons_[e_] = lon2;
DEBUG_ASSERT(j + l < numlats); DEBUG_ASSERT(j + l < numlats);
self->las[iter->e] = lats[j + l]; lats_[e_] = lats[j + l];
iter->e++; e_++;
} }
} }
if (iter->e != iter->nv) { if (e_ != nv_) {
/* Fewer counted points in the sub-area than the number of data values */ /* 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); 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) */ /* 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 { else {
/* TODO: A gap exists! Not all values can be mapped. Inconsistent grid or error in calculating num. points! */ /* TODO: A gap exists! Not all values can be mapped. Inconsistent grid or error in calculating num. points! */
@ -300,9 +221,13 @@ static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
return err; return err;
} }
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args) int GaussianReduced::init(grib_handle* h, grib_arguments* args)
{ {
int ret = GRIB_SUCCESS, j, is_global = 0; int ret = GRIB_SUCCESS;
if ((ret = Gen::init(h, args)) != GRIB_SUCCESS)
return ret;
int j, is_global = 0;
double lat_first = 0, lon_first = 0, lat_last = 0, lon_last = 0; double lat_first = 0, lon_first = 0, lat_last = 0, lon_last = 0;
double angular_precision = 1.0 / 1000000.0; double angular_precision = 1.0 / 1000000.0;
double* lats; double* lats;
@ -311,31 +236,30 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
long* pl; long* pl;
long max_pl = 0; long max_pl = 0;
long nj = 0, order = 0, i; long nj = 0, order = 0, i;
long row_count = 0; long row_count = 0;
long angleSubdivisions = 0; long angleSubdivisions = 0;
const grib_context* c = h->context; 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, carg_++);
const char* slat_first = grib_arguments_get_name(h, args, self->carg++); const char* slon_first = grib_arguments_get_name(h, args, carg_++);
const char* slon_first = grib_arguments_get_name(h, args, self->carg++); const char* slat_last = grib_arguments_get_name(h, args, carg_++);
const char* slat_last = grib_arguments_get_name(h, args, self->carg++); const char* slon_last = grib_arguments_get_name(h, args, carg_++);
const char* slon_last = grib_arguments_get_name(h, args, self->carg++); const char* sorder = grib_arguments_get_name(h, args, carg_++);
const char* sorder = grib_arguments_get_name(h, args, self->carg++); const char* spl = grib_arguments_get_name(h, args, carg_++);
const char* spl = grib_arguments_get_name(h, args, self->carg++); const char* snj = grib_arguments_get_name(h, args, carg_++);
const char* snj = grib_arguments_get_name(h, args, self->carg++);
self->angleOfRotation = 0; angleOfRotation_ = 0;
self->isRotated = 0; isRotated_ = 0;
self->southPoleLat = 0; southPoleLat_ = 0;
self->southPoleLon = 0; southPoleLon_ = 0;
self->disableUnrotate = 0; /* unrotate enabled by default */ disableUnrotate_ = 0; /* unrotate enabled by default */
ret = grib_get_long(h, "isRotatedGrid", &self->isRotated); ret = grib_get_long(h, "isRotatedGrid", &isRotated_);
if (ret == GRIB_SUCCESS && self->isRotated) { if (ret == GRIB_SUCCESS && isRotated_) {
if ((ret = grib_get_double_internal(h, "angleOfRotation", &self->angleOfRotation))) if ((ret = grib_get_double_internal(h, "angleOfRotation", &angleOfRotation_)))
return ret; return ret;
if ((ret = grib_get_double_internal(h, "latitudeOfSouthernPoleInDegrees", &self->southPoleLat))) if ((ret = grib_get_double_internal(h, "latitudeOfSouthernPoleInDegrees", &southPoleLat_)))
return ret; return ret;
if ((ret = grib_get_double_internal(h, "longitudeOfSouthernPoleInDegrees", &self->southPoleLon))) if ((ret = grib_get_double_internal(h, "longitudeOfSouthernPoleInDegrees", &southPoleLon_)))
return ret; return ret;
} }
@ -379,11 +303,11 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
grib_get_long_array_internal(h, spl, pl, &plsize); grib_get_long_array_internal(h, spl, pl, &plsize);
self->las = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double)); lats_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (!self->las) if (!lats_)
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
self->los = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double)); lons_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (!self->los) if (!lons_)
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
while (lon_last < 0) while (lon_last < 0)
@ -402,11 +326,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); is_global = is_gaussian_global(lat_first, lat_last, lon_first, lon_last, max_pl, lats, angular_precision);
if (!is_global) { if (!is_global) {
/*sub area*/ /*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 { else {
/*global*/ /*global*/
iter->e = 0; e_ = 0;
if (h->context->debug) { if (h->context->debug) {
const size_t np = count_global_points(pl, plsize); const size_t np = count_global_points(pl, plsize);
fprintf(stderr, "ECCODES DEBUG grib_iterator_class_gaussian_reduced: global num points=%zu\n", np); fprintf(stderr, "ECCODES DEBUG grib_iterator_class_gaussian_reduced: global num points=%zu\n", np);
@ -415,37 +339,39 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
for (j = 0; j < plsize; j++) { for (j = 0; j < plsize; j++) {
row_count = pl[j]; row_count = pl[j];
for (i = 0; i < row_count; i++) { 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)");*/ /*grib_context_log(h->context,GRIB_LOG_ERROR, "Failed to initialise reduced Gaussian iterator (global)");*/
/*return GRIB_WRONG_GRID;*/ /*return GRIB_WRONG_GRID;*/
/*Try now as NON-global*/ /*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) if (ret != GRIB_SUCCESS)
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Failed to initialise iterator (global)", ITER); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Failed to initialise iterator (global)", ITER);
goto finalise; goto finalise;
} }
self->los[iter->e] = (i * 360.0) / row_count; lons_[e_] = (i * 360.0) / row_count;
self->las[iter->e] = lats[j]; lats_[e_] = lats[j];
iter->e++; e_++;
} }
} }
} }
finalise: finalise:
iter->e = -1; e_ = -1;
grib_context_free(h->context, lats); grib_context_free(h->context, lats);
grib_context_free(h->context, pl); grib_context_free(h->context, pl);
return ret; return ret;
} }
static int destroy(grib_iterator* iter) int GaussianReduced::destroy()
{ {
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter; const grib_context* c = h_->context;
const grib_context* c = iter->h->context;
grib_context_free(c, self->las); grib_context_free(c, lats_);
grib_context_free(c, self->los); grib_context_free(c, lons_);
return GRIB_SUCCESS;
return Gen::destroy();
} }
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,47 @@
/*
* (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"
namespace eccodes::geo_iterator {
class GaussianReduced : public Gen
{
public:
GaussianReduced() :
Gen() { class_name_ = "gaussian_reduced"; }
Iterator* create() const override { return new GaussianReduced(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int destroy() override;
private:
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);
};
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,120 @@
/*
* (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_gen.h"
namespace eccodes::geo_iterator {
int Gen::init(grib_handle* h, grib_arguments* args)
{
int err = GRIB_SUCCESS;
if ((err = Iterator::init(h, args)) != GRIB_SUCCESS)
return err;
size_t dli = 0;
const char* s_rawData = NULL;
const char* s_numPoints = NULL;
long numberOfPoints = 0;
carg_ = 1;
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_++);
data_ = NULL;
if ((err = grib_get_size(h, s_rawData, &dli)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(h, s_numPoints, &numberOfPoints)) != GRIB_SUCCESS)
return err;
// 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 (flags_ & GRIB_GEOITERATOR_NO_VALUES) {
// Iterator's number of values taken from the Grid Section
nv_ = numberOfPoints;
}
else {
// Check for consistency between the Grid and Data Sections
if (numberOfPoints != dli) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Geoiterator: %s != size(%s) (%ld!=%ld)",
s_numPoints, s_rawData, numberOfPoints, dli);
return GRIB_WRONG_GRID;
}
nv_ = dli;
}
if (nv_ == 0) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Geoiterator: size(%s) is %ld", s_rawData, dli);
return GRIB_WRONG_GRID;
}
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
data_ = (double*)grib_context_malloc(h->context, (nv_) * sizeof(double));
if ((err = grib_get_double_array_internal(h, s_rawData, data_, &(nv_)))) {
return err;
}
}
e_ = -1;
return err;
}
int Gen::reset()
{
e_ = -1;
return 0;
}
int Gen::destroy()
{
const grib_context* c = h_->context;
grib_context_free(c, data_);
return Iterator::destroy();
}
bool Gen::has_next() const
{
if (flags_ == 0 && data_ == NULL)
return false;
if (e_ >= (long)(nv_ - 1))
return false;
return true;
}
int Gen::previous(double*, double*, double*) const
{
return GRIB_NOT_IMPLEMENTED;
}
int Gen::next(double*, double*, double*) const
{
return GRIB_NOT_IMPLEMENTED;
}
// int Gen::get(double* lat, double* lon, double* val)
//{
// if (e_ >= (long)(nv_ - 1))
// return GRIB_END_OF_ITERATION;
// e_++;
// if (lat) *lat = 0;
// if (lon) *lon = 0;
// if (val) *val = 0;
// return 1;
//}
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,41 @@
/*
* (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"
namespace eccodes::geo_iterator {
class Gen : public Iterator
{
public:
Gen() :
Iterator() { class_name_ = "gen"; }
Iterator* create() const override { return new Gen(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int previous(double*, double*, double*) const override;
int reset() override;
int destroy() override;
bool has_next() const override;
protected:
int carg_;
double* lats_;
double* lons_;
private:
const char* missingValue_;
// int get(double*, double*, double*);
};
} // namespace eccodes::geo_iterator

View File

@ -8,93 +8,22 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. * 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 <algorithm>
#include <bitset> #include <bitset>
#include <cmath> #include <cmath>
#include <cstdint> #include <cstdint>
#include <tuple>
#include <vector> #include <vector>
/* eccodes::geo_iterator::Healpix _grib_iterator_healpix;
This is used by make_class.pl eccodes::geo_iterator::Iterator* grib_iterator_healpix = &_grib_iterator_healpix;
START_CLASS_DEF namespace eccodes::geo_iterator {
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 */
#define ITER "HEALPix Geoiterator" #define ITER "HEALPix Geoiterator"
constexpr double RAD2DEG = 57.29577951308232087684; // 180 over pi constexpr double RAD2DEG = 57.29577951308232087684; // 180 over pi
namespace
{
struct CodecFijNest struct CodecFijNest
{ {
static constexpr uint64_t __masks[] = { 0x00000000ffffffff, 0x0000ffff0000ffff, 0x00ff00ff00ff00ff, static constexpr uint64_t __masks[] = { 0x00000000ffffffff, 0x0000ffff0000ffff, 0x00ff00ff00ff00ff,
@ -211,7 +140,7 @@ static std::vector<double> HEALPix_longitudes(size_t N, size_t i)
return longitudes; return longitudes;
} }
static int iterate_healpix(grib_iterator_healpix* self, long N, bool nested) int Healpix::iterate_healpix(long N)
{ {
size_t Ny = 4 * static_cast<size_t>(N) - 1; size_t Ny = 4 * static_cast<size_t>(N) - 1;
auto Nd = static_cast<double>(N); auto Nd = static_cast<double>(N);
@ -241,7 +170,7 @@ static int iterate_healpix(grib_iterator_healpix* self, long N, bool nested)
// Equator // Equator
latitudes[2 * N - 1] = 0.0; latitudes[2 * N - 1] = 0.0;
if (nested) { if (nested_) {
if (!is_power_of_2(N)) { if (!is_power_of_2(N)) {
grib_context* c = grib_context_get_default(); 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); grib_context_log(c, GRIB_LOG_ERROR, "%s: For nested ordering, Nside must be a power of 2", ITER);
@ -309,21 +238,22 @@ static int iterate_healpix(grib_iterator_healpix* self, long N, bool nested)
ring_to_nest[r] = to_nest(f, ring, Nside, phi, ring & 1); ring_to_nest[r] = to_nest(f, ring, Nside, phi, ring & 1);
} }
for (size_t i = 0, j=0; i < Ny; i++) {
// 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];
++j;
}
}
} else {
for (size_t i = 0, j = 0; i < Ny; i++) { for (size_t i = 0, j = 0; i < Ny; i++) {
// Compute the longitudes at a given latitude // Compute the longitudes at a given latitude
for (double longitude : HEALPix_longitudes(N, i)) { for (double longitude : HEALPix_longitudes(N, i)) {
self->lons[j] = longitude; Assert(ring_to_nest.at(j) < Npix);
self->lats[j] = latitudes[i]; lons_[ring_to_nest.at(j)] = longitude;
lats_[ring_to_nest.at(j)] = latitudes[i];
++j;
}
}
}
else {
for (size_t i = 0, j = 0; i < Ny; i++) {
// Compute the longitudes at a given latitude
for (double longitude : HEALPix_longitudes(N, i)) {
lons_[j] = longitude;
lats_[j] = latitudes[i];
++j; ++j;
} }
} }
@ -332,15 +262,14 @@ static int iterate_healpix(grib_iterator_healpix* self, long N, bool nested)
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
} // anonymous namespace int Healpix::init(grib_handle* h, grib_arguments* args)
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
{ {
int err = 0; int err = GRIB_SUCCESS;
grib_iterator_healpix* self = (grib_iterator_healpix*)iter; if ((err = Gen::init(h, args)) != GRIB_SUCCESS)
return err;
const char* snside = grib_arguments_get_name(h, args, self->carg++); const char* snside = grib_arguments_get_name(h, args, carg_++);
const char* sorder = grib_arguments_get_name(h, args, self->carg++); const char* sorder = grib_arguments_get_name(h, args, carg_++);
long N = 0; long N = 0;
if ((err = grib_get_long_internal(h, snside, &N)) != GRIB_SUCCESS) { if ((err = grib_get_long_internal(h, snside, &N)) != GRIB_SUCCESS) {
@ -351,18 +280,20 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
return GRIB_WRONG_GRID; return GRIB_WRONG_GRID;
} }
char ordering[32] = {0,}; char ordering[32] = {
0,
};
size_t slen = sizeof(ordering); size_t slen = sizeof(ordering);
if ((err = grib_get_string_internal(h, sorder, ordering, &slen)) != GRIB_SUCCESS) { if ((err = grib_get_string_internal(h, sorder, ordering, &slen)) != GRIB_SUCCESS) {
return err; return err;
} }
self->nested = STR_EQUAL(ordering, "nested"); nested_ = STR_EQUAL(ordering, "nested");
if (!STR_EQUAL(ordering, "ring") && !self->nested) { if (!STR_EQUAL(ordering, "ring") && !nested_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Only orderingConvention=(ring|nested) are supported", ITER); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Only orderingConvention=(ring|nested) are supported", ITER);
return GRIB_GEOCALCULUS_PROBLEM; 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); // grib_context_log(h->context, GRIB_LOG_ERROR, "%s: For orderingConvention=nested, N must be greater than 1", ITER);
// return GRIB_GEOCALCULUS_PROBLEM; // return GRIB_GEOCALCULUS_PROBLEM;
// } // }
@ -372,58 +303,56 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
return GRIB_WRONG_GRID; 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)", 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; return GRIB_WRONG_GRID;
} }
self->lats = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double)); lats_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (self->lats == nullptr) { if (lats_ == nullptr) {
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
self->lons = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double)); lons_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (self->lons == nullptr) { if (lons_ == nullptr) {
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
try { try {
err = iterate_healpix(self, N, self->nested); err = iterate_healpix(N);
} }
catch (...) { catch (...) {
return GRIB_INTERNAL_ERROR; return GRIB_INTERNAL_ERROR;
} }
iter->e = -1; e_ = -1;
return err; return err;
} }
static int next(grib_iterator* iter, double* lat, double* lon, double* val) int Healpix::next(double* lat, double* lon, double* val) const
{ {
auto* self = (grib_iterator_healpix*)iter; if (e_ >= static_cast<long>(nv_ - 1)) {
if (iter->e >= static_cast<long>(iter->nv - 1)) {
return 0; return 0;
} }
iter->e++; e_++;
*lat = self->lats[iter->e]; *lat = lats_[e_];
*lon = self->lons[iter->e]; *lon = lons_[e_];
if (val != nullptr && iter->data != nullptr) { if (val != nullptr && data_ != nullptr) {
*val = iter->data[iter->e]; *val = data_[e_];
} }
return 1; return 1;
} }
static int destroy(grib_iterator* iter) int Healpix::destroy()
{ {
auto* self = (grib_iterator_healpix*)iter; grib_context_free(context_, lats_);
const auto* c = iter->h->context; grib_context_free(context_, lons_);
grib_context_free(c, self->lats); return Gen::destroy();
grib_context_free(c, self->lons);
return GRIB_SUCCESS;
} }
} // namespace eccodes::geo_iterator

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"
namespace eccodes::geo_iterator {
class Healpix : public Gen
{
public:
Healpix() :
Gen() { class_name_ = "healpix"; }
Iterator* create() const override { return new Healpix(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int destroy() override;
private:
long Nsides_;
bool nested_;
int iterate_healpix(long N);
};
} // namespace eccodes::geo_iterator

View File

@ -8,107 +8,39 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. * 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>
/* eccodes::geo_iterator::LambertAzimuthalEqualArea _grib_iterator_lambert_azimuthal_equal_area{};
This is used by make_class.pl eccodes::geo_iterator::Iterator* grib_iterator_lambert_azimuthal_equal_area = &_grib_iterator_lambert_azimuthal_equal_area;
START_CLASS_DEF namespace eccodes::geo_iterator {
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_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" #define ITER "Lambert azimuthal equal area Geoiterator"
static int next(grib_iterator* iter, double* lat, double* lon, double* val) int LambertAzimuthalEqualArea::next(double* lat, double* lon, double* val) const
{ {
grib_iterator_lambert_azimuthal_equal_area* self = (grib_iterator_lambert_azimuthal_equal_area*)iter; if ((long)e_ >= (long)(nv_ - 1))
if ((long)iter->e >= (long)(iter->nv - 1))
return 0; return 0;
iter->e++; e_++;
*lat = self->lats[iter->e]; *lat = lats_[e_];
*lon = self->lons[iter->e]; *lon = lons_[e_];
if (val && iter->data) { if (val && data_) {
*val = iter->data[iter->e]; *val = data_[e_];
} }
return 1; return 1;
} }
#ifndef M_PI #ifndef M_PI
#define M_PI 3.14159265358979323846 /* Whole pie */ #define M_PI 3.14159265358979323846 /* Whole pie */
#endif #endif
#ifndef M_PI_2 #ifndef M_PI_2
#define M_PI_2 1.57079632679489661923 /* Half a pie */ #define M_PI_2 1.57079632679489661923 /* Half a pie */
#endif #endif
#ifndef M_PI_4 #ifndef M_PI_4
#define M_PI_4 0.78539816339744830962 /* Quarter of a pie */ #define M_PI_4 0.78539816339744830962 /* Quarter of a pie */
#endif #endif
#define RAD2DEG 57.29577951308232087684 /* 180 over pi */ #define RAD2DEG 57.29577951308232087684 /* 180 over pi */
@ -160,8 +92,7 @@ static double pj_qsfn(double sinphi, double e, double one_es)
} }
#define EPS10 1.e-10 #define EPS10 1.e-10
static int init_oblate(grib_handle* h, int LambertAzimuthalEqualArea::init_oblate(grib_handle* h,
grib_iterator_lambert_azimuthal_equal_area* self,
size_t nv, long nx, long ny, size_t nv, long nx, long ny,
double Dx, double Dy, double earthMinorAxisInMetres, double earthMajorAxisInMetres, double Dx, double Dy, double earthMinorAxisInMetres, double earthMajorAxisInMetres,
double latFirstInRadians, double lonFirstInRadians, double latFirstInRadians, double lonFirstInRadians,
@ -175,7 +106,9 @@ static int init_oblate(grib_handle* h,
double Q__qp = 0, Q__rq = 0, Q__cosb1, Q__sinb1, Q__dd, Q__xmf, Q__ymf, t; double Q__qp = 0, Q__rq = 0, Q__cosb1, Q__sinb1, Q__dd, Q__xmf, Q__ymf, t;
/* double Q__mmf = 0; */ /* double Q__mmf = 0; */
double e, es, temp, one_es; double e, es, temp, one_es;
double APA[3] = {0,}; double APA[3] = {
0,
};
double xFirst, yFirst; double xFirst, yFirst;
Dx = iScansNegatively == 0 ? Dx / 1000 : -Dx / 1000; Dx = iScansNegatively == 0 ? Dx / 1000 : -Dx / 1000;
@ -202,7 +135,7 @@ static int init_oblate(grib_handle* h,
else else
Q->mode = OBLIQ; Q->mode = OBLIQ;
*/ */
Q__qp = pj_qsfn(1.0, e, one_es); Q__qp = pj_qsfn(1.0, e, one_es);
/* Q__mmf = 0.5 / one_es; ---- TODO(masn): do I need this? */ /* Q__mmf = 0.5 / one_es; ---- TODO(masn): do I need this? */
pj_authset(es, APA); /* sets up APA array */ pj_authset(es, APA); /* sets up APA array */
Q__rq = sqrt(0.5 * Q__qp); Q__rq = sqrt(0.5 * Q__qp);
@ -211,10 +144,11 @@ static int init_oblate(grib_handle* h,
Q__cosb1 = sqrt(1.0 - Q__sinb1 * Q__sinb1); Q__cosb1 = sqrt(1.0 - Q__sinb1 * Q__sinb1);
if (Q__cosb1 == 0) { if (Q__cosb1 == 0) {
Q__dd = 1.0; Q__dd = 1.0;
} else { }
else {
Q__dd = cos(standardParallelInRadians) / (sqrt(1. - es * sinphi_ * sinphi_) * Q__rq * Q__cosb1); Q__dd = cos(standardParallelInRadians) / (sqrt(1. - es * sinphi_ * sinphi_) * Q__rq * Q__cosb1);
} }
Q__ymf = (Q__xmf = Q__rq) / Q__dd; Q__ymf = (Q__xmf = Q__rq) / Q__dd;
Q__xmf *= Q__dd; Q__xmf *= Q__dd;
sinb = q / Q__qp; sinb = q / Q__qp;
@ -231,18 +165,18 @@ static int init_oblate(grib_handle* h,
x0 = Q__xmf * b * cosb * sinlam; x0 = Q__xmf * b * cosb * sinlam;
/* Allocate latitude and longitude arrays */ /* Allocate latitude and longitude arrays */
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lons_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
lats = self->lats; lats = lats_;
lons = self->lons; lons = lons_;
/* Populate the lat and lon arrays */ /* Populate the lat and lon arrays */
{ {
@ -291,8 +225,7 @@ static int init_oblate(grib_handle* h,
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
static int init_sphere(grib_handle* h, int LambertAzimuthalEqualArea::init_sphere(grib_handle* h,
grib_iterator_lambert_azimuthal_equal_area* self,
size_t nv, long nx, long ny, size_t nv, long nx, long ny,
double Dx, double Dy, double radius, double Dx, double Dy, double radius,
double latFirstInRadians, double lonFirstInRadians, double latFirstInRadians, double lonFirstInRadians,
@ -315,20 +248,20 @@ static int init_sphere(grib_handle* h,
cosphi1 = cos(phi1); cosphi1 = cos(phi1);
sinphi1 = sin(phi1); sinphi1 = sin(phi1);
Dx = iScansNegatively == 0 ? Dx / 1000 : -Dx / 1000; Dx = iScansNegatively == 0 ? Dx / 1000 : -Dx / 1000;
Dy = jScansPositively == 1 ? Dy / 1000 : -Dy / 1000; Dy = jScansPositively == 1 ? Dy / 1000 : -Dy / 1000;
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lons_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
lats = self->lats; lats = lats_;
lons = self->lons; lons = lons_;
/* compute xFirst,yFirst in metres */ /* compute xFirst,yFirst in metres */
sinphi = sin(latFirstInRadians); sinphi = sin(latFirstInRadians);
@ -398,9 +331,12 @@ static int init_sphere(grib_handle* h,
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args) int LambertAzimuthalEqualArea::init(grib_handle* h, grib_arguments* args)
{ {
int err = 0; int err = 0;
if ((err = Gen::init(h, args)) != GRIB_SUCCESS)
return err;
int is_oblate = 0; int is_oblate = 0;
double lonFirstInDegrees, latFirstInDegrees, lonFirstInRadians, latFirstInRadians, radius = 0; double lonFirstInDegrees, latFirstInDegrees, lonFirstInRadians, latFirstInRadians, radius = 0;
long nx, ny; long nx, ny;
@ -411,22 +347,19 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
long jScansPositively, jPointsAreConsecutive; long jScansPositively, jPointsAreConsecutive;
double earthMajorAxisInMetres = 0, earthMinorAxisInMetres = 0; double earthMajorAxisInMetres = 0, earthMinorAxisInMetres = 0;
grib_iterator_lambert_azimuthal_equal_area* self = const char* sradius = grib_arguments_get_name(h, args, carg_++);
(grib_iterator_lambert_azimuthal_equal_area*)iter; const char* snx = grib_arguments_get_name(h, args, carg_++);
const char* sny = grib_arguments_get_name(h, args, carg_++);
const char* sradius = grib_arguments_get_name(h, args, self->carg++); const char* slatFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* snx = grib_arguments_get_name(h, args, self->carg++); const char* slonFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sny = grib_arguments_get_name(h, args, self->carg++); const char* sstandardParallel = grib_arguments_get_name(h, args, carg_++);
const char* slatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* scentralLongitude = grib_arguments_get_name(h, args, carg_++);
const char* slonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* sDx = grib_arguments_get_name(h, args, carg_++);
const char* sstandardParallel = grib_arguments_get_name(h, args, self->carg++); const char* sDy = grib_arguments_get_name(h, args, carg_++);
const char* scentralLongitude = grib_arguments_get_name(h, args, self->carg++); const char* siScansNegatively = grib_arguments_get_name(h, args, carg_++);
const char* sDx = grib_arguments_get_name(h, args, self->carg++); const char* sjScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* sDy = grib_arguments_get_name(h, args, self->carg++); const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++); const char* salternativeRowScanning = grib_arguments_get_name(h, args, 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 double d2r = acos(0.0) / 90.0; const double d2r = acos(0.0) / 90.0;
is_oblate = grib_is_earth_oblate(h); is_oblate = grib_is_earth_oblate(h);
@ -443,8 +376,8 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
if ((err = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS) if ((err = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS)
return err; return err;
if (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, iter->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; return GRIB_WRONG_GRID;
} }
if ((err = grib_get_double_internal(h, slatFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS) if ((err = grib_get_double_internal(h, slatFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS)
@ -474,14 +407,14 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
standardParallelInRadians = standardParallelInDegrees * d2r; standardParallelInRadians = standardParallelInDegrees * d2r;
if (is_oblate) { if (is_oblate) {
err = init_oblate(h, self, iter->nv, nx, ny, err = init_oblate(h, nv_, nx, ny,
Dx, Dy, earthMinorAxisInMetres, earthMajorAxisInMetres, Dx, Dy, earthMinorAxisInMetres, earthMajorAxisInMetres,
latFirstInRadians, lonFirstInRadians, latFirstInRadians, lonFirstInRadians,
centralLongitudeInRadians, standardParallelInRadians, centralLongitudeInRadians, standardParallelInRadians,
iScansNegatively, jScansPositively, jPointsAreConsecutive); iScansNegatively, jScansPositively, jPointsAreConsecutive);
} }
else { else {
err = init_sphere(h, self, iter->nv, nx, ny, err = init_sphere(h, nv_, nx, ny,
Dx, Dy, radius, Dx, Dy, radius,
latFirstInRadians, lonFirstInRadians, latFirstInRadians, lonFirstInRadians,
centralLongitudeInRadians, standardParallelInRadians, centralLongitudeInRadians, standardParallelInRadians,
@ -489,17 +422,19 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
} }
if (err) return err; if (err) return err;
iter->e = -1; e_ = -1;
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
static int destroy(grib_iterator* iter) int LambertAzimuthalEqualArea::destroy()
{ {
grib_iterator_lambert_azimuthal_equal_area* self = (grib_iterator_lambert_azimuthal_equal_area*)iter; const grib_context* c = h_->context;
const grib_context* c = iter->h->context;
grib_context_free(c, self->lats); grib_context_free(c, lats_);
grib_context_free(c, self->lons); grib_context_free(c, lons_);
return GRIB_SUCCESS;
return Gen::destroy();
} }
} // namespace eccodes::geo_iterator

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"
namespace eccodes::geo_iterator {
class LambertAzimuthalEqualArea : public Gen
{
public:
LambertAzimuthalEqualArea() :
Gen() { class_name_ = "lambert_azimuthal_equal_area"; }
Iterator* create() const override { return new LambertAzimuthalEqualArea(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int destroy() override;
private:
long Nj_;
int init_sphere(grib_handle*,
size_t, long, long,
double, double, double,
double, double,
double, double,
long, long, long);
int init_oblate(grib_handle*,
size_t, long, long,
double, double, double, double,
double, double,
double, double,
long, long, long);
};
} // namespace eccodes::geo_iterator

View File

@ -8,101 +8,36 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. * 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> #include <cmath>
/* eccodes::geo_iterator::LambertConformal _grib_iterator_lambert_conformal{};
This is used by make_class.pl eccodes::geo_iterator::Iterator* grib_iterator_lambert_conformal = &_grib_iterator_lambert_conformal;
START_CLASS_DEF namespace eccodes::geo_iterator {
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 */ #define ITER "Lambert conformal Geoiterator"
/*
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 */
#define ITER "Lambert conformal Geoiterator"
#define EPSILON 1.0e-10 #define EPSILON 1.0e-10
#ifndef M_PI #ifndef M_PI
#define M_PI 3.14159265358979323846 // Whole pie #define M_PI 3.14159265358979323846 // Whole pie
#endif #endif
#ifndef M_PI_2 #ifndef M_PI_2
#define M_PI_2 1.57079632679489661923 // Half a pie #define M_PI_2 1.57079632679489661923 // Half a pie
#endif #endif
#ifndef M_PI_4 #ifndef M_PI_4
#define M_PI_4 0.78539816339744830962 // Quarter of a pie #define M_PI_4 0.78539816339744830962 // Quarter of a pie
#endif #endif
#define RAD2DEG 57.29577951308232087684 // 180 over pi #define RAD2DEG 57.29577951308232087684 // 180 over pi
#define DEG2RAD 0.01745329251994329576 // pi over 180 #define DEG2RAD 0.01745329251994329576 // pi over 180
// Adjust longitude (in radians) to range -180 to 180 // Adjust longitude (in radians) to range -180 to 180
static double adjust_lon_radians(double lon) static double adjust_lon_radians(double lon)
{ {
if (lon > M_PI) lon -= 2 * M_PI; if (lon > M_PI) lon -= 2 * M_PI;
if (lon < -M_PI) lon += 2 * M_PI; if (lon < -M_PI) lon += 2 * M_PI;
return lon; return lon;
} }
@ -114,8 +49,8 @@ static double adjust_lon_radians(double lon)
// calculate phi on the left side. Substitute the calculated phi) into the right side, // calculate phi on the left side. Substitute the calculated phi) into the right side,
// calculate a new phi, etc., until phi does not change significantly from the preceding trial value of phi // calculate a new phi, etc., until phi does not change significantly from the preceding trial value of phi
static double compute_phi( static double compute_phi(
double eccent, // Spheroid eccentricity double eccent, // Spheroid eccentricity
double ts, // Constant value t double ts, // Constant value t
int* error) int* error)
{ {
double eccnth, phi, con, dphi, sinpi; double eccnth, phi, con, dphi, sinpi;
@ -145,9 +80,9 @@ static double compute_m(double eccent, double sinphi, double cosphi)
// Compute the constant small t for use in the forward computations // Compute the constant small t for use in the forward computations
static double compute_t( static double compute_t(
double eccent, // Eccentricity of the spheroid double eccent, // Eccentricity of the spheroid
double phi, // Latitude phi double phi, // Latitude phi
double sinphi) // Sine of the latitude double sinphi) // Sine of the latitude
{ {
double con = eccent * sinphi; double con = eccent * sinphi;
double com = 0.5 * eccent; double com = 0.5 * eccent;
@ -162,25 +97,25 @@ static double calculate_eccentricity(double minor, double major)
} }
static void xy2lonlat(double radius, double n, double f, double rho0_bare, double LoVInRadians, static void xy2lonlat(double radius, double n, double f, double rho0_bare, double LoVInRadians,
double x, double y, double x, double y,
double* lonDeg, double* latDeg) double* lonDeg, double* latDeg)
{ {
DEBUG_ASSERT(radius > 0); DEBUG_ASSERT(radius > 0);
DEBUG_ASSERT(n != 0.0); DEBUG_ASSERT(n != 0.0);
x /= radius; x /= radius;
y /= radius; y /= radius;
y = rho0_bare - y; y = rho0_bare - y;
double rho = hypot(x, y); double rho = hypot(x, y);
if (rho != 0.0) { if (rho != 0.0) {
if (n < 0.0) { if (n < 0.0) {
rho = -rho; rho = -rho;
x = -x; x = -x;
y = -y; y = -y;
} }
double latRadians = 2. * atan(pow(f / rho, 1.0/n)) - M_PI_2; double latRadians = 2. * atan(pow(f / rho, 1.0 / n)) - M_PI_2;
double lonRadians = atan2(x, y) / n; double lonRadians = atan2(x, y) / n;
*lonDeg = (lonRadians + LoVInRadians) * RAD2DEG; *lonDeg = (lonRadians + LoVInRadians) * RAD2DEG;
*latDeg = latRadians * RAD2DEG; *latDeg = latRadians * RAD2DEG;
} }
else { else {
*lonDeg = 0.0; *lonDeg = 0.0;
@ -188,29 +123,29 @@ static void xy2lonlat(double radius, double n, double f, double rho0_bare, doubl
} }
} }
static int init_sphere(const grib_handle* h, int LambertConformal::init_sphere(const grib_handle* h,
grib_iterator_lambert_conformal* self, size_t nv, long nx, long ny,
size_t nv, long nx, long ny, double LoVInDegrees,
double LoVInDegrees, double Dx, double Dy, double radius,
double Dx, double Dy, double radius, double latFirstInRadians, double lonFirstInRadians,
double latFirstInRadians, double lonFirstInRadians, double LoVInRadians, double Latin1InRadians, double Latin2InRadians,
double LoVInRadians, double Latin1InRadians, double Latin2InRadians, double LaDInRadians)
double LaDInRadians)
{ {
double n, x, y; double n, x, y;
if (fabs(Latin1InRadians - Latin2InRadians) < 1E-09) { if (fabs(Latin1InRadians - Latin2InRadians) < 1E-09) {
n = sin(Latin1InRadians); n = sin(Latin1InRadians);
} else { }
else {
n = log(cos(Latin1InRadians) / cos(Latin2InRadians)) / n = log(cos(Latin1InRadians) / cos(Latin2InRadians)) /
log(tan(M_PI_4 + Latin2InRadians / 2.0) / tan(M_PI_4 + Latin1InRadians / 2.0)); log(tan(M_PI_4 + Latin2InRadians / 2.0) / tan(M_PI_4 + Latin1InRadians / 2.0));
} }
double f = (cos(Latin1InRadians) * pow(tan(M_PI_4 + Latin1InRadians / 2.0), n)) / n; double f = (cos(Latin1InRadians) * pow(tan(M_PI_4 + Latin1InRadians / 2.0), n)) / n;
double rho = radius * f * pow(tan(M_PI_4 + latFirstInRadians / 2.0), -n); double rho = radius * f * pow(tan(M_PI_4 + latFirstInRadians / 2.0), -n);
double rho0_bare = f * pow(tan(M_PI_4 + LaDInRadians / 2.0), -n); double rho0_bare = f * pow(tan(M_PI_4 + LaDInRadians / 2.0), -n);
double rho0 = radius * rho0_bare; // scaled double rho0 = radius * rho0_bare; // scaled
double lonDiff = lonFirstInRadians - LoVInRadians; double lonDiff = lonFirstInRadians - LoVInRadians;
// Adjust longitude to range -180 to 180 // Adjust longitude to range -180 to 180
if (lonDiff > M_PI) if (lonDiff > M_PI)
@ -225,13 +160,13 @@ static int init_sphere(const grib_handle* h,
// Dy = jScansPositively == 1 ? Dy : -Dy; // Dy = jScansPositively == 1 ? Dy : -Dy;
// Allocate latitude and longitude arrays // Allocate latitude and longitude arrays
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lons_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
@ -241,11 +176,11 @@ static int init_sphere(const grib_handle* h,
y = y0 + j * Dy; y = y0 + j * Dy;
for (long i = 0; i < nx; i++) { for (long i = 0; i < nx; i++) {
const long index = i + j * nx; const long index = i + j * nx;
x = x0 + i * Dx; x = x0 + i * Dx;
xy2lonlat(radius, n, f, rho0_bare, LoVInRadians, x, y, &lonDeg, &latDeg); xy2lonlat(radius, n, f, rho0_bare, LoVInRadians, x, y, &lonDeg, &latDeg);
lonDeg = normalise_longitude_in_degrees(lonDeg); lonDeg = normalise_longitude_in_degrees(lonDeg);
self->lons[index] = lonDeg; lons_[index] = lonDeg;
self->lats[index] = latDeg; lats_[index] = latDeg;
} }
} }
@ -269,8 +204,8 @@ static int init_sphere(const grib_handle* h,
lonDeg = LoVInDegrees + (angle / n) * RAD2DEG; lonDeg = LoVInDegrees + (angle / n) * RAD2DEG;
latDeg = (2.0 * atan(pow(radius * f / rho, 1.0 / n)) - M_PI_2) * RAD2DEG; latDeg = (2.0 * atan(pow(radius * f / rho, 1.0 / n)) - M_PI_2) * RAD2DEG;
lonDeg = normalise_longitude_in_degrees(lonDeg); lonDeg = normalise_longitude_in_degrees(lonDeg);
self->lons[index] = lonDeg; lons_[index] = lonDeg;
self->lats[index] = latDeg; lats_[index] = latDeg;
} }
} }
#endif #endif
@ -278,32 +213,31 @@ static int init_sphere(const grib_handle* h,
} }
// Oblate spheroid // Oblate spheroid
static int init_oblate(const grib_handle* h, int LambertConformal::init_oblate(const grib_handle* h,
grib_iterator_lambert_conformal* self, size_t nv, long nx, long ny,
size_t nv, long nx, long ny, double LoVInDegrees,
double LoVInDegrees, double Dx, double Dy,
double Dx, double Dy, double earthMinorAxisInMetres, double earthMajorAxisInMetres,
double earthMinorAxisInMetres, double earthMajorAxisInMetres, double latFirstInRadians, double lonFirstInRadians,
double latFirstInRadians, double lonFirstInRadians, double LoVInRadians, double Latin1InRadians, double Latin2InRadians,
double LoVInRadians, double Latin1InRadians, double Latin2InRadians, double LaDInRadians)
double LaDInRadians)
{ {
int i, j, err = 0; int i, j, err = 0;
double x0, y0, x, y, latRad, lonRad, latDeg, lonDeg, sinphi, ts, rh1, theta; double x0, y0, x, y, latRad, lonRad, latDeg, lonDeg, sinphi, ts, rh1, theta;
double false_easting; // x offset in meters double false_easting; // x offset in meters
double false_northing; // y offset in meters double false_northing; // y offset in meters
double ns; // ratio of angle between meridian double ns; // ratio of angle between meridian
double F; // flattening of ellipsoid double F; // flattening of ellipsoid
double rh; // height above ellipsoid double rh; // height above ellipsoid
double sin_po; // sin value double sin_po; // sin value
double cos_po; // cos value double cos_po; // cos value
double con; // temporary variable double con; // temporary variable
double ms1; // small m 1 double ms1; // small m 1
double ms2; // small m 2 double ms2; // small m 2
double ts0; // small t 0 double ts0; // small t 0
double ts1; // small t 1 double ts1; // small t 1
double ts2; // small t 2 double ts2; // small t 2
double e = calculate_eccentricity(earthMinorAxisInMetres, earthMajorAxisInMetres); double e = calculate_eccentricity(earthMinorAxisInMetres, earthMajorAxisInMetres);
@ -322,10 +256,11 @@ static int init_oblate(const grib_handle* h,
if (fabs(Latin1InRadians - Latin2InRadians) > EPSILON) { if (fabs(Latin1InRadians - Latin2InRadians) > EPSILON) {
ns = log(ms1 / ms2) / log(ts1 / ts2); ns = log(ms1 / ms2) / log(ts1 / ts2);
} else { }
else {
ns = con; ns = con;
} }
F = ms1 / (ns * pow(ts1, ns)); F = ms1 / (ns * pow(ts1, ns));
rh = earthMajorAxisInMetres * F * pow(ts0, ns); rh = earthMajorAxisInMetres * F * pow(ts0, ns);
// Forward projection: convert lat,lon to x,y // Forward projection: convert lat,lon to x,y
@ -334,7 +269,8 @@ static int init_oblate(const grib_handle* h,
sinphi = sin(latFirstInRadians); sinphi = sin(latFirstInRadians);
ts = compute_t(e, latFirstInRadians, sinphi); ts = compute_t(e, latFirstInRadians, sinphi);
rh1 = earthMajorAxisInMetres * F * pow(ts, ns); rh1 = earthMajorAxisInMetres * F * pow(ts, ns);
} else { }
else {
con = latFirstInRadians * ns; con = latFirstInRadians * ns;
if (con <= 0) { if (con <= 0) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Point cannot be projected: latFirstInRadians=%g", ITER, latFirstInRadians); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Point cannot be projected: latFirstInRadians=%g", ITER, latFirstInRadians);
@ -349,13 +285,13 @@ static int init_oblate(const grib_handle* h,
y0 = -y0; y0 = -y0;
// Allocate latitude and longitude arrays // Allocate latitude and longitude arrays
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lons_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
@ -370,8 +306,8 @@ static int init_oblate(const grib_handle* h,
double _x, _y; double _x, _y;
x = i * Dx; x = i * Dx;
// Inverse projection to convert from x,y to lat,lon // Inverse projection to convert from x,y to lat,lon
_x = x - false_easting; _x = x - false_easting;
_y = rh - y + false_northing; _y = rh - y + false_northing;
rh1 = sqrt(_x * _x + _y * _y); rh1 = sqrt(_x * _x + _y * _y);
con = 1.0; con = 1.0;
if (ns <= 0) { if (ns <= 0) {
@ -388,54 +324,57 @@ static int init_oblate(const grib_handle* h,
if (err) { if (err) {
grib_context_log(h->context, GRIB_LOG_ERROR, grib_context_log(h->context, GRIB_LOG_ERROR,
"%s: Failed to compute the latitude angle, phi2, for the inverse", ITER); "%s: Failed to compute the latitude angle, phi2, for the inverse", ITER);
grib_context_free(h->context, self->lats); grib_context_free(h->context, lats_);
grib_context_free(h->context, self->lons); grib_context_free(h->context, lons_);
return err; return err;
} }
} else { }
else {
latRad = -M_PI_2; latRad = -M_PI_2;
} }
lonRad = adjust_lon_radians(theta / ns + LoVInRadians); lonRad = adjust_lon_radians(theta / ns + LoVInRadians);
if (i == 0 && j == 0) { if (i == 0 && j == 0) {
DEBUG_ASSERT(fabs(latFirstInRadians - latRad) <= EPSILON); DEBUG_ASSERT(fabs(latFirstInRadians - latRad) <= EPSILON);
} }
latDeg = latRad * RAD2DEG; // Convert to degrees latDeg = latRad * RAD2DEG; // Convert to degrees
lonDeg = normalise_longitude_in_degrees(lonRad * RAD2DEG); lonDeg = normalise_longitude_in_degrees(lonRad * RAD2DEG);
self->lons[index] = lonDeg; lons_[index] = lonDeg;
self->lats[index] = latDeg; lats_[index] = latDeg;
} }
} }
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args) int LambertConformal::init(grib_handle* h, grib_arguments* args)
{ {
int err = 0, is_oblate = 0; int err = GRIB_SUCCESS;
if ((err = Gen::init(h, args)) != GRIB_SUCCESS)
return err;
int is_oblate = 0;
long nx, ny, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning; long nx, ny, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning;
double LoVInDegrees, LaDInDegrees, Latin1InDegrees, Latin2InDegrees, latFirstInDegrees, double LoVInDegrees, LaDInDegrees, Latin1InDegrees, Latin2InDegrees, latFirstInDegrees,
lonFirstInDegrees, Dx, Dy, radius = 0; lonFirstInDegrees, Dx, Dy, radius = 0;
double latFirstInRadians, lonFirstInRadians, LoVInRadians, Latin1InRadians, Latin2InRadians, double latFirstInRadians, lonFirstInRadians, LoVInRadians, Latin1InRadians, Latin2InRadians,
LaDInRadians; LaDInRadians;
double earthMajorAxisInMetres=0, earthMinorAxisInMetres=0; double earthMajorAxisInMetres = 0, earthMinorAxisInMetres = 0;
grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)iter; const char* sradius = grib_arguments_get_name(h, args, carg_++);
const char* snx = grib_arguments_get_name(h, args, carg_++);
const char* sradius = grib_arguments_get_name(h, args, self->carg++); const char* sny = grib_arguments_get_name(h, args, carg_++);
const char* snx = grib_arguments_get_name(h, args, self->carg++); const char* sLoVInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sny = grib_arguments_get_name(h, args, self->carg++); const char* sLaDInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLoVInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* sLatin1InDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLaDInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* sLatin2InDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLatin1InDegrees = grib_arguments_get_name(h, args, self->carg++); const char* slatFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLatin2InDegrees = grib_arguments_get_name(h, args, self->carg++); const char* slonFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* slatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* slonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
// Dx and Dy are in Metres // Dx and Dy are in Metres
const char* sDx = 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, self->carg++); const char* sDy = grib_arguments_get_name(h, args, carg_++);
const char* siScansNegatively = 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, self->carg++); const char* sjScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++); const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* salternativeRowScanning = grib_arguments_get_name(h, args, self->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, snx, &nx)) != GRIB_SUCCESS) return err;
if ((err = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS) return err; if ((err = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS) return err;
@ -445,12 +384,13 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
if (is_oblate) { if (is_oblate) {
if ((err = grib_get_double_internal(h, "earthMinorAxisInMetres", &earthMinorAxisInMetres)) != GRIB_SUCCESS) return err; if ((err = grib_get_double_internal(h, "earthMinorAxisInMetres", &earthMinorAxisInMetres)) != GRIB_SUCCESS) return err;
if ((err = grib_get_double_internal(h, "earthMajorAxisInMetres", &earthMajorAxisInMetres)) != GRIB_SUCCESS) return err; if ((err = grib_get_double_internal(h, "earthMajorAxisInMetres", &earthMajorAxisInMetres)) != GRIB_SUCCESS) return err;
} else { }
else {
if ((err = grib_get_double_internal(h, sradius, &radius)) != GRIB_SUCCESS) return err; if ((err = grib_get_double_internal(h, sradius, &radius)) != GRIB_SUCCESS) return err;
} }
if (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, iter->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; return GRIB_WRONG_GRID;
} }
@ -497,14 +437,15 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
LoVInRadians = LoVInDegrees * DEG2RAD; LoVInRadians = LoVInDegrees * DEG2RAD;
if (is_oblate) { if (is_oblate) {
err = init_oblate(h, self, iter->nv, nx, ny, err = init_oblate(h, nv_, nx, ny,
LoVInDegrees, LoVInDegrees,
Dx, Dy, earthMinorAxisInMetres, earthMajorAxisInMetres, Dx, Dy, earthMinorAxisInMetres, earthMajorAxisInMetres,
latFirstInRadians, lonFirstInRadians, latFirstInRadians, lonFirstInRadians,
LoVInRadians, Latin1InRadians, Latin2InRadians, LoVInRadians, Latin1InRadians, Latin2InRadians,
LaDInRadians); LaDInRadians);
} else { }
err = init_sphere(h, self, iter->nv, nx, ny, else {
err = init_sphere(h, nv_, nx, ny,
LoVInDegrees, LoVInDegrees,
Dx, Dy, radius, Dx, Dy, radius,
latFirstInRadians, lonFirstInRadians, latFirstInRadians, lonFirstInRadians,
@ -512,37 +453,37 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
} }
if (err) return err; if (err) return err;
iter->e = -1; e_ = -1;
// Apply the scanning mode flags which may require data array to be transformed // 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, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning,
iter->nv, nx, ny); nv_, nx, ny);
return err; return err;
} }
static int next(grib_iterator* iter, double* lat, double* lon, double* val) int LambertConformal::next(double* lat, double* lon, double* val) const
{ {
grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)iter; if ((long)e_ >= (long)(nv_ - 1))
if ((long)iter->e >= (long)(iter->nv - 1))
return 0; return 0;
iter->e++; e_++;
*lat = self->lats[iter->e]; *lat = lats_[e_];
*lon = self->lons[iter->e]; *lon = lons_[e_];
if (val && iter->data) { if (val && data_) {
*val = iter->data[iter->e]; *val = data_[e_];
} }
return 1; return 1;
} }
static int destroy(grib_iterator* i) int LambertConformal::destroy()
{ {
grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)i; const grib_context* c = h_->context;
const grib_context* c = i->h->context;
grib_context_free(c, self->lats); grib_context_free(c, lats_);
grib_context_free(c, self->lons); grib_context_free(c, lons_);
return GRIB_SUCCESS;
return Gen::destroy();
} }
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,49 @@
/*
* (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"
namespace eccodes::geo_iterator {
class LambertConformal : public Gen
{
public:
LambertConformal() :
Gen() { class_name_ = "lambert_conformal"; }
Iterator* create() const override { return new LambertConformal(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int destroy() override;
private:
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);
};
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,157 @@
/*
* (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"
eccodes::geo_iterator::Latlon _grib_iterator_latlon{};
eccodes::geo_iterator::Iterator* grib_iterator_latlon = &_grib_iterator_latlon;
namespace eccodes::geo_iterator {
int Latlon::next(double* lat, double* lon, double* val) const
{
/* 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 = lats_[(long)floor(e_ / Ni_)];
ret_lon = lons_[(long)e_ % Ni_];
if (data_)
ret_val = data_[e_];
}
else {
/* Adjacent points in j (y) direction is consecutive */
ret_lon = lons_[(long)e_ / Nj_];
ret_lat = lats_[(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 Latlon::init(grib_handle* h, grib_arguments* args)
{
int err = 0;
if ((err = Regular::init(h, args)) != GRIB_SUCCESS)
return err;
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++) {
lats_[lai] = lat1;
lat1 -= jdir;
}
/* ECC-1406: Due to rounding, errors can accumulate.
* So we ensure the last latitude is latitudeOfLastGridPointInDegrees
*/
lats_[Nj_ - 1] = lat2;
e_ = -1;
return err;
}
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,28 @@
/*
* (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"
namespace eccodes::geo_iterator {
class Latlon : public Regular
{
public:
Latlon() :
Regular() { class_name_ = "latlon"; }
Iterator* create() const override { return new Latlon(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
};
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,141 @@
/*
* (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"
eccodes::geo_iterator::LatlonReduced _grib_iterator_latlon_reduced{};
eccodes::geo_iterator::Iterator* grib_iterator_latlon_reduced = &_grib_iterator_latlon_reduced;
namespace eccodes::geo_iterator {
int LatlonReduced::next(double* lat, double* lon, double* val) const
{
if ((long)e_ >= (long)(nv_ - 1))
return 0;
e_++;
*lat = lats_[e_];
*lon = lons_[e_];
if (val && data_) {
*val = data_[e_];
}
return 1;
}
int LatlonReduced::init(grib_handle* h, grib_arguments* args)
{
int ret = GRIB_SUCCESS;
if ((ret = Gen::init(h, args)) != GRIB_SUCCESS)
return ret;
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);
lats_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
lons_ = (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++) {
lats_[k] = laf;
lons_[k] = tlof;
tlof += idirinc;
k++;
}
laf += jdirinc;
}
e_ = -1;
grib_context_free(h->context, pl);
return ret;
}
int LatlonReduced::destroy()
{
const grib_context* c = h_->context;
grib_context_free(c, lats_);
grib_context_free(c, lons_);
return Gen::destroy();
}
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,29 @@
/*
* (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"
namespace eccodes::geo_iterator {
class LatlonReduced : public Gen
{
public:
LatlonReduced() :
Gen() { class_name_ = "latlon_reduced"; }
Iterator* create() const override { return new LatlonReduced(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int destroy() override;
};
} // namespace eccodes::geo_iterator

View File

@ -8,92 +8,26 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. * 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>
/* eccodes::geo_iterator::Mercator _grib_iterator_mercator{};
This is used by make_class.pl eccodes::geo_iterator::Iterator* grib_iterator_mercator = &_grib_iterator_mercator;
START_CLASS_DEF namespace eccodes::geo_iterator {
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 */ #define ITER "Mercator Geoiterator"
/*
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 */
#define ITER "Mercator Geoiterator"
#define EPSILON 1.0e-10 #define EPSILON 1.0e-10
#ifndef M_PI #ifndef M_PI
#define M_PI 3.14159265358979323846 /* Whole pie */ #define M_PI 3.14159265358979323846 /* Whole pie */
#endif #endif
#ifndef M_PI_2 #ifndef M_PI_2
#define M_PI_2 1.57079632679489661923 /* Half a pie */ #define M_PI_2 1.57079632679489661923 /* Half a pie */
#endif #endif
#ifndef M_PI_4 #ifndef M_PI_4
#define M_PI_4 0.78539816339744830962 /* Quarter of a pie */ #define M_PI_4 0.78539816339744830962 /* Quarter of a pie */
#endif #endif
#define RAD2DEG 57.29577951308232087684 /* 180 over pi */ #define RAD2DEG 57.29577951308232087684 /* 180 over pi */
@ -102,7 +36,7 @@ static void init_class(grib_iterator_class* c)
/* Adjust longitude (in radians) to range -180 to 180 */ /* Adjust longitude (in radians) to range -180 to 180 */
static double adjust_lon_radians(double lon) static double adjust_lon_radians(double lon)
{ {
if (lon > M_PI) lon -= 2 * M_PI; if (lon > M_PI) lon -= 2 * M_PI;
if (lon < -M_PI) lon += 2 * M_PI; if (lon < -M_PI) lon += 2 * M_PI;
return lon; return lon;
} }
@ -148,14 +82,13 @@ static double compute_t(
return (tan(0.5 * (M_PI_2 - phi)) / con); return (tan(0.5 * (M_PI_2 - phi)) / con);
} }
static int init_mercator(grib_handle* h, int Mercator::init_mercator(grib_handle* h,
grib_iterator_mercator* self, size_t nv, long nx, long ny,
size_t nv, long nx, long ny, double DiInMetres, double DjInMetres,
double DiInMetres, double DjInMetres, double earthMinorAxisInMetres, double earthMajorAxisInMetres,
double earthMinorAxisInMetres, double earthMajorAxisInMetres, double latFirstInRadians, double lonFirstInRadians,
double latFirstInRadians, double lonFirstInRadians, double latLastInRadians, double lonLastInRadians,
double latLastInRadians, double lonLastInRadians, double LaDInRadians, double orientationInRadians)
double LaDInRadians, double orientationInRadians)
{ {
int i, j, err = 0; int i, j, err = 0;
double x0, y0, x, y, latRad, lonRad, latDeg, lonDeg, sinphi, ts; double x0, y0, x, y, latRad, lonRad, latDeg, lonDeg, sinphi, ts;
@ -184,13 +117,13 @@ static int init_mercator(grib_handle* h,
y0 = -y0; y0 = -y0;
/* Allocate latitude and longitude arrays */ /* Allocate latitude and longitude arrays */
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lats_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double)); lons_ = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!lons_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
@ -211,26 +144,29 @@ static int init_mercator(grib_handle* h,
latRad = compute_phi(e, ts, &err); latRad = compute_phi(e, ts, &err);
if (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_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, lats_);
grib_context_free(h->context, self->lons); grib_context_free(h->context, lons_);
return err; return err;
} }
lonRad = adjust_lon_radians(orientationInRadians + _x / (earthMajorAxisInMetres * m1)); lonRad = adjust_lon_radians(orientationInRadians + _x / (earthMajorAxisInMetres * m1));
if (i == 0 && j == 0) { if (i == 0 && j == 0) {
DEBUG_ASSERT(fabs(latFirstInRadians - latRad) <= EPSILON); DEBUG_ASSERT(fabs(latFirstInRadians - latRad) <= EPSILON);
} }
latDeg = latRad * RAD2DEG; /* Convert to degrees */ latDeg = latRad * RAD2DEG; /* Convert to degrees */
lonDeg = normalise_longitude_in_degrees(lonRad * RAD2DEG); lonDeg = normalise_longitude_in_degrees(lonRad * RAD2DEG);
self->lons[index] = lonDeg; lons_[index] = lonDeg;
self->lats[index] = latDeg; lats_[index] = latDeg;
} }
} }
return GRIB_SUCCESS; return GRIB_SUCCESS;
} }
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args) int Mercator::init(grib_handle* h, grib_arguments* args)
{ {
int err = 0; int err = GRIB_SUCCESS;
if ((err = Gen::init(h, args)) != GRIB_SUCCESS)
return err;
long ni, nj, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning; long ni, nj, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning;
double latFirstInDegrees, lonFirstInDegrees, LaDInDegrees; double latFirstInDegrees, lonFirstInDegrees, LaDInDegrees;
double latLastInDegrees, lonLastInDegrees, orientationInDegrees, DiInMetres, DjInMetres, radius = 0; double latLastInDegrees, lonLastInDegrees, orientationInDegrees, DiInMetres, DjInMetres, radius = 0;
@ -238,24 +174,22 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
LaDInRadians, orientationInRadians; LaDInRadians, orientationInRadians;
double earthMajorAxisInMetres = 0, earthMinorAxisInMetres = 0; double earthMajorAxisInMetres = 0, earthMinorAxisInMetres = 0;
grib_iterator_mercator* self = (grib_iterator_mercator*)iter; const char* sRadius = grib_arguments_get_name(h, args, carg_++);
const char* sNi = grib_arguments_get_name(h, args, carg_++);
const char* sRadius = grib_arguments_get_name(h, args, self->carg++); const char* sNj = grib_arguments_get_name(h, args, carg_++);
const char* sNi = grib_arguments_get_name(h, args, self->carg++); const char* sLatFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sNj = grib_arguments_get_name(h, args, self->carg++); const char* sLonFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* sLaDInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* sLatLastInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLaDInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* sLonLastInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLatLastInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* sOrientationInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* sLonLastInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sOrientationInDegrees = grib_arguments_get_name(h, args, self->carg++);
/* Dx and Dy are in Metres */ /* Dx and Dy are in Metres */
const char* sDi = 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, self->carg++); const char* sDj = grib_arguments_get_name(h, args, carg_++);
const char* siScansNegatively = 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, self->carg++); const char* sjScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++); const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* sAlternativeRowScanning = grib_arguments_get_name(h, args, self->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, sNi, &ni)) != GRIB_SUCCESS) return err;
if ((err = grib_get_long_internal(h, sNj, &nj)) != GRIB_SUCCESS) return err; if ((err = grib_get_long_internal(h, sNj, &nj)) != GRIB_SUCCESS) return err;
@ -269,8 +203,8 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
earthMinorAxisInMetres = earthMajorAxisInMetres = radius; earthMinorAxisInMetres = earthMajorAxisInMetres = radius;
} }
if (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, iter->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; return GRIB_WRONG_GRID;
} }
@ -307,44 +241,44 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
LaDInRadians = LaDInDegrees * DEG2RAD; LaDInRadians = LaDInDegrees * DEG2RAD;
orientationInRadians = orientationInDegrees * DEG2RAD; orientationInRadians = orientationInDegrees * DEG2RAD;
err = init_mercator(h, self, iter->nv, ni, nj, err = init_mercator(h, nv_, ni, nj,
DiInMetres, DjInMetres, earthMinorAxisInMetres, earthMajorAxisInMetres, DiInMetres, DjInMetres, earthMinorAxisInMetres, earthMajorAxisInMetres,
latFirstInRadians, lonFirstInRadians, latFirstInRadians, lonFirstInRadians,
latLastInRadians, lonLastInRadians, latLastInRadians, lonLastInRadians,
LaDInRadians, orientationInRadians); LaDInRadians, orientationInRadians);
if (err) return err; if (err) return err;
iter->e = -1; e_ = -1;
/* Apply the scanning mode flags which may require data array to be transformed */ /* 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, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning,
iter->nv, ni, nj); nv_, ni, nj);
return err; return err;
} }
static int next(grib_iterator* iter, double* lat, double* lon, double* val) int Mercator::next(double* lat, double* lon, double* val) const
{ {
grib_iterator_mercator* self = (grib_iterator_mercator*)iter; if ((long)e_ >= (long)(nv_ - 1))
if ((long)iter->e >= (long)(iter->nv - 1))
return 0; return 0;
iter->e++; e_++;
*lat = self->lats[iter->e]; *lat = lats_[e_];
*lon = self->lons[iter->e]; *lon = lons_[e_];
if (val && iter->data) { if (val && data_) {
*val = iter->data[iter->e]; *val = data_[e_];
} }
return 1; return 1;
} }
static int destroy(grib_iterator* iter) int Mercator::destroy()
{ {
grib_iterator_mercator* self = (grib_iterator_mercator*)iter; const grib_context* c = h_->context;
const grib_context* c = iter->h->context;
grib_context_free(c, self->lats); grib_context_free(c, lats_);
grib_context_free(c, self->lons); grib_context_free(c, lons_);
return GRIB_SUCCESS;
return Gen::destroy();
} }
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,39 @@
/*
* (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"
namespace eccodes::geo_iterator {
class Mercator : public Gen
{
public:
Mercator() { class_name_ = "mercator"; }
Iterator* create() const override { return new Mercator(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int destroy() override;
private:
long Nj_;
int init_mercator(grib_handle*,
size_t, long, long,
double, double,
double, double,
double, double,
double, double,
double, double);
};
} // namespace eccodes::geo_iterator

View File

@ -8,93 +8,25 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. * 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>
/* eccodes::geo_iterator::PolarStereographic _grib_iterator_polar_stereographic{};
This is used by make_class.pl eccodes::geo_iterator::Iterator* grib_iterator_polar_stereographic = &_grib_iterator_polar_stereographic;
START_CLASS_DEF namespace eccodes::geo_iterator {
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 */
#define ITER "Polar stereographic Geoiterator" #define ITER "Polar stereographic Geoiterator"
static int next(grib_iterator* iter, double* lat, double* lon, double* val) int PolarStereographic::next(double* lat, double* lon, double* val) const
{ {
grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)iter; if ((long)e_ >= (long)(nv_ - 1))
if ((long)iter->e >= (long)(iter->nv - 1))
return 0; return 0;
iter->e++; e_++;
*lat = self->lats[iter->e]; *lat = lats_[e_];
*lon = self->lons[iter->e]; *lon = lons_[e_];
if (val && iter->data) { if (val && data_) {
*val = iter->data[iter->e]; *val = data_[e_];
} }
return 1; return 1;
} }
@ -112,14 +44,17 @@ typedef struct proj_data_t
double false_easting; /* x offset in meters */ double false_easting; /* x offset in meters */
} proj_data_t; } proj_data_t;
#define RAD2DEG 57.29577951308232087684 /* 180 over pi */ #define RAD2DEG 57.29577951308232087684 /* 180 over pi */
#define DEG2RAD 0.01745329251994329576 /* pi over 180 */ #define DEG2RAD 0.01745329251994329576 /* pi over 180 */
#define PI_OVER_2 1.5707963267948966 /* half pi */ #define PI_OVER_2 1.5707963267948966 /* half pi */
#define EPSILON 1.0e-10 #define EPSILON 1.0e-10
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args) int PolarStereographic::init(grib_handle* h, grib_arguments* args)
{ {
int ret = 0; int ret = GRIB_SUCCESS;
if ((ret = Gen::init(h, args)) != GRIB_SUCCESS)
return ret;
double *lats, *lons; /* arrays for latitudes and longitudes */ double *lats, *lons; /* arrays for latitudes and longitudes */
double lonFirstInDegrees, latFirstInDegrees, radius; double lonFirstInDegrees, latFirstInDegrees, radius;
double x, y, Dx, Dy; double x, y, Dx, Dy;
@ -132,25 +67,27 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
double ts; /* value of small t */ double ts; /* value of small t */
double height; /* height above ellipsoid */ double height; /* height above ellipsoid */
double x0, y0, lonFirst, latFirst; double x0, y0, lonFirst, latFirst;
proj_data_t fwd_proj_data = {0,}; proj_data_t fwd_proj_data = {
proj_data_t inv_proj_data = {0,}; 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, carg_++);
const char* s_nx = grib_arguments_get_name(h, args, carg_++);
const char* s_radius = grib_arguments_get_name(h, args, self->carg++); const char* s_ny = grib_arguments_get_name(h, args, carg_++);
const char* s_nx = grib_arguments_get_name(h, args, self->carg++); const char* s_latFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* s_ny = grib_arguments_get_name(h, args, self->carg++); const char* s_lonFirstInDegrees = grib_arguments_get_name(h, args, carg_++);
const char* s_latFirstInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* s_southPoleOnPlane = grib_arguments_get_name(h, args, carg_++);
const char* s_lonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++); const char* s_centralLongitude = grib_arguments_get_name(h, args, carg_++);
const char* s_southPoleOnPlane = grib_arguments_get_name(h, args, self->carg++); const char* s_centralLatitude = grib_arguments_get_name(h, args, carg_++);
const char* s_centralLongitude = grib_arguments_get_name(h, args, self->carg++); const char* s_Dx = grib_arguments_get_name(h, args, carg_++);
const char* s_centralLatitude = grib_arguments_get_name(h, args, self->carg++); const char* s_Dy = grib_arguments_get_name(h, args, carg_++);
const char* s_Dx = grib_arguments_get_name(h, args, self->carg++); const char* s_iScansNegatively = grib_arguments_get_name(h, args, carg_++);
const char* s_Dy = grib_arguments_get_name(h, args, self->carg++); const char* s_jScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* s_iScansNegatively = grib_arguments_get_name(h, args, self->carg++); const char* s_jPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* s_jScansPositively = grib_arguments_get_name(h, args, self->carg++); const char* s_alternativeRowScanning = grib_arguments_get_name(h, args, 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++);
if (grib_is_earth_oblate(h)) { if (grib_is_earth_oblate(h)) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Only supported for spherical earth.", ITER); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Only supported for spherical earth.", ITER);
@ -164,8 +101,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) if ((ret = grib_get_long_internal(h, s_ny, &ny)) != GRIB_SUCCESS)
return ret; return ret;
if (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, iter->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; return GRIB_WRONG_GRID;
} }
if ((ret = grib_get_double_internal(h, s_latFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS) if ((ret = grib_get_double_internal(h, s_latFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS)
@ -243,18 +180,18 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
inv_proj_data.mcs = cos(con1); inv_proj_data.mcs = cos(con1);
inv_proj_data.tcs = tan(0.5 * (PI_OVER_2 - con1)); inv_proj_data.tcs = tan(0.5 * (PI_OVER_2 - con1));
} }
self->lats = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double)); lats_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (!self->lats) { if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, iter->nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv_ * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
self->lons = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double)); lons_ = (double*)grib_context_malloc(h->context, nv_ * sizeof(double));
if (!self->lats) { if (!lons_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, iter->nv * sizeof(double)); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, nv_ * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
lats = self->lats; lats = lats_;
lons = self->lons; lons = lons_;
/* These will be processed later in transform_iterator_data() */ /* These will be processed later in transform_iterator_data() */
/* Dx = iScansNegatively == 0 ? Dx : -Dx; */ /* Dx = iScansNegatively == 0 ? Dx : -Dx; */
/* Dy = jScansPositively == 1 ? Dy : -Dy; */ /* Dy = jScansPositively == 1 ? Dy : -Dy; */
@ -294,81 +231,83 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
y += Dy; y += Dy;
} }
// /*standardParallel = (southPoleOnPlane == 1) ? -90 : +90;*/ // /*standardParallel = (southPoleOnPlane == 1) ? -90 : +90;*/
// if (jPointsAreConsecutive) // if (jPointsAreConsecutive)
// { // {
// x=xFirst; // x=xFirst;
// for (i=0;i<nx;i++) { // for (i=0;i<nx;i++) {
// y=yFirst; // y=yFirst;
// for (j=0;j<ny;j++) { // for (j=0;j<ny;j++) {
// rho=sqrt(x*x+y*y); // rho=sqrt(x*x+y*y);
// if (rho == 0) { // if (rho == 0) {
// /* indeterminate case */ // /* indeterminate case */
// *lats = standardParallel; // *lats = standardParallel;
// *lons = centralLongitude; // *lons = centralLongitude;
// } // }
// else { // else {
// c=2*atan2(rho,(2.0*radius)); // c=2*atan2(rho,(2.0*radius));
// cosc=cos(c); // cosc=cos(c);
// sinc=sin(c); // sinc=sin(c);
// *lats = asin( cosc*sinphi1 + y*sinc*cosphi1/rho ) * RAD2DEG; // *lats = asin( cosc*sinphi1 + y*sinc*cosphi1/rho ) * RAD2DEG;
// *lons = (lambda0+atan2(x*sinc, rho*cosphi1*cosc - y*sinphi1*sinc)) * RAD2DEG; // *lons = (lambda0+atan2(x*sinc, rho*cosphi1*cosc - y*sinphi1*sinc)) * RAD2DEG;
// } // }
// while (*lons<0) *lons += 360; // while (*lons<0) *lons += 360;
// while (*lons>360) *lons -= 360; // while (*lons>360) *lons -= 360;
// lons++; // lons++;
// lats++; // lats++;
// y+=Dy; // y+=Dy;
// } // }
// x+=Dx; // x+=Dx;
// } // }
// } // }
// else // else
// { // {
// y=yFirst; // y=yFirst;
// for (j=0;j<ny;j++) { // for (j=0;j<ny;j++) {
// x=xFirst; // x=xFirst;
// for (i=0;i<nx;i++) { // for (i=0;i<nx;i++) {
// /* int index =i+j*nx; */ // /* int index =i+j*nx; */
// rho=sqrt(x*x+y*y); // rho=sqrt(x*x+y*y);
// if (rho == 0) { // if (rho == 0) {
// /* indeterminate case */ // /* indeterminate case */
// *lats = standardParallel; // *lats = standardParallel;
// *lons = centralLongitude; // *lons = centralLongitude;
// } // }
// else { // else {
// c=2*atan2(rho,(2.0*radius)); // c=2*atan2(rho,(2.0*radius));
// cosc=cos(c); // cosc=cos(c);
// sinc=sin(c); // sinc=sin(c);
// *lats = asin( cosc*sinphi1 + y*sinc*cosphi1/rho ) * RAD2DEG; // *lats = asin( cosc*sinphi1 + y*sinc*cosphi1/rho ) * RAD2DEG;
// *lons = (lambda0+atan2(x*sinc, rho*cosphi1*cosc - y*sinphi1*sinc)) * RAD2DEG; // *lons = (lambda0+atan2(x*sinc, rho*cosphi1*cosc - y*sinphi1*sinc)) * RAD2DEG;
// } // }
// while (*lons<0) *lons += 360; // while (*lons<0) *lons += 360;
// while (*lons>360) *lons -= 360; // while (*lons>360) *lons -= 360;
// lons++; // lons++;
// lats++; // lats++;
// x+=Dx; // x+=Dx;
// } // }
// y+=Dy; // y+=Dy;
// } // }
// } // }
iter->e = -1; e_ = -1;
/* Apply the scanning mode flags which may require data array to be transformed */ /* 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, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning,
iter->nv, nx, ny); nv_, nx, ny);
return ret; return ret;
} }
static int destroy(grib_iterator* i) int PolarStereographic::destroy()
{ {
grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)i; const grib_context* c = h_->context;
const grib_context* c = i->h->context;
grib_context_free(c, self->lats); grib_context_free(c, lats_);
grib_context_free(c, self->lons); grib_context_free(c, lons_);
return GRIB_SUCCESS;
return Gen::destroy();
} }
} // namespace eccodes::geo_iterator

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"
namespace eccodes::geo_iterator {
class PolarStereographic : public Gen
{
public:
PolarStereographic() :
Gen() { class_name_ = "polar_stereographic"; }
Iterator* create() const override { return new PolarStereographic(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int destroy() override;
private:
long Nj_;
};
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,162 @@
/*
* (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"
eccodes::geo_iterator::Regular _grib_iterator_regular{};
eccodes::geo_iterator::Iterator* grib_iterator_regular = &_grib_iterator_regular;
namespace eccodes::geo_iterator {
#define ITER "Regular grid Geoiterator"
int Regular::next(double* lat, double* lon, double* val) const
{
if ((long)e_ >= (long)(nv_ - 1))
return 0;
e_++;
*lat = lats_[(long)floor(e_ / Ni_)];
*lon = lons_[(long)e_ % Ni_];
if (val && data_) {
*val = data_[e_];
}
return 1;
}
int Regular::previous(double* lat, double* lon, double* val) const
{
if (e_ < 0)
return 0;
*lat = lats_[(long)floor(e_ / Ni_)];
*lon = lons_[e_ % Ni_];
if (val && data_) {
*val = data_[e_];
}
e_--;
return 1;
}
int Regular::destroy()
{
const grib_context* c = h_->context;
grib_context_free(c, lats_);
grib_context_free(c, lons_);
return Gen::destroy();
}
int Regular::init(grib_handle* h, grib_arguments* args)
{
int ret = GRIB_SUCCESS;
if ((ret = Gen::init(h, args)) != GRIB_SUCCESS)
return ret;
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;
lats_ = (double*)grib_context_malloc(h->context, Nj * sizeof(double));
lons_ = (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++) {
lons_[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);
}
lons_[Ni - 1] = lon2;
return ret;
}
} // namespace eccodes::geo_iterator

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"
namespace eccodes::geo_iterator {
class Regular : public Gen
{
public:
Regular() { class_name_ = "regular"; }
Iterator* create() const override { return new Regular(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int previous(double*, double*, double*) const override;
int destroy() override;
protected:
long Ni_;
long Nj_;
long iScansNegatively_;
long isRotated_;
double angleOfRotation_;
double southPoleLat_;
double southPoleLon_;
long jPointsAreConsecutive_;
long disableUnrotate_;
};
} // namespace eccodes::geo_iterator

View File

@ -8,93 +8,25 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. * 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>
/* eccodes::geo_iterator::SpaceView _grib_iterator_space_view{};
This is used by make_class.pl eccodes::geo_iterator::Iterator* grib_iterator_space_view = &_grib_iterator_space_view;
START_CLASS_DEF namespace eccodes::geo_iterator {
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 */
#define ITER "Space view Geoiterator" #define ITER "Space view Geoiterator"
static int next(grib_iterator* iter, double* lat, double* lon, double* val) int SpaceView::next(double* lat, double* lon, double* val) const
{ {
grib_iterator_space_view* self = (grib_iterator_space_view*)iter; if ((long)e_ >= (long)(nv_ - 1))
if ((long)iter->e >= (long)(iter->nv - 1))
return 0; return 0;
iter->e++; e_++;
*lat = self->lats[iter->e]; *lat = lats_[e_];
*lon = self->lons[iter->e]; *lon = lons_[e_];
if (val && iter->data) { if (val && data_) {
*val = iter->data[iter->e]; *val = data_[e_];
} }
return 1; return 1;
} }
@ -142,12 +74,15 @@ static int next(grib_iterator* iter, double* lat, double* lon, double* val)
#define RAD2DEG 57.29577951308232087684 /* 180 over pi */ #define RAD2DEG 57.29577951308232087684 /* 180 over pi */
#define DEG2RAD 0.01745329251994329576 /* pi over 180 */ #define DEG2RAD 0.01745329251994329576 /* pi over 180 */
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args) int SpaceView::init(grib_handle* h, grib_arguments* args)
{ {
/* REFERENCE:
* LRIT/HRIT Global Specification (CGMS 03, Issue 2.6, 12.08.1999)
*/
int ret = GRIB_SUCCESS; int ret = GRIB_SUCCESS;
if ((ret = Gen::init(h, args)) != GRIB_SUCCESS)
return ret;
/* REFERENCE:
* LRIT/HRIT Global Specification (CGMS 03, Issue 2.6, 12.08.1999)
*/
double *lats, *lons; /* arrays of latitudes and longitudes */ double *lats, *lons; /* arrays of latitudes and longitudes */
double latOfSubSatellitePointInDegrees, lonOfSubSatellitePointInDegrees; double latOfSubSatellitePointInDegrees, lonOfSubSatellitePointInDegrees;
double orientationInDegrees, nrInRadiusOfEarth; double orientationInDegrees, nrInRadiusOfEarth;
@ -163,31 +98,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; double factor_1, factor_2, tmp1, Sd, Sn, Sxy, S1, S2, S3;
int x0, y0, ix, iy; int x0, y0, ix, iy;
double *s_x, *c_x; /* arrays storing sin and cos values */ 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* siScansNegatively = grib_arguments_get_name(h, args, carg_++);
const char* sEarthIsOblate = grib_arguments_get_name(h, args, self->carg++); const char* sjScansPositively = grib_arguments_get_name(h, args, carg_++);
const char* sMajorAxisInMetres = grib_arguments_get_name(h, args, self->carg++); const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, carg_++);
const char* sMinorAxisInMetres = grib_arguments_get_name(h, args, self->carg++); const char* sAlternativeRowScanning = grib_arguments_get_name(h, args, 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++);
if ((ret = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS) if ((ret = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS)
return ret; return ret;
@ -207,8 +140,8 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
return ret; return ret;
} }
if (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, iter->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; return GRIB_WRONG_GRID;
} }
if ((ret = grib_get_double_internal(h, sLatOfSubSatellitePointInDegrees, &latOfSubSatellitePointInDegrees)) != GRIB_SUCCESS) if ((ret = grib_get_double_internal(h, sLatOfSubSatellitePointInDegrees, &latOfSubSatellitePointInDegrees)) != GRIB_SUCCESS)
@ -288,18 +221,18 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
rx = angular_size / dx; rx = angular_size / dx;
ry = (r_pol / r_eq) * angular_size / dy; ry = (r_pol / r_eq) * angular_size / dy;
self->lats = (double*)grib_context_malloc(h->context, array_size); lats_ = (double*)grib_context_malloc(h->context, array_size);
if (!self->lats) { if (!lats_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, array_size); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, array_size);
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
self->lons = (double*)grib_context_malloc(h->context, array_size); lons_ = (double*)grib_context_malloc(h->context, array_size);
if (!self->lats) { if (!lons_) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, array_size); grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Error allocating %zu bytes", ITER, array_size);
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
lats = self->lats; lats = lats_;
lons = self->lons; lons = lons_;
if (!iScansNegatively) { if (!iScansNegatively) {
xp = xp - x0; xp = xp - x0;
@ -373,17 +306,19 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
} }
grib_context_free(h->context, s_x); grib_context_free(h->context, s_x);
grib_context_free(h->context, c_x); grib_context_free(h->context, c_x);
iter->e = -1; e_ = -1;
return ret; return ret;
} }
static int destroy(grib_iterator* iter) int SpaceView::destroy()
{ {
grib_iterator_space_view* self = (grib_iterator_space_view*)iter; const grib_context* c = h_->context;
const grib_context* c = iter->h->context;
grib_context_free(c, self->lats); grib_context_free(c, lats_);
grib_context_free(c, self->lons); grib_context_free(c, lons_);
return GRIB_SUCCESS;
return Gen::destroy();
} }
} // namespace eccodes::geo_iterator

View File

@ -0,0 +1,31 @@
/*
* (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"
namespace eccodes::geo_iterator {
class SpaceView : public Gen
{
public:
SpaceView() { class_name_ = "space_view"; }
Iterator* create() const override { return new SpaceView(); }
int init(grib_handle*, grib_arguments*) override;
int next(double*, double*, double*) const override;
int destroy() override;
private:
long Nj_;
};
} // namespace eccodes::geo_iterator

View File

@ -248,7 +248,14 @@ typedef struct grib_codetable grib_codetable;
typedef struct grib_smart_table grib_smart_table; typedef struct grib_smart_table grib_smart_table;
class grib_accessor; class grib_accessor;
typedef struct grib_iterator_class grib_iterator_class; namespace eccodes::geo_iterator {
class Iterator;
}
typedef struct grib_iterator {
eccodes::geo_iterator::Iterator* iterator;
} grib_iterator;
typedef struct grib_nearest_class grib_nearest_class; typedef struct grib_nearest_class grib_nearest_class;
typedef struct grib_dumper grib_dumper; typedef struct grib_dumper grib_dumper;
typedef struct grib_dumper_class grib_dumper_class; typedef struct grib_dumper_class grib_dumper_class;
@ -267,14 +274,6 @@ typedef int (*nearest_find_proc)(grib_nearest* nearest, grib_handle* h,
double* distances, int* indexes, size_t* len); double* distances, int* indexes, size_t* len);
typedef int (*nearest_destroy_proc)(grib_nearest* nearest); 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_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); typedef int (*grib_unpack_proc)(grib_handle* h, const void* in, size_t inlen, double* out, size_t* outlen);
@ -459,20 +458,6 @@ struct grib_section
size_t padding; 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 struct grib_nearest_class
{ {
@ -529,17 +514,6 @@ struct grib_dumper_class
dumper_footer_proc footer; 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 struct grib_nearest
{ {
grib_arguments* args; /** args of iterator */ grib_arguments* args; /** args of iterator */
@ -1280,6 +1254,7 @@ typedef struct j2k_encode_helper
} }
#include "accessor/grib_accessor.h" #include "accessor/grib_accessor.h"
#include "accessor/grib_accessors_list.h" #include "accessor/grib_accessors_list.h"
#include "geo_iterator/grib_iterator.h"
#endif #endif
#endif #endif

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;
}

View File

@ -14,22 +14,53 @@
***************************************************************************/ ***************************************************************************/
#include "grib_api_internal.h" #include "grib_api_internal.h"
#include "geo_iterator/grib_iterator.h"
#include "accessor/grib_accessor_class_iterator.h"
/* This file is generated by ./make_class.pl */ /* This file is generated by ./make_class.pl */
#include "grib_iterator_class.h" #include "grib_iterator_class.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_grib_accessor_class_c)
{
if (once == 0) {
omp_init_nest_lock(&mutex);
once = 1;
}
}
}
#endif
struct table_entry struct table_entry
{ {
const char* type; const char* type;
grib_iterator_class** cclass; eccodes::geo_iterator::Iterator** iterator;
}; };
static const struct table_entry table[] = { static const struct table_entry table[] = {
/* This file is generated by ./make_class.pl */ /* This file is generated by ./make_class.pl */
#include "grib_iterator_factory.h" #include "grib_iterator_factory.h"
}; };
grib_iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsigned long flags, int* error)
eccodes::geo_iterator::Iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsigned long flags, int* error)
{ {
size_t i = 0, num_table_entries = 0; size_t i = 0, num_table_entries = 0;
const char* type = (char*)grib_arguments_get_name(h, args, 0); const char* type = (char*)grib_arguments_get_name(h, args, 0);
@ -38,16 +69,20 @@ grib_iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsig
num_table_entries = sizeof(table) / sizeof(table[0]); num_table_entries = sizeof(table) / sizeof(table[0]);
for (i = 0; i < num_table_entries; i++) { for (i = 0; i < num_table_entries; i++) {
if (strcmp(type, table[i].type) == 0) { if (strcmp(type, table[i].type) == 0) {
grib_iterator_class* c = *(table[i].cclass); eccodes::geo_iterator::Iterator* builder = *(table[i].iterator);
grib_iterator* it = (grib_iterator*)grib_context_malloc_clear(h->context, c->size); eccodes::geo_iterator::Iterator* it = builder->create();
it->cclass = c; it->flags_ = flags;
it->flags = flags;
*error = grib_iterator_init(it, h, args); GRIB_MUTEX_INIT_ONCE(&once, &init_mutex);
GRIB_MUTEX_LOCK(&mutex);
*error = it->init(h, args);
GRIB_MUTEX_UNLOCK(&mutex);
if (*error == GRIB_SUCCESS) if (*error == GRIB_SUCCESS)
return it; return it;
grib_context_log(h->context, GRIB_LOG_ERROR, "Geoiterator factory: Error instantiating iterator %s (%s)", grib_context_log(h->context, GRIB_LOG_ERROR, "Geoiterator factory: Error instantiating iterator %s (%s)",
table[i].type, grib_get_error_message(*error)); table[i].type, grib_get_error_message(*error));
grib_iterator_delete(it); gribIteratorDelete(it);
return NULL; return NULL;
} }
} }
@ -56,3 +91,119 @@ grib_iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsig
return NULL; return NULL;
} }
int grib_get_data(const grib_handle* h, double* lats, double* lons, double* values)
{
int err = 0;
eccodes::geo_iterator::Iterator* iter = NULL;
double *lat, *lon, *val;
iter = eccodes::geo_iterator::gribIteratorNew(h, 0, &err);
if (!iter || err != GRIB_SUCCESS)
return err;
lat = lats;
lon = lons;
val = values;
while (iter->next(lat++, lon++, val++)) {}
gribIteratorDelete(iter);
return err;
}
/*
* Return pointer to data at (i,j) (Fortran convention)
*/
static double* pointer_to_data(unsigned int i, unsigned int j,
long iScansNegatively, long jScansPositively,
long jPointsAreConsecutive, long alternativeRowScanning,
unsigned int nx, unsigned int ny, double* data)
{
/* Regular grid */
if (nx > 0 && ny > 0) {
if (i >= nx || j >= ny)
return NULL;
j = (jScansPositively) ? j : ny - 1 - j;
i = ((alternativeRowScanning) && (j % 2 == 1)) ? nx - 1 - i : i;
i = (iScansNegatively) ? nx - 1 - i : i;
return (jPointsAreConsecutive) ? data + j + i * ny : data + i + nx * j;
}
/* Reduced or other data not on a grid */
return NULL;
}
/* Apply the scanning mode flags which may require data array to be transformed
* to standard west-to-east (+i) south-to-north (+j) mode.
* The data array passed in should have 'numPoints' elements.
*/
int transform_iterator_data(grib_context* context, double* data,
long iScansNegatively, long jScansPositively,
long jPointsAreConsecutive, long alternativeRowScanning,
size_t numPoints, long nx, long ny)
{
double* data2;
double *pData0, *pData1, *pData2;
long ix, iy;
if (!iScansNegatively && jScansPositively && !jPointsAreConsecutive && !alternativeRowScanning) {
/* Already +i and +j. No need to change */
return GRIB_SUCCESS;
}
if (!data) return GRIB_SUCCESS;
if (!context) context = grib_context_get_default();
if (!iScansNegatively && !jScansPositively && !jPointsAreConsecutive && !alternativeRowScanning &&
nx > 0 && ny > 0) {
/* Regular grid +i -j: convert from we:ns to we:sn */
size_t row_size = ((size_t)nx) * sizeof(double);
data2 = (double*)grib_context_malloc(context, row_size);
if (!data2) {
grib_context_log(context, GRIB_LOG_ERROR, "Geoiterator data: Error allocating %ld bytes", row_size);
return GRIB_OUT_OF_MEMORY;
}
for (iy = 0; iy < ny / 2; iy++) {
memcpy(data2, data + ((size_t)iy) * nx, row_size);
memcpy(data + iy * nx, data + (ny - 1 - iy) * ((size_t)nx), row_size);
memcpy(data + (ny - 1 - iy) * ((size_t)nx), data2, row_size);
}
grib_context_free(context, data2);
return GRIB_SUCCESS;
}
if (nx < 1 || ny < 1) {
grib_context_log(context, GRIB_LOG_ERROR, "Geoiterator data: Invalid values for Nx and/or Ny");
return GRIB_GEOCALCULUS_PROBLEM;
}
data2 = (double*)grib_context_malloc(context, numPoints * sizeof(double));
if (!data2) {
grib_context_log(context, GRIB_LOG_ERROR, "Geoiterator data: Error allocating %ld bytes", numPoints * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
pData0 = data2;
for (iy = 0; iy < ny; iy++) {
long deltaX = 0;
pData1 = pointer_to_data(0, iy, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning, nx, ny, data);
if (!pData1) {
grib_context_free(context, data2);
return GRIB_GEOCALCULUS_PROBLEM;
}
pData2 = pointer_to_data(1, iy, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning, nx, ny, data);
if (!pData2) {
grib_context_free(context, data2);
return GRIB_GEOCALCULUS_PROBLEM;
}
deltaX = pData2 - pData1;
for (ix = 0; ix < nx; ix++) {
*pData0++ = *pData1;
pData1 += deltaX;
}
}
memcpy(data, data2, ((size_t)numPoints) * sizeof(double));
grib_context_free(context, data2);
return GRIB_SUCCESS;
}

View File

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

View File

@ -1,251 +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
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)
*/
static double* pointer_to_data(unsigned int i, unsigned int j,
long iScansNegatively, long jScansPositively,
long jPointsAreConsecutive, long alternativeRowScanning,
unsigned int nx, unsigned int ny, double* data)
{
/* Regular grid */
if (nx > 0 && ny > 0) {
if (i >= nx || j >= ny)
return NULL;
j = (jScansPositively) ? j : ny - 1 - j;
i = ((alternativeRowScanning) && (j % 2 == 1)) ? nx - 1 - i : i;
i = (iScansNegatively) ? nx - 1 - i : i;
return (jPointsAreConsecutive) ? data + j + i * ny : data + i + nx * j;
}
/* Reduced or other data not on a grid */
return NULL;
}
/* Apply the scanning mode flags which may require data array to be transformed
* to standard west-to-east (+i) south-to-north (+j) mode.
* The data array passed in should have 'numPoints' elements.
*/
int transform_iterator_data(grib_context* context, double* data,
long iScansNegatively, long jScansPositively,
long jPointsAreConsecutive, long alternativeRowScanning,
size_t numPoints, long nx, long ny)
{
double* data2;
double *pData0, *pData1, *pData2;
long ix, iy;
if (!iScansNegatively && jScansPositively && !jPointsAreConsecutive && !alternativeRowScanning) {
/* Already +i and +j. No need to change */
return GRIB_SUCCESS;
}
if (!data) return GRIB_SUCCESS;
if (!context) context = grib_context_get_default();
if (!iScansNegatively && !jScansPositively && !jPointsAreConsecutive && !alternativeRowScanning &&
nx > 0 && ny > 0) {
/* Regular grid +i -j: convert from we:ns to we:sn */
size_t row_size = ((size_t)nx) * sizeof(double);
data2 = (double*)grib_context_malloc(context, row_size);
if (!data2) {
grib_context_log(context, GRIB_LOG_ERROR, "Geoiterator data: Error allocating %ld bytes", row_size);
return GRIB_OUT_OF_MEMORY;
}
for (iy = 0; iy < ny / 2; iy++) {
memcpy(data2, data + ((size_t)iy) * nx, row_size);
memcpy(data + iy * nx, data + (ny - 1 - iy) * ((size_t)nx), row_size);
memcpy(data + (ny - 1 - iy) * ((size_t)nx), data2, row_size);
}
grib_context_free(context, data2);
return GRIB_SUCCESS;
}
if (nx < 1 || ny < 1) {
grib_context_log(context, GRIB_LOG_ERROR, "Geoiterator data: Invalid values for Nx and/or Ny");
return GRIB_GEOCALCULUS_PROBLEM;
}
data2 = (double*)grib_context_malloc(context, numPoints * sizeof(double));
if (!data2) {
grib_context_log(context, GRIB_LOG_ERROR, "Geoiterator data: Error allocating %ld bytes", numPoints * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
pData0 = data2;
for (iy = 0; iy < ny; iy++) {
long deltaX = 0;
pData1 = pointer_to_data(0, iy, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning, nx, ny, data);
if (!pData1) {
grib_context_free(context, data2);
return GRIB_GEOCALCULUS_PROBLEM;
}
pData2 = pointer_to_data(1, iy, iScansNegatively, jScansPositively, jPointsAreConsecutive, alternativeRowScanning, nx, ny, data);
if (!pData2) {
grib_context_free(context, data2);
return GRIB_GEOCALCULUS_PROBLEM;
}
deltaX = pData2 - pData1;
for (ix = 0; ix < nx; ix++) {
*pData0++ = *pData1;
pData1 += deltaX;
}
}
memcpy(data, data2, ((size_t)numPoints) * sizeof(double));
grib_context_free(context, data2);
return GRIB_SUCCESS;
}
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
{
grib_iterator_gen* self = (grib_iterator_gen*)iter;
size_t dli = 0;
int err = GRIB_SUCCESS;
const char* s_rawData = NULL;
const char* s_numPoints = NULL;
long numberOfPoints = 0;
self->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++);
iter->data = NULL;
iter->h = h; /* We may not need to keep them */
iter->args = args;
if ((err = grib_get_size(h, s_rawData, &dli)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(h, s_numPoints, &numberOfPoints)) != GRIB_SUCCESS)
return err;
// 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) {
// Iterator's number of values taken from the Grid Section
iter->nv = numberOfPoints;
} else {
// Check for consistency between the Grid and Data Sections
if (numberOfPoints != dli) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Geoiterator: %s != size(%s) (%ld!=%ld)",
s_numPoints, s_rawData, numberOfPoints, dli);
return GRIB_WRONG_GRID;
}
iter->nv = dli;
}
if (iter->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) {
// 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));
if ((err = grib_get_double_array_internal(h, s_rawData, iter->data, &(iter->nv)))) {
return err;
}
}
iter->e = -1;
return err;
}
static int reset(grib_iterator* iter)
{
iter->e = -1;
return 0;
}
static int destroy(grib_iterator* iter)
{
const grib_context* c = iter->h->context;
grib_context_free(c, iter->data);
return GRIB_SUCCESS;
}
static long has_next(grib_iterator* iter)
{
if (iter->flags == 0 && iter->data == NULL)
return 0;
if (iter->e >= (long)(iter->nv - 1))
return 0;
return 1;
}

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;
}

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;
}

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;
}

View File

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