diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index aa38b6ad3..aae376d94 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,6 +12,7 @@ include_directories( "${CMAKE_CURRENT_SOURCE_DIR}/accessor" "${CMAKE_CURRENT_SOURCE_DIR}/geo_iterator" + "${CMAKE_CURRENT_SOURCE_DIR}/geo_nearest" ) list( APPEND eccodes_src_files @@ -325,32 +326,31 @@ list( APPEND eccodes_src_files grib_expression_class_double.cc grib_expression_class_string.cc grib_expression_class_sub_string.cc - grib_nearest.cc - grib_nearest_class.cc - grib_nearest_class_gen.cc - grib_nearest_class_healpix.cc - grib_nearest_class_regular.cc - grib_nearest_class_reduced.cc - grib_nearest_class_latlon_reduced.cc - grib_nearest_class_lambert_conformal.cc - grib_nearest_class_lambert_azimuthal_equal_area.cc - grib_nearest_class_mercator.cc - grib_nearest_class_polar_stereographic.cc - grib_nearest_class_space_view.cc - geo_iterator/grib_iterator_class_polar_stereographic.cc - geo_iterator/grib_iterator_class_lambert_azimuthal_equal_area.cc - geo_iterator/grib_iterator_class_lambert_conformal.cc - geo_iterator/grib_iterator_class_mercator.cc + grib_nearest_factory.cc + geo_nearest/grib_nearest.cc + geo_nearest/grib_nearest_class_gen.cc + geo_nearest/grib_nearest_class_healpix.cc + geo_nearest/grib_nearest_class_lambert_azimuthal_equal_area.cc + geo_nearest/grib_nearest_class_lambert_conformal.cc + geo_nearest/grib_nearest_class_latlon_reduced.cc + geo_nearest/grib_nearest_class_mercator.cc + geo_nearest/grib_nearest_class_polar_stereographic.cc + geo_nearest/grib_nearest_class_reduced.cc + geo_nearest/grib_nearest_class_regular.cc + geo_nearest/grib_nearest_class_space_view.cc geo_iterator/grib_iterator.cc - grib_iterator_class.cc geo_iterator/grib_iterator_class_gaussian.cc geo_iterator/grib_iterator_class_gaussian_reduced.cc - geo_iterator/grib_iterator_class_latlon_reduced.cc geo_iterator/grib_iterator_class_gen.cc + geo_iterator/grib_iterator_class_healpix.cc + geo_iterator/grib_iterator_class_lambert_azimuthal_equal_area.cc + geo_iterator/grib_iterator_class_lambert_conformal.cc geo_iterator/grib_iterator_class_latlon.cc + geo_iterator/grib_iterator_class_latlon_reduced.cc + geo_iterator/grib_iterator_class_mercator.cc + geo_iterator/grib_iterator_class_polar_stereographic.cc geo_iterator/grib_iterator_class_regular.cc geo_iterator/grib_iterator_class_space_view.cc - geo_iterator/grib_iterator_class_healpix.cc grib_expression.cc codes_util.cc grib_util.cc @@ -363,10 +363,8 @@ list( APPEND eccodes_src_files eccodes_prototypes.h grib_dumper_class.h grib_dumper_factory.h - grib_iterator_class.h grib_iterator_factory.h - grib_nearest_class.h - grib_nearest_factory.h + grib_iterator_factory.cc grib_yacc.h md5.h md5.cc diff --git a/src/accessor/grib_accessor_class_nearest.cc b/src/accessor/grib_accessor_class_nearest.cc index f8d97d00a..31617c009 100644 --- a/src/accessor/grib_accessor_class_nearest.cc +++ b/src/accessor/grib_accessor_class_nearest.cc @@ -25,34 +25,3 @@ void grib_accessor_nearest_t::dump(grib_dumper* dumper) grib_dump_label(dumper, this, NULL); } -#if defined(HAVE_GEOGRAPHY) -grib_nearest* grib_nearest_new(const grib_handle* ch, int* error) -{ - grib_handle* h = (grib_handle*)ch; - grib_accessor* a = NULL; - grib_accessor_nearest_t* na = NULL; - grib_nearest* n = NULL; - *error = GRIB_NOT_IMPLEMENTED; - a = grib_find_accessor(h, "NEAREST"); - na = (grib_accessor_nearest_t*)a; - - if (!a) - return NULL; - - n = grib_nearest_factory(h, na->args_, error); - - if (n) - *error = GRIB_SUCCESS; - - return n; -} -#else -grib_nearest* grib_nearest_new(const grib_handle* ch, int* error) -{ - *error = GRIB_FUNCTIONALITY_NOT_ENABLED; - grib_context_log(ch->context, GRIB_LOG_ERROR, - "Nearest neighbour functionality not enabled. Please rebuild with -DENABLE_GEOGRAPHY=ON"); - - return NULL; -} -#endif diff --git a/src/accessor/grib_accessor_class_nearest.h b/src/accessor/grib_accessor_class_nearest.h index 9a012c349..15ed88e73 100644 --- a/src/accessor/grib_accessor_class_nearest.h +++ b/src/accessor/grib_accessor_class_nearest.h @@ -11,6 +11,7 @@ #pragma once #include "grib_accessor_class_gen.h" +#include "geo_nearest/grib_nearest.h" class grib_accessor_nearest_t : public grib_accessor_gen_t { @@ -23,8 +24,5 @@ public: private: grib_arguments* args_ = nullptr; - - friend grib_nearest* grib_nearest_new(const grib_handle* ch, int* error); + friend eccodes::geo_nearest::Nearest* eccodes::geo_nearest::gribNearestNew(const grib_handle* ch, int* error); }; - -// grib_nearest* grib_nearest_new(const grib_handle* ch, int* error); diff --git a/src/eccodes_prototypes.h b/src/eccodes_prototypes.h index 90592dd72..87aa35e69 100644 --- a/src/eccodes_prototypes.h +++ b/src/eccodes_prototypes.h @@ -792,23 +792,23 @@ grib_expression* new_string_expression(grib_context* c, const char* value); grib_expression* new_sub_string_expression(grib_context* c, const char* value, size_t start, size_t length); /* grib_nearest.cc */ -int grib_nearest_find(grib_nearest* nearest, const grib_handle* h, double inlat, double inlon, unsigned long flags, double* outlats, double* outlons, double* values, double* distances, int* indexes, size_t* len); -int grib_nearest_init(grib_nearest* i, grib_handle* h, grib_arguments* args); -int grib_nearest_delete(grib_nearest* i); -int grib_nearest_get_radius(grib_handle* h, double* radiusInKm); -void grib_binary_search(const double xx[], const size_t n, double x, size_t* ju, size_t* jl); -int grib_nearest_find_multiple(const grib_handle* h, int is_lsm, const double* inlats, const double* inlons, long npoints, double* outlats, double* outlons, double* values, double* distances, int* indexes); -int grib_nearest_find_generic(grib_nearest* nearest, grib_handle* h, double inlat, double inlon, unsigned long flags, - const char* values_keyname, - double** out_lats, - int* out_lats_count, - double** out_lons, - int* out_lons_count, - double** out_distances, - double* outlats, double* outlons, double* values, double* distances, int* indexes, size_t* len); +//int grib_nearest_find(grib_nearest* nearest, const grib_handle* h, double inlat, double inlon, unsigned long flags, double* outlats, double* outlons, double* values, double* distances, int* indexes, size_t* len); +//int grib_nearest_init(grib_nearest* i, grib_handle* h, grib_arguments* args); +//int grib_nearest_delete(grib_nearest* i); +//int grib_nearest_get_radius(grib_handle* h, double* radiusInKm); +//void grib_binary_search(const double xx[], const size_t n, double x, size_t* ju, size_t* jl); +//int grib_nearest_find_multiple(const grib_handle* h, int is_lsm, const double* inlats, const double* inlons, long npoints, double* outlats, double* outlons, double* values, double* distances, int* indexes); +//int grib_nearest_find_generic(grib_nearest* nearest, grib_handle* h, double inlat, double inlon, unsigned long flags, +// const char* values_keyname, +// double** out_lats, +// int* out_lats_count, +// double** out_lons, +// int* out_lons_count, +// double** out_distances, +// double* outlats, double* outlons, double* values, double* distances, int* indexes, size_t* len); /* grib_nearest_class.cc */ -grib_nearest* grib_nearest_factory(grib_handle* h, grib_arguments* args, int* error); +//eccodes::geo_nearest::Nearest* grib_nearest_factory(grib_handle* h, grib_arguments* args, int* error); /* grib_iterator.cc */ int grib_get_data(const grib_handle* h, double* lats, double* lons, double* values); diff --git a/src/grib_nearest.cc b/src/geo_nearest/grib_nearest.cc similarity index 81% rename from src/grib_nearest.cc rename to src/geo_nearest/grib_nearest.cc index d020fb269..81feae303 100644 --- a/src/grib_nearest.cc +++ b/src/geo_nearest/grib_nearest.cc @@ -8,95 +8,250 @@ * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. */ -/** -* Author: Enrico Fucile -* date: 31/10/2007 -* -*/ +#include "grib_nearest.h" +#include "grib_nearest_factory.h" +#include "accessor/grib_accessor_class_nearest.h" -#include "grib_api_internal.h" - -/* Note: The 'values' argument can be NULL in which case the data section will not be decoded - * See ECC-499 - */ -int grib_nearest_find( - grib_nearest* nearest, const grib_handle* ch, - double inlat, double inlon, - unsigned long flags, - double* outlats, double* outlons, - double* values, double* distances, int* indexes, size_t* len) +struct PointStore { - grib_handle* h = (grib_handle*)ch; - grib_nearest_class* c = NULL; - if (!nearest) - return GRIB_INVALID_ARGUMENT; - c = nearest->cclass; - Assert(flags <= (GRIB_NEAREST_SAME_GRID | GRIB_NEAREST_SAME_DATA | GRIB_NEAREST_SAME_POINT)); + double m_lat; + double m_lon; + double m_dist; + double m_value; + int m_index; +}; - while (c) { - grib_nearest_class* s = c->super ? *(c->super) : NULL; - if (c->find) { - int ret = c->find(nearest, h, inlat, inlon, flags, outlats, outlons, values, distances, indexes, len); - if (ret != GRIB_SUCCESS) { - if (inlon > 0) - inlon -= 360; - else - inlon += 360; - ret = c->find(nearest, h, inlat, inlon, flags, outlats, outlons, values, distances, indexes, len); - } - return ret; - } - c = s; +/* Generic implementation of nearest for Lambert, Polar stereo, Mercator etc */ +static int compare_doubles(const void* a, const void* b, int ascending) +{ + /* ascending is a boolean: 0 or 1 */ + const double* arg1 = (const double*)a; + const double* arg2 = (const double*)b; + if (ascending) { + if (*arg1 < *arg2) + return -1; /*Smaller values come before larger ones*/ } - Assert(0); + else { + if (*arg1 > *arg2) + return -1; /*Larger values come before smaller ones*/ + } + if (*arg1 == *arg2) + return 0; + else + return 1; +} + +static int compare_doubles_ascending(const void* a, const void* b) +{ + return compare_doubles(a, b, 1); +} + +/* Comparison function to sort points by distance */ +static int compare_points(const void* a, const void* b) +{ + const PointStore* pA = (const PointStore*)a; + const PointStore* pB = (const PointStore*)b; + + if (pA->m_dist < pB->m_dist) return -1; + if (pA->m_dist > pB->m_dist) return 1; return 0; } -/* For this one, ALL init are called */ -static int init_nearest(grib_nearest_class* c, grib_nearest* i, grib_handle* h, grib_arguments* args) +namespace eccodes::geo_nearest { + +int Nearest::init(grib_handle* h, grib_arguments* args) { - if (c) { - int ret = GRIB_SUCCESS; - grib_nearest_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_nearest(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_nearest_init(grib_nearest* i, grib_handle* h, grib_arguments* args) -{ - return init_nearest(i->cclass, i, h, args); + //h_ = h; + return GRIB_SUCCESS; } /* For this one, ALL destroy are called */ - -int grib_nearest_delete(grib_nearest* i) +int Nearest::destroy() { - grib_nearest_class* c = NULL; - if (!i) - return GRIB_INVALID_ARGUMENT; - c = i->cclass; - while (c) { - grib_nearest_class* s = c->super ? *(c->super) : NULL; - if (c->destroy) - c->destroy(i); - c = s; - } - return 0; + delete this; + return GRIB_SUCCESS; } +int Nearest::grib_nearest_find_generic( + grib_handle* h, + double inlat, double inlon, unsigned long flags, + + const char* values_keyname, + double** out_lats, + int* out_lats_count, + double** out_lons, + int* out_lons_count, + double** out_distances, + + double* outlats, double* outlons, + double* values, double* distances, int* indexes, size_t* len) +{ + int ret = 0; + size_t i = 0, nvalues = 0, nneighbours = 0; + double radiusInKm; + grib_iterator* iter = NULL; + double lat = 0, lon = 0; + + /* array of candidates for nearest neighbours */ + PointStore* neighbours = NULL; + + inlon = normalise_longitude_in_degrees(inlon); + + if ((ret = grib_get_size(h, values_keyname, &nvalues)) != GRIB_SUCCESS) + return ret; + values_count_ = nvalues; + + if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS) + return ret; + + neighbours = (PointStore*)grib_context_malloc(h->context, nvalues * sizeof(PointStore)); + for (i = 0; i < nvalues; ++i) { + neighbours[i].m_dist = 1e10; /* set all distances to large number to begin with */ + neighbours[i].m_lat = 0; + neighbours[i].m_lon = 0; + neighbours[i].m_value = 0; + neighbours[i].m_index = 0; + } + + /* GRIB_NEAREST_SAME_GRID not yet implemented */ + { + double the_value = 0; + double min_dist = 1e10; + size_t the_index = 0; + int ilat = 0, ilon = 0; + size_t idx_upper = 0, idx_lower = 0; + double lat1 = 0, lat2 = 0; /* inlat will be between these */ + const double LAT_DELTA = 10.0; /* in degrees */ + + *out_lons_count = (int)nvalues; /* Maybe overestimate but safe */ + *out_lats_count = (int)nvalues; + + if (*out_lats) + grib_context_free(h->context, *out_lats); + *out_lats = (double*)grib_context_malloc(h->context, nvalues * sizeof(double)); + if (!*out_lats) + return GRIB_OUT_OF_MEMORY; + + if (*out_lons) + grib_context_free(h->context, *out_lons); + *out_lons = (double*)grib_context_malloc(h->context, nvalues * sizeof(double)); + if (!*out_lons) + return GRIB_OUT_OF_MEMORY; + + iter = grib_iterator_new(h, 0, &ret); + if (ret) { + free(neighbours); + return ret; + } + /* First pass: collect all latitudes and longitudes */ + while (grib_iterator_next(iter, &lat, &lon, &the_value)) { + ++the_index; + Assert(ilat < *out_lats_count); + Assert(ilon < *out_lons_count); + (*out_lats)[ilat++] = lat; + (*out_lons)[ilon++] = lon; + } + + /* See between which 2 latitudes our point lies */ + qsort(*out_lats, nvalues, sizeof(double), &compare_doubles_ascending); + grib_binary_search(*out_lats, *out_lats_count - 1, inlat, &idx_upper, &idx_lower); + lat2 = (*out_lats)[idx_upper]; + lat1 = (*out_lats)[idx_lower]; + Assert(lat1 <= lat2); + + /* Second pass: Iterate again and collect candidate neighbours */ + grib_iterator_reset(iter); + the_index = 0; + i = 0; + while (grib_iterator_next(iter, &lat, &lon, &the_value)) { + if (lat > lat2 + LAT_DELTA || lat < lat1 - LAT_DELTA) { + /* Ignore latitudes too far from our point */ + } + else { + double dist = geographic_distance_spherical(radiusInKm, inlon, inlat, lon, lat); + if (dist < min_dist) + min_dist = dist; + /*printf("Candidate: lat=%.5f lon=%.5f dist=%f Idx=%ld Val=%f\n",lat,lon,dist,the_index,the_value);*/ + /* store this candidate point */ + neighbours[i].m_dist = dist; + neighbours[i].m_index = (int)the_index; + neighbours[i].m_lat = lat; + neighbours[i].m_lon = lon; + neighbours[i].m_value = the_value; + i++; + } + ++the_index; + } + nneighbours = i; + /* Sort the candidate neighbours in ascending order of distance */ + /* The first 4 entries will now be the closest 4 neighbours */ + qsort(neighbours, nneighbours, sizeof(PointStore), &compare_points); + + grib_iterator_delete(iter); + } + h_ = h; + + /* Sanity check for sorting */ +#ifdef DEBUG + for (i = 0; i < nneighbours - 1; ++i) { + Assert(neighbours[i].m_dist <= neighbours[i + 1].m_dist); + } +#endif + + /* GRIB_NEAREST_SAME_XXX not yet implemented */ + if (!*out_distances) { + *out_distances = (double*)grib_context_malloc(h->context, 4 * sizeof(double)); + } + (*out_distances)[0] = neighbours[0].m_dist; + (*out_distances)[1] = neighbours[1].m_dist; + (*out_distances)[2] = neighbours[2].m_dist; + (*out_distances)[3] = neighbours[3].m_dist; + + for (i = 0; i < 4; ++i) { + distances[i] = neighbours[i].m_dist; + outlats[i] = neighbours[i].m_lat; + outlons[i] = neighbours[i].m_lon; + indexes[i] = neighbours[i].m_index; + if (values) { + values[i] = neighbours[i].m_value; + } + /*printf("(%f,%f) i=%d d=%f v=%f\n",outlats[i],outlons[i],indexes[i],distances[i],values[i]);*/ + } + + free(neighbours); + return GRIB_SUCCESS; +} + +eccodes::geo_nearest::Nearest* gribNearestNew(const grib_handle* ch, int* error) +{ + *error = GRIB_NOT_IMPLEMENTED; + + grib_handle* h = (grib_handle*)ch; + grib_accessor* a = grib_find_accessor(h, "NEAREST"); + grib_accessor_nearest_t* n = (grib_accessor_nearest_t*)a; + + if (!a) + return NULL; + + eccodes::geo_nearest::Nearest* nearest = grib_nearest_factory(h, n->args_, error); + + if (nearest) + *error = GRIB_SUCCESS; + + return nearest; +} + +int gribNearestDelete(eccodes::geo_nearest::Nearest* i) +{ + if (i) + i->destroy(); + return GRIB_SUCCESS; +} + + +} // namespace eccodes::geo_nearest + + /* Get the radius in kilometres for nearest neighbour distance calculations */ /* For an ellipsoid, approximate the radius using the average of the semimajor and semiminor axes */ int grib_nearest_get_radius(grib_handle* h, double* radiusInKm) @@ -105,8 +260,6 @@ int grib_nearest_get_radius(grib_handle* h, double* radiusInKm) long lRadiusInMetres; double result = 0; const char* s_radius = "radius"; - const char* s_minor = "earthMinorAxisInMetres"; - const char* s_major = "earthMajorAxisInMetres"; if ((err = grib_get_long(h, s_radius, &lRadiusInMetres)) == GRIB_SUCCESS) { if (grib_is_missing(h, s_radius, &err) || lRadiusInMetres == GRIB_MISSING_LONG) { @@ -117,6 +270,8 @@ int grib_nearest_get_radius(grib_handle* h, double* radiusInKm) } else { double minor = 0, major = 0; + const char* s_minor = "earthMinorAxisInMetres"; + const char* s_major = "earthMajorAxisInMetres"; if ((err = grib_get_double_internal(h, s_minor, &minor)) != GRIB_SUCCESS) return err; if ((err = grib_get_double_internal(h, s_major, &major)) != GRIB_SUCCESS) return err; if (grib_is_missing(h, s_minor, &err)) return GRIB_GEOCALCULUS_PROBLEM; @@ -145,6 +300,13 @@ void grib_binary_search(const double xx[], const size_t n, double x, size_t* ju, } } + +/* + * C API Implementation + * codes_iterator_* wrappers are in eccodes.h and eccodes.cc + * grib_iterator_* declarations are in grib_api.h + */ + int grib_nearest_find_multiple( const grib_handle* h, int is_lsm, const double* inlats, const double* inlons, long npoints, @@ -244,198 +406,65 @@ int grib_nearest_find_multiple( return ret; } - -/* Generic implementation of nearest for Lambert, Polar stereo, Mercator etc */ -static int compare_doubles(const void* a, const void* b, int ascending) -{ - /* ascending is a boolean: 0 or 1 */ - double* arg1 = (double*)a; - double* arg2 = (double*)b; - if (ascending) { - if (*arg1 < *arg2) - return -1; /*Smaller values come before larger ones*/ - } - else { - if (*arg1 > *arg2) - return -1; /*Larger values come before smaller ones*/ - } - if (*arg1 == *arg2) - return 0; - else - return 1; -} - -static int compare_doubles_ascending(const void* a, const void* b) -{ - return compare_doubles(a, b, 1); -} - -typedef struct PointStore -{ - double m_lat; - double m_lon; - double m_dist; - double m_value; - int m_index; -} PointStore; - -/* Comparison function to sort points by distance */ -static int compare_points(const void* a, const void* b) -{ - PointStore* pA = (PointStore*)a; - PointStore* pB = (PointStore*)b; - - if (pA->m_dist < pB->m_dist) return -1; - if (pA->m_dist > pB->m_dist) return 1; - return 0; -} - -int grib_nearest_find_generic( - grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - - const char* values_keyname, - double** out_lats, - int* out_lats_count, - double** out_lons, - int* out_lons_count, - double** out_distances, - +/* Note: The 'values' argument can be NULL in which case the data section will not be decoded + * See ECC-499 + */ +int grib_nearest_find( + grib_nearest* nearest, const grib_handle* ch, + double inlat, double inlon, + unsigned long flags, double* outlats, double* outlons, double* values, double* distances, int* indexes, size_t* len) { - int ret = 0; - size_t i = 0, nvalues = 0, nneighbours = 0; - double radiusInKm; - grib_iterator* iter = NULL; - double lat = 0, lon = 0; + grib_handle* h = (grib_handle*)ch; + if (!nearest) + return GRIB_INVALID_ARGUMENT; + Assert(flags <= (GRIB_NEAREST_SAME_GRID | GRIB_NEAREST_SAME_DATA | GRIB_NEAREST_SAME_POINT)); - /* array of candidates for nearest neighbours */ - PointStore* neighbours = NULL; - - inlon = normalise_longitude_in_degrees(inlon); - - if ((ret = grib_get_size(h, values_keyname, &nvalues)) != GRIB_SUCCESS) - return ret; - nearest->values_count = nvalues; - - if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS) - return ret; - - neighbours = (PointStore*)grib_context_malloc(h->context, nvalues * sizeof(PointStore)); - for (i = 0; i < nvalues; ++i) { - neighbours[i].m_dist = 1e10; /* set all distances to large number to begin with */ - neighbours[i].m_lat = 0; - neighbours[i].m_lon = 0; - neighbours[i].m_value = 0; - neighbours[i].m_index = 0; + int ret = nearest->nearest->find(h, inlat, inlon, flags, outlats, outlons, values, distances, indexes, len); + if (ret != GRIB_SUCCESS) { + if (inlon > 0) + inlon -= 360; + else + inlon += 360; + ret = nearest->nearest->find(h, inlat, inlon, flags, outlats, outlons, values, distances, indexes, len); } + return ret; +} - /* GRIB_NEAREST_SAME_GRID not yet implemented */ - { - double the_value = 0; - double min_dist = 1e10; - size_t the_index = 0; - int ilat = 0, ilon = 0; - size_t idx_upper = 0, idx_lower = 0; - double lat1 = 0, lat2 = 0; /* inlat will be between these */ - const double LAT_DELTA = 10.0; /* in degrees */ +int grib_nearest_init(grib_nearest* i, grib_handle* h, grib_arguments* args) +{ + return i->nearest->init(h, args); +} - *out_lons_count = (int)nvalues; /* Maybe overestimate but safe */ - *out_lats_count = (int)nvalues; - - if (*out_lats) - grib_context_free(h->context, *out_lats); - *out_lats = (double*)grib_context_malloc(h->context, nvalues * sizeof(double)); - if (!*out_lats) - return GRIB_OUT_OF_MEMORY; - - if (*out_lons) - grib_context_free(h->context, *out_lons); - *out_lons = (double*)grib_context_malloc(h->context, nvalues * sizeof(double)); - if (!*out_lons) - return GRIB_OUT_OF_MEMORY; - - iter = grib_iterator_new(h, 0, &ret); - if (ret) { - free(neighbours); - return ret; - } - /* First pass: collect all latitudes and longitudes */ - while (grib_iterator_next(iter, &lat, &lon, &the_value)) { - ++the_index; - Assert(ilat < *out_lats_count); - Assert(ilon < *out_lons_count); - (*out_lats)[ilat++] = lat; - (*out_lons)[ilon++] = lon; - } - - /* See between which 2 latitudes our point lies */ - qsort(*out_lats, nvalues, sizeof(double), &compare_doubles_ascending); - grib_binary_search(*out_lats, *out_lats_count - 1, inlat, &idx_upper, &idx_lower); - lat2 = (*out_lats)[idx_upper]; - lat1 = (*out_lats)[idx_lower]; - Assert(lat1 <= lat2); - - /* Second pass: Iterate again and collect candidate neighbours */ - grib_iterator_reset(iter); - the_index = 0; - i = 0; - while (grib_iterator_next(iter, &lat, &lon, &the_value)) { - if (lat > lat2 + LAT_DELTA || lat < lat1 - LAT_DELTA) { - /* Ignore latitudes too far from our point */ - } - else { - double dist = geographic_distance_spherical(radiusInKm, inlon, inlat, lon, lat); - if (dist < min_dist) - min_dist = dist; - /*printf("Candidate: lat=%.5f lon=%.5f dist=%f Idx=%ld Val=%f\n",lat,lon,dist,the_index,the_value);*/ - /* store this candidate point */ - neighbours[i].m_dist = dist; - neighbours[i].m_index = (int)the_index; - neighbours[i].m_lat = lat; - neighbours[i].m_lon = lon; - neighbours[i].m_value = the_value; - i++; - } - ++the_index; - } - nneighbours = i; - /* Sort the candidate neighbours in ascending order of distance */ - /* The first 4 entries will now be the closest 4 neighbours */ - qsort(neighbours, nneighbours, sizeof(PointStore), &compare_points); - - grib_iterator_delete(iter); +int grib_nearest_delete(grib_nearest* i) +{ + if (i) { + grib_context* c = grib_context_get_default(); + gribNearestDelete(i->nearest); + grib_context_free(c, i); } - nearest->h = h; - - /* Sanity check for sorting */ -#ifdef DEBUG - for (i = 0; i < nneighbours - 1; ++i) { - Assert(neighbours[i].m_dist <= neighbours[i + 1].m_dist); - } -#endif - - /* GRIB_NEAREST_SAME_XXX not yet implemented */ - if (!*out_distances) { - *out_distances = (double*)grib_context_malloc(h->context, 4 * sizeof(double)); - } - (*out_distances)[0] = neighbours[0].m_dist; - (*out_distances)[1] = neighbours[1].m_dist; - (*out_distances)[2] = neighbours[2].m_dist; - (*out_distances)[3] = neighbours[3].m_dist; - - for (i = 0; i < 4; ++i) { - distances[i] = neighbours[i].m_dist; - outlats[i] = neighbours[i].m_lat; - outlons[i] = neighbours[i].m_lon; - indexes[i] = neighbours[i].m_index; - if (values) { - values[i] = neighbours[i].m_value; - } - /*printf("(%f,%f) i=%d d=%f v=%f\n",outlats[i],outlons[i],indexes[i],distances[i],values[i]);*/ - } - - free(neighbours); return GRIB_SUCCESS; } + +#if defined(HAVE_GEOGRAPHY) +grib_nearest* grib_nearest_new(const grib_handle* ch, int* error) +{ + grib_nearest* i = (grib_nearest*)grib_context_malloc_clear(ch->context, sizeof(grib_nearest)); + i->nearest = eccodes::geo_nearest::gribNearestNew(ch, error); + if (!i->nearest) { + grib_context_free(ch->context, i); + return NULL; + } + return i; +} +#else +grib_nearest* grib_nearest_new(const grib_handle* ch, int* error) +{ + *error = GRIB_FUNCTIONALITY_NOT_ENABLED; + grib_context_log(ch->context, GRIB_LOG_ERROR, + "Nearest neighbour functionality not enabled. Please rebuild with -DENABLE_GEOGRAPHY=ON"); + + return NULL; +} +#endif diff --git a/src/geo_nearest/grib_nearest.h b/src/geo_nearest/grib_nearest.h new file mode 100644 index 000000000..74cf9a0ff --- /dev/null +++ b/src/geo_nearest/grib_nearest.h @@ -0,0 +1,48 @@ +/* + * (C) Copyright 2005- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * + * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by + * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. + */ + +#pragma once + +#include "grib_api_internal.h" + +namespace eccodes::geo_nearest { + +class Nearest { +public: + virtual ~Nearest() {} + virtual int init(grib_handle*, grib_arguments*) = 0; + virtual int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) = 0; + virtual int destroy() = 0; + virtual Nearest* create() = 0; + +protected: + int grib_nearest_find_generic(grib_handle*, double, double, unsigned long, + const char*, + double**, + int*, + double**, + int*, + double**, + double*, double*, double*, double*, int*, size_t*); + + grib_handle* h_ = nullptr; + double* values_ = nullptr; + size_t values_count_ = 0; + unsigned long flags_ = 0; + const char* class_name_ = nullptr; +}; + +Nearest* gribNearestNew(const grib_handle*, int*); +int gribNearestDelete(Nearest*); + +} // namespace eccodes::geo_nearest + +int grib_nearest_get_radius(grib_handle* h, double* radiusInKm); +void grib_binary_search(const double xx[], const size_t n, double x, size_t* ju, size_t* jl); diff --git a/src/geo_nearest/grib_nearest_class_gen.cc b/src/geo_nearest/grib_nearest_class_gen.cc new file mode 100644 index 000000000..16e31508e --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_gen.cc @@ -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. + */ + +#include "grib_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +int Gen::init(grib_handle* h, grib_arguments* args) +{ + int ret = GRIB_SUCCESS; + if ((ret = Nearest::init(h, args) != GRIB_SUCCESS)) + return ret; + + cargs_ = 1; + + values_key_ = grib_arguments_get_name(h, args, cargs_++); + radius_ = grib_arguments_get_name(h, args, cargs_++); + values_ = NULL; + + return ret; +} + +int Gen::destroy() +{ + grib_context* c = grib_context_get_default(); + if (values_) + grib_context_free(c, values_); + + return Nearest::destroy(); +} + +int Gen::find(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, double* values, + double* distances, int* indexes, size_t* len) +{ + return GRIB_NOT_IMPLEMENTED; +} + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_gen.h b/src/geo_nearest/grib_nearest_class_gen.h new file mode 100644 index 000000000..8ef4fa05c --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_gen.h @@ -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_nearest.h" + +namespace eccodes::geo_nearest { + +class Gen : public Nearest { +public: + Gen() { class_name_ = "gen"; } + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +protected: + int cargs_ = 0; + const char* values_key_ = nullptr; + +private: + const char* radius_ = nullptr; +}; + +int grib_nearest_find_generic( + Nearest* nearest, grib_handle* h, + double inlat, double inlon, unsigned long flags, + + const char* values_keyname, + double** out_lats, + int* out_lats_count, + double** out_lons, + int* out_lons_count, + double** out_distances, + + double* outlats, double* outlons, + double* values, double* distances, int* indexes, size_t* len); + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_healpix.cc b/src/geo_nearest/grib_nearest_class_healpix.cc new file mode 100644 index 000000000..d55f220d4 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_healpix.cc @@ -0,0 +1,65 @@ +/* + * (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_nearest_class_healpix.h" + +eccodes::geo_nearest::Healpix _grib_nearest_healpix{}; +eccodes::geo_nearest::Healpix* grib_nearest_healpix = &_grib_nearest_healpix; + +namespace eccodes::geo_nearest { + +int Healpix::init(grib_handle* h, grib_arguments* args) +{ + int ret = GRIB_SUCCESS; + if ((ret = Gen::init(h, args) != GRIB_SUCCESS)) + return ret; + + Ni_ = grib_arguments_get_name(h, args, cargs_++); + Nj_ = grib_arguments_get_name(h, args, cargs_++); + i_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + j_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + + return ret; +} + +int Healpix::find(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, + double* values, double* distances, int* indexes, size_t* len) +{ + return grib_nearest_find_generic( + h, inlat, inlon, flags, /* inputs */ + + values_key_, /* outputs to set the 'self' object */ + &(lats_), + &(lats_count_), + &(lons_), + &(lons_count_), + &(distances_), + + outlats, outlons, /* outputs of the find function */ + values, distances, indexes, len); +} + +int Healpix::destroy() +{ + grib_context* c = grib_context_get_default(); + + if (lats_) grib_context_free(c, lats_); + if (lons_) grib_context_free(c, lons_); + if (i_) grib_context_free(c, i_); + if (j_) grib_context_free(c, j_); + if (k_) grib_context_free(c, k_); + if (distances_) grib_context_free(c, distances_); + + return GRIB_SUCCESS; +} + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_healpix.h b/src/geo_nearest/grib_nearest_class_healpix.h new file mode 100644 index 000000000..709f8cb56 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_healpix.h @@ -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_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +class Healpix : public Gen { +public: + Healpix() { + class_name_ = "healpix"; + } + Nearest* create() override { return new Healpix(); }; + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +private: + double* lats_ = nullptr; + int lats_count_ = 0; + double* lons_ = nullptr; + int lons_count_ = 0; + double* distances_ = nullptr; + int* k_ = nullptr; + int* i_ = nullptr; + int* j_ = nullptr; + const char* Ni_ = nullptr; + const char* Nj_ = nullptr; +}; + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_lambert_azimuthal_equal_area.cc b/src/geo_nearest/grib_nearest_class_lambert_azimuthal_equal_area.cc new file mode 100644 index 000000000..d2d4d4523 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_lambert_azimuthal_equal_area.cc @@ -0,0 +1,64 @@ +/* + * (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_nearest_class_lambert_azimuthal_equal_area.h" + +eccodes::geo_nearest::LambertAzimuthalEqualArea _grib_nearest_lambert_azimuthal_equal_area{}; +eccodes::geo_nearest::LambertAzimuthalEqualArea* grib_nearest_lambert_azimuthal_equal_area = &_grib_nearest_lambert_azimuthal_equal_area; + +namespace eccodes::geo_nearest { + +int LambertAzimuthalEqualArea::init(grib_handle* h, grib_arguments* args) +{ + int ret = GRIB_SUCCESS; + if ((ret = Gen::init(h, args) != GRIB_SUCCESS)) + return ret; + + Ni_ = grib_arguments_get_name(h, args, cargs_++); + Nj_ = grib_arguments_get_name(h, args, cargs_++); + i_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + j_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + return ret; +} + +int LambertAzimuthalEqualArea::find(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, + double* values, double* distances, int* indexes, size_t* len) +{ + return grib_nearest_find_generic( + h, inlat, inlon, flags, /* inputs */ + + values_key_, /* outputs to set the 'self' object */ + &(lats_), + &(lats_count_), + &(lons_), + &(lons_count_), + &(distances_), + + outlats, outlons, /* outputs of the find function */ + values, distances, indexes, len); +} + +int LambertAzimuthalEqualArea::destroy() +{ + grib_context* c = grib_context_get_default(); + + if (lats_) grib_context_free(c, lats_); + if (lons_) grib_context_free(c, lons_); + if (i_) grib_context_free(c, i_); + if (j_) grib_context_free(c, j_); + if (k_) grib_context_free(c, k_); + if (distances_) grib_context_free(c, distances_); + + return Gen::destroy(); +} + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_lambert_azimuthal_equal_area.h b/src/geo_nearest/grib_nearest_class_lambert_azimuthal_equal_area.h new file mode 100644 index 000000000..ae3a2403a --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_lambert_azimuthal_equal_area.h @@ -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_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +class LambertAzimuthalEqualArea : public Gen { +public: + LambertAzimuthalEqualArea() { + class_name_ = "lambert_azimuthal_equal_area"; + } + Nearest* create() override { return new LambertAzimuthalEqualArea(); } + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +private: + double* lats_ = nullptr; + int lats_count_ = 0; + double* lons_ = nullptr; + int lons_count_ = 0; + double* distances_ = nullptr; + int* k_ = nullptr; + int* i_ = nullptr; + int* j_ = nullptr; + const char* Ni_ = nullptr; + const char* Nj_ = nullptr; +}; + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_lambert_conformal.cc b/src/geo_nearest/grib_nearest_class_lambert_conformal.cc new file mode 100644 index 000000000..eb0241020 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_lambert_conformal.cc @@ -0,0 +1,65 @@ +/* + * (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_nearest_class_lambert_conformal.h" + +eccodes::geo_nearest::LambertConformal _grib_nearest_lambert_conformal{}; +eccodes::geo_nearest::LambertConformal* grib_nearest_lambert_conformal = &_grib_nearest_lambert_conformal; + +namespace eccodes::geo_nearest { + +int LambertConformal::init(grib_handle* h, grib_arguments* args) +{ + int ret = GRIB_SUCCESS; + if ((ret = Gen::init(h, args) != GRIB_SUCCESS)) + return ret; + + Ni_ = grib_arguments_get_name(h, args, cargs_++); + Nj_ = grib_arguments_get_name(h, args, cargs_++); + i_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + j_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + + return ret; +} + +int LambertConformal::find(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, + double* values, double* distances, int* indexes, size_t* len) +{ + return grib_nearest_find_generic( + h, inlat, inlon, flags, /* inputs */ + + values_key_, /* outputs to set the 'self' object */ + &(lats_), + &(lats_count_), + &(lons_), + &(lons_count_), + &(distances_), + + outlats, outlons, /* outputs of the find function */ + values, distances, indexes, len); +} + +int LambertConformal::destroy() +{ + grib_context* c = grib_context_get_default(); + + if (lats_) grib_context_free(c, lats_); + if (lons_) grib_context_free(c, lons_); + if (i_) grib_context_free(c, i_); + if (j_) grib_context_free(c, j_); + if (k_) grib_context_free(c, k_); + if (distances_) grib_context_free(c, distances_); + + return Gen::destroy(); +} + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_lambert_conformal.h b/src/geo_nearest/grib_nearest_class_lambert_conformal.h new file mode 100644 index 000000000..fc1970399 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_lambert_conformal.h @@ -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_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +class LambertConformal : public Gen { +public: + LambertConformal() { + class_name_ = "lambert_conformal"; + } + Nearest* create() override { return new LambertConformal(); }; + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +private: + double* lats_ = nullptr; + int lats_count_ = 0; + double* lons_ = nullptr; + int lons_count_ = 0; + double* distances_ = nullptr; + int* k_ = nullptr; + int* i_ = nullptr; + int* j_ = nullptr; + const char* Ni_ = nullptr; + const char* Nj_ = nullptr; +}; + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_latlon_reduced.cc b/src/geo_nearest/grib_nearest_class_latlon_reduced.cc new file mode 100644 index 000000000..024e5d3a9 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_latlon_reduced.cc @@ -0,0 +1,360 @@ +/* + * (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_nearest_class_latlon_reduced.h" + +eccodes::geo_nearest::LatlonReduced _grib_nearest_latlon_reduced{}; +eccodes::geo_nearest::LatlonReduced* grib_nearest_latlon_reduced = &_grib_nearest_latlon_reduced; + +namespace eccodes::geo_nearest { + +int LatlonReduced::init(grib_handle* h, grib_arguments* args) +{ + int ret = GRIB_SUCCESS; + if ((ret = Gen::init(h, args) != GRIB_SUCCESS)) + return ret; + + Nj_ = grib_arguments_get_name(h, args, cargs_++); + pl_ = grib_arguments_get_name(h, args, cargs_++); + lonFirst_ = grib_arguments_get_name(h, args, cargs_++); + lonLast_ = grib_arguments_get_name(h, args, cargs_++); + j_ = (size_t*)grib_context_malloc(h->context, 2 * sizeof(size_t)); + if (!j_) + return GRIB_OUT_OF_MEMORY; + k_ = (size_t*)grib_context_malloc(h->context, 4 * sizeof(size_t)); + if (!k_) + return GRIB_OUT_OF_MEMORY; + + return ret; +} + +int LatlonReduced::find(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, double* values, + double* distances, int* indexes, size_t* len) +{ + int err = 0; + double lat1, lat2, lon1, lon2; + int is_global = 1; + + if (grib_get_double(h, "longitudeFirstInDegrees", &lon1) == GRIB_SUCCESS && + grib_get_double(h, "longitudeLastInDegrees", &lon2) == GRIB_SUCCESS && + grib_get_double(h, "latitudeFirstInDegrees", &lat1) == GRIB_SUCCESS && + grib_get_double(h, "latitudeLastInDegrees", &lat2) == GRIB_SUCCESS) + { + const double difflat = fabs(lat1-lat2); + if (difflat < 180 || lon1 != 0 || lon2 < 359) { + is_global = 0; /* subarea */ + } + } + + if (is_global) { + err = find_global(h, inlat, inlon, flags, + outlats, outlons, values, + distances, indexes, len); + } + else + { + int lons_count = 0; /*dummy*/ + err = grib_nearest_find_generic( + h, inlat, inlon, flags, + values_key_, + &(lats_), + &(lats_count_), + &(lons_), + &(lons_count), + &(distances_), + outlats, outlons, + values, distances, indexes, len); + } + return err; +} + +int LatlonReduced::find_global(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, double* values, + double* distances, int* indexes, size_t* len) +{ + int ret = 0, kk = 0, ii = 0, jj = 0; + int j = 0; + long* pla = NULL; + long* pl = NULL; + size_t nvalues = 0; + grib_iterator* iter = NULL; + double lat = 0, lon = 0; + double radiusInKm; + int ilat = 0, ilon = 0; + + if ((ret = grib_get_size(h, values_key_, &nvalues)) != GRIB_SUCCESS) + return ret; + values_count_ = nvalues; + + if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS) + return ret; + + /* Compute lat/lon info, create iterator etc if it's the 1st time or different grid. + * This is for performance: if the grid has not changed, we only do this once + * and reuse for other messages */ + if (!h_ || (flags & GRIB_NEAREST_SAME_GRID) == 0) { + double olat = 1.e10; + long n = 0; + + ilat = 0; + ilon = 0; + if (grib_is_missing(h, Nj_, &ret)) { + grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Nj_); + return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; + } + + if ((ret = grib_get_long(h, Nj_, &n)) != GRIB_SUCCESS) + return ret; + lats_count_ = n; + + if (lats_) + grib_context_free(h->context, lats_); + lats_ = (double*)grib_context_malloc(h->context, + lats_count_ * sizeof(double)); + if (!lats_) + return GRIB_OUT_OF_MEMORY; + + if (lons_) + grib_context_free(h->context, lons_); + lons_ = (double*)grib_context_malloc(h->context, + values_count_ * sizeof(double)); + if (!lons_) + return GRIB_OUT_OF_MEMORY; + + iter = grib_iterator_new(h, GRIB_GEOITERATOR_NO_VALUES, &ret); + if (ret) { + grib_context_log(h->context, GRIB_LOG_ERROR, "Unable to create iterator"); + return ret; + } + while (grib_iterator_next(iter, &lat, &lon, NULL)) { + if (ilat < lats_count_ && olat != lat) { + lats_[ilat++] = lat; + olat = lat; + } + lons_[ilon++] = lon; + } + lats_count_ = ilat; + grib_iterator_delete(iter); + } + h_ = h; + + /* Compute distances if it's the 1st time or different point or different grid. + * This is for performance: if the grid and the input point have not changed + * we only do this once and reuse for other messages */ + if (!distances_ || (flags & GRIB_NEAREST_SAME_POINT) == 0 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { + double* lons = NULL; + int nlon = 0; + size_t plsize = 0; + long nplm1 = 0; + int nearest_lons_found = 0; + double lon_first, lon_last; + int islocal = 0; + long plmax; + double dimin; + + if ((ret = grib_get_double(h, lonFirst_, &lon_first)) != GRIB_SUCCESS) { + grib_context_log(h->context, GRIB_LOG_ERROR, + "grib_nearest_latlon_reduced.find(): unable to get %s %s\n", lonFirst_, + grib_get_error_message(ret)); + return ret; + } + if ((ret = grib_get_double(h, lonLast_, &lon_last)) != GRIB_SUCCESS) { + grib_context_log(h->context, GRIB_LOG_ERROR, + "grib_nearest_latlon_reduced.find(): unable to get %s %s\n", lonLast_, + grib_get_error_message(ret)); + return ret; + } + + plsize = lats_count_; + if ((ret = grib_get_size(h, pl_, &plsize)) != GRIB_SUCCESS) + return ret; + pla = (long*)grib_context_malloc(h->context, plsize * sizeof(long)); + if (!pla) + return GRIB_OUT_OF_MEMORY; + if ((ret = grib_get_long_array(h, pl_, pla, &plsize)) != GRIB_SUCCESS) + return ret; + + pl = pla; + while ((*pl) == 0) { + pl++; + } + + plmax = pla[0]; + for (j = 0; j < plsize; j++) + if (plmax < pla[j]) + plmax = pla[j]; + dimin = 360.0 / plmax; + + if (360 - fabs(lon_last - lon_first) < 2 * dimin) { + islocal = 0; + } + else { + islocal = 1; + } + + if (islocal) + for (j = 0; j < plsize; j++) + pla[j]--; + + /* printf("XXXX islocal=%d\n",islocal); */ + while (inlon < 0) + inlon += 360; + while (inlon > 360) + inlon -= 360; + + ilat = lats_count_; + if (lats_[ilat - 1] > lats_[0]) { + if (inlat < lats_[0] || inlat > lats_[ilat - 1]) + return GRIB_OUT_OF_AREA; + } + else { + if (inlat > lats_[0] || inlat < lats_[ilat - 1]) + return GRIB_OUT_OF_AREA; + } + + if (!distances_) + distances_ = (double*)grib_context_malloc(h->context, 4 * sizeof(double)); + if (!distances_) + return GRIB_OUT_OF_MEMORY; + + grib_binary_search(lats_, ilat - 1, inlat, + &(j_[0]), &(j_[1])); + + nlon = 0; + for (jj = 0; jj < j_[0]; jj++) + nlon += pl[jj]; + nplm1 = pl[j_[0]] - 1; + + lons = lons_ + nlon; + + nearest_lons_found = 0; + if (lons[nplm1] > lons[0]) { + if (inlon < lons[0] || inlon > lons[nplm1]) { + if (lons[nplm1] - lons[0] - 360 <= + lons[nplm1] - lons[nplm1 - 1]) { + k_[0] = 0; + k_[1] = nplm1; + nearest_lons_found = 1; + } + else + return GRIB_OUT_OF_AREA; + } + } + else { + if (inlon > lons[0] || inlon < lons[nplm1]) { + if (lons[0] - lons[nplm1] - 360 <= + lons[0] - lons[1]) { + k_[0] = 0; + k_[1] = nplm1; + nearest_lons_found = 1; + } + else + return GRIB_OUT_OF_AREA; + } + } + + if (!nearest_lons_found) { + grib_binary_search(lons, pl[j_[0]] - 1, inlon, + &(k_[0]), &(k_[1])); + } + k_[0] += nlon; + k_[1] += nlon; + + nlon = 0; + for (jj = 0; jj < j_[1]; jj++) + nlon += pl[jj]; + nplm1 = pl[j_[1]] - 1; + + lons = lons_ + nlon; + + nearest_lons_found = 0; + if (lons[nplm1] > lons[0]) { + if (inlon < lons[0] || inlon > lons[nplm1]) { + if (lons[nplm1] - lons[0] - 360 <= + lons[nplm1] - lons[nplm1 - 1]) { + k_[2] = 0; + k_[3] = nplm1; + nearest_lons_found = 1; + } + else + return GRIB_OUT_OF_AREA; + } + } + else { + if (inlon > lons[0] || inlon < lons[nplm1]) { + if (lons[0] - lons[nplm1] - 360 <= + lons[0] - lons[1]) { + k_[2] = 0; + k_[3] = nplm1; + nearest_lons_found = 1; + } + else + return GRIB_OUT_OF_AREA; + } + } + + if (!nearest_lons_found) { + grib_binary_search(lons, pl[j_[1]] - 1, inlon, + &(k_[2]), &(k_[3])); + } + + k_[2] += nlon; + k_[3] += nlon; + + kk = 0; + for (jj = 0; jj < 2; jj++) { + for (ii = 0; ii < 2; ii++) { + distances_[kk] = geographic_distance_spherical(radiusInKm, inlon, inlat, + lons_[k_[kk]], lats_[j_[jj]]); + kk++; + } + } + + grib_context_free(h->context, pla); + } + + kk = 0; + for (jj = 0; jj < 2; jj++) { + for (ii = 0; ii < 2; ii++) { + distances[kk] = distances_[kk]; + outlats[kk] = lats_[j_[jj]]; + outlons[kk] = lons_[k_[kk]]; + if (values) { /* ECC-499 */ + grib_get_double_element_internal(h, values_key_, k_[kk], &(values[kk])); + } + indexes[kk] = (int)k_[kk]; + kk++; + } + } + + return GRIB_SUCCESS; +} + +int LatlonReduced::destroy() +{ + grib_context* c = grib_context_get_default(); + if (lats_) + grib_context_free(c, lats_); + if (lons_) + grib_context_free(c, lons_); + if (j_) + grib_context_free(c, j_); + if (k_) + grib_context_free(c, k_); + if (distances_) + grib_context_free(c, distances_); + + return Gen::destroy(); +} + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_latlon_reduced.h b/src/geo_nearest/grib_nearest_class_latlon_reduced.h new file mode 100644 index 000000000..743f7d4a3 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_latlon_reduced.h @@ -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_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +class LatlonReduced : public Gen { +public: + LatlonReduced() { + class_name_ = "latlon_reduced"; + } + Nearest* create() override { return new LatlonReduced(); } + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +private: + double* lats_ = nullptr; + int lats_count_ = 0; + double* lons_ = nullptr; + double* distances_ = nullptr; + size_t* k_ = nullptr; + size_t* j_ = nullptr; + const char* Nj_ = nullptr; + const char* pl_ = nullptr; + const char* lonFirst_ = nullptr; + const char* lonLast_ = nullptr; + + + int find_global(grib_handle*, + double, double, unsigned long, + double*, double*, double*, + double*, int*, size_t*); +}; + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_mercator.cc b/src/geo_nearest/grib_nearest_class_mercator.cc new file mode 100644 index 000000000..e494af5b2 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_mercator.cc @@ -0,0 +1,65 @@ +/* + * (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_nearest_class_mercator.h" + +eccodes::geo_nearest::Mercator _grib_nearest_mercator{}; +eccodes::geo_nearest::Mercator* grib_nearest_mercator = &_grib_nearest_mercator; + +namespace eccodes::geo_nearest { + +int Mercator::init(grib_handle* h, grib_arguments* args) +{ + int ret = GRIB_SUCCESS; + if ((ret = Gen::init(h, args) != GRIB_SUCCESS)) + return ret; + + Ni_ = grib_arguments_get_name(h, args, cargs_++); + Nj_ = grib_arguments_get_name(h, args, cargs_++); + lats_ = lons_ = distances_ = NULL; + lats_count_ = lons_count_ = 0; + i_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + j_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + return ret; +} + +int Mercator::find(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, + double* values, double* distances, int* indexes, size_t* len) +{ + return grib_nearest_find_generic( + h, inlat, inlon, flags, /* inputs */ + + values_key_, /* outputs to set the 'self' object */ + &(lats_), + &(lats_count_), + &(lons_), + &(lons_count_), + &(distances_), + + outlats, outlons, /* outputs of the find function */ + values, distances, indexes, len); +} + +int Mercator::destroy() +{ + grib_context* c = grib_context_get_default(); + + if (lats_) grib_context_free(c, lats_); + if (lons_) grib_context_free(c, lons_); + if (i_) grib_context_free(c, i_); + if (j_) grib_context_free(c, j_); + if (k_) grib_context_free(c, k_); + if (distances_) grib_context_free(c, distances_); + return GRIB_SUCCESS; +} + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_mercator.h b/src/geo_nearest/grib_nearest_class_mercator.h new file mode 100644 index 000000000..5bdb50a2b --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_mercator.h @@ -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_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +class Mercator : public Gen { +public: + Mercator() { + class_name_ = "mercator"; + } + Nearest* create() override { return new Mercator(); }; + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +private: + double* lats_ = nullptr; + int lats_count_ = 0; + double* lons_ = nullptr; + int lons_count_ = 0; + double* distances_ = nullptr; + int* k_ = nullptr; + int* i_ = nullptr; + int* j_ = nullptr; + const char* Ni_ = nullptr; + const char* Nj_ = nullptr; +}; + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_polar_stereographic.cc b/src/geo_nearest/grib_nearest_class_polar_stereographic.cc new file mode 100644 index 000000000..5f1da4b35 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_polar_stereographic.cc @@ -0,0 +1,63 @@ +/* + * (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_nearest_class_polar_stereographic.h" + +eccodes::geo_nearest::PolarStereographic _grib_nearest_polar_stereographic{}; +eccodes::geo_nearest::PolarStereographic* grib_nearest_polar_stereographic = &_grib_nearest_polar_stereographic; + +namespace eccodes::geo_nearest { + +int PolarStereographic::init(grib_handle* h, grib_arguments* args) +{ + int ret = GRIB_SUCCESS; + if ((ret = Gen::init(h, args) != GRIB_SUCCESS)) + return ret; + + Ni_ = grib_arguments_get_name(h, args, cargs_++); + Nj_ = grib_arguments_get_name(h, args, cargs_++); + i_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + j_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + return ret; +} + +int PolarStereographic::find(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, + double* values, double* distances, int* indexes, size_t* len) +{ + return grib_nearest_find_generic( + h, inlat, inlon, flags, /* inputs */ + + values_key_, /* outputs to set the 'self' object */ + &(lats_), + &(lats_count_), + &(lons_), + &(lons_count_), + &(distances_), + + outlats, outlons, /* outputs of the find function */ + values, distances, indexes, len); +} + +int PolarStereographic::destroy() +{ + grib_context* c = grib_context_get_default(); + + if (lats_) grib_context_free(c, lats_); + if (lons_) grib_context_free(c, lons_); + if (i_) grib_context_free(c, i_); + if (j_) grib_context_free(c, j_); + if (k_) grib_context_free(c, k_); + if (distances_) grib_context_free(c, distances_); + return GRIB_SUCCESS; +} + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_polar_stereographic.h b/src/geo_nearest/grib_nearest_class_polar_stereographic.h new file mode 100644 index 000000000..0af9ebb4c --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_polar_stereographic.h @@ -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_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +class PolarStereographic : public Gen { +public: + PolarStereographic() { + class_name_ = "polar_stereographic"; + } + Nearest* create() override { return new PolarStereographic(); } + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +private: + double* lats_ = nullptr; + int lats_count_ = 0; + double* lons_ = nullptr; + int lons_count_ = 0; + double* distances_ = nullptr; + int* k_ = nullptr; + int* i_ = nullptr; + int* j_ = nullptr; + const char* Ni_ = nullptr; + const char* Nj_ = nullptr; +}; + +} // namespace eccodes::geo_nearest diff --git a/src/grib_nearest_class_reduced.cc b/src/geo_nearest/grib_nearest_class_reduced.cc similarity index 51% rename from src/grib_nearest_class_reduced.cc rename to src/geo_nearest/grib_nearest_class_reduced.cc index ffc8eaad8..f5c976282 100644 --- a/src/grib_nearest_class_reduced.cc +++ b/src/geo_nearest/grib_nearest_class_reduced.cc @@ -8,119 +8,43 @@ * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. */ -#include "grib_api_internal.h" +#include "grib_nearest_class_reduced.h" -/* - This is used by make_class.pl +eccodes::geo_nearest::Reduced _grib_nearest_reduced{}; +eccodes::geo_nearest::Reduced* grib_nearest_reduced = &_grib_nearest_reduced; - START_CLASS_DEF - CLASS = nearest - SUPER = grib_nearest_class_gen - IMPLEMENTS = init;destroy;find - MEMBERS = double* lats - MEMBERS = int lats_count - MEMBERS = double* lons - MEMBERS = double* distances - MEMBERS = size_t* k - MEMBERS = size_t* j - MEMBERS = const char* Nj - MEMBERS = const char* pl - MEMBERS = long global - MEMBERS = double lon_first - MEMBERS = double lon_last - MEMBERS = int legacy - MEMBERS = int rotated - 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_reduced{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; - /* Members defined in reduced */ - double* lats; - int lats_count; - double* lons; - double* distances; - size_t* k; - size_t* j; - const char* Nj; - const char* pl; - long global; - double lon_first; - double lon_last; - int legacy; - int rotated; -} grib_nearest_reduced; - -extern grib_nearest_class* grib_nearest_class_gen; - -static grib_nearest_class _grib_nearest_class_reduced = { - &grib_nearest_class_gen, /* super */ - "reduced", /* name */ - sizeof(grib_nearest_reduced), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_reduced = &_grib_nearest_class_reduced; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ +namespace eccodes::geo_nearest { #define NUM_NEIGHBOURS 4 -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) +int Reduced::init(grib_handle* h, grib_arguments* args) { - grib_nearest_reduced* self = (grib_nearest_reduced*)nearest; - self->Nj = grib_arguments_get_name(h, args, self->cargs++); - self->pl = grib_arguments_get_name(h, args, self->cargs++); - self->j = (size_t*)grib_context_malloc(h->context, 2 * sizeof(size_t)); - self->legacy = -1; - self->rotated = -1; - if (!self->j) + int ret = GRIB_SUCCESS; + if ((ret = Gen::init(h, args) != GRIB_SUCCESS)) + return ret; + + Nj_ = grib_arguments_get_name(h, args, cargs_++); + pl_ = grib_arguments_get_name(h, args, cargs_++); + j_ = (size_t*)grib_context_malloc(h->context, 2 * sizeof(size_t)); + legacy_ = -1; + rotated_ = -1; + if (!j_) return GRIB_OUT_OF_MEMORY; - self->k = (size_t*)grib_context_malloc(h->context, NUM_NEIGHBOURS * sizeof(size_t)); - if (!self->k) + k_ = (size_t*)grib_context_malloc(h->context, NUM_NEIGHBOURS * sizeof(size_t)); + if (!k_) return GRIB_OUT_OF_MEMORY; - grib_get_long(h, "global", &self->global); - if (!self->global) { + grib_get_long(h, "global", &global_); + if (!global_) { int err; /*TODO longitudeOfFirstGridPointInDegrees from the def file*/ - if ((err = grib_get_double(h, "longitudeOfFirstGridPointInDegrees", &self->lon_first)) != GRIB_SUCCESS) { + if ((err = grib_get_double(h, "longitudeOfFirstGridPointInDegrees", &lon_first_)) != GRIB_SUCCESS) { grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_reduced: Unable to get longitudeOfFirstGridPointInDegrees %s\n", grib_get_error_message(err)); return err; } /*TODO longitudeOfLastGridPointInDegrees from the def file*/ - if ((err = grib_get_double(h, "longitudeOfLastGridPointInDegrees", &self->lon_last)) != GRIB_SUCCESS) { + if ((err = grib_get_double(h, "longitudeOfLastGridPointInDegrees", &lon_last_)) != GRIB_SUCCESS) { grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_reduced: Unable to get longitudeOfLastGridPointInDegrees %s\n", grib_get_error_message(err)); @@ -128,16 +52,11 @@ static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) } } - return 0; + return ret; } typedef void (*get_reduced_row_proc)(long pl, double lon_first, double lon_last, long* npoints, long* ilon_first, long* ilon_last); -static int find_global(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, double* values, - double* distances, int* indexes, size_t* len); - static int is_legacy(grib_handle* h, int* legacy) { int err = 0; @@ -159,21 +78,20 @@ static int is_rotated(grib_handle* h, int* rotated) return GRIB_SUCCESS; } -static int find(grib_nearest* nearest, grib_handle* h, +int Reduced::find(grib_handle* h, double inlat, double inlon, unsigned long flags, double* outlats, double* outlons, double* values, double* distances, int* indexes, size_t* len) { int err = 0; - grib_nearest_reduced* self = (grib_nearest_reduced*)nearest; - if (self->rotated == -1 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { - err = is_rotated(h, &(self->rotated)); + if (rotated_ == -1 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { + err = is_rotated(h, &(rotated_)); if (err) return err; } - if (self->global && self->rotated == 0) { - err = find_global(nearest, h, + if (global_ && rotated_ == 0) { + err = find_global( h, inlat, inlon, flags, outlats, outlons, values, distances, indexes, len); @@ -186,13 +104,13 @@ static int find(grib_nearest* nearest, grib_handle* h, int lons_count = 0; /*dummy*/ err = grib_nearest_find_generic( - nearest, h, inlat, inlon, flags, - self->values_key, - &(self->lats), - &(self->lats_count), - &(self->lons), + h, inlat, inlon, flags, + values_key_, + &(lats_), + &(lats_count_), + &(lons_), &(lons_count), - &(self->distances), + &(distances_), outlats, outlons, values, distances, indexes, len); } @@ -200,12 +118,11 @@ static int find(grib_nearest* nearest, grib_handle* h, } /* Old implementation in src/deprecated/grib_nearest_class_reduced.old */ -static int find_global(grib_nearest* nearest, grib_handle* h, +int Reduced::find_global(grib_handle* h, double inlat, double inlon, unsigned long flags, double* outlats, double* outlons, double* values, double* distances, int* indexes, size_t* len) { - grib_nearest_reduced* self = (grib_nearest_reduced*)nearest; int err = 0, kk = 0, ii = 0; size_t jj = 0; long* pla = NULL; @@ -217,17 +134,17 @@ static int find_global(grib_nearest* nearest, grib_handle* h, int ilat = 0, ilon = 0; get_reduced_row_proc get_reduced_row_func = &grib_get_reduced_row; - if (self->legacy == -1 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { - err = is_legacy(h, &(self->legacy)); + if (legacy_ == -1 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { + err = is_legacy(h, &(legacy_)); if (err) return err; } - if (self->legacy == 1) { + if (legacy_ == 1) { get_reduced_row_func = &grib_get_reduced_row_legacy; } - if ((err = grib_get_size(h, self->values_key, &nvalues)) != GRIB_SUCCESS) + if ((err = grib_get_size(h, values_key_, &nvalues)) != GRIB_SUCCESS) return err; - nearest->values_count = nvalues; + values_count_ = nvalues; if ((err = grib_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS) return err; @@ -235,31 +152,31 @@ static int find_global(grib_nearest* nearest, grib_handle* h, /* Compute lat/lon info, create iterator etc if it's the 1st time or different grid. * This is for performance: if the grid has not changed, we only do this once * and reuse for other messages */ - if (!nearest->h || (flags & GRIB_NEAREST_SAME_GRID) == 0) { + if (!h_ || (flags & GRIB_NEAREST_SAME_GRID) == 0) { double olat = 1.e10; long n = 0; ilat = 0; ilon = 0; - if (grib_is_missing(h, self->Nj, &err)) { - grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->Nj); + if (grib_is_missing(h, Nj_, &err)) { + grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Nj_); return err ? err : GRIB_GEOCALCULUS_PROBLEM; } - if ((err = grib_get_long(h, self->Nj, &n)) != GRIB_SUCCESS) + if ((err = grib_get_long(h, Nj_, &n)) != GRIB_SUCCESS) return err; - self->lats_count = n; + lats_count_ = n; - if (self->lats) - grib_context_free(h->context, self->lats); - self->lats = (double*)grib_context_malloc(h->context, self->lats_count * sizeof(double)); - if (!self->lats) + if (lats_) + grib_context_free(h->context, lats_); + lats_ = (double*)grib_context_malloc(h->context, lats_count_ * sizeof(double)); + if (!lats_) return GRIB_OUT_OF_MEMORY; - if (self->lons) - grib_context_free(h->context, self->lons); - self->lons = (double*)grib_context_malloc(h->context, nearest->values_count * sizeof(double)); - if (!self->lons) + if (lons_) + grib_context_free(h->context, lons_); + lons_ = (double*)grib_context_malloc(h->context, values_count_ * sizeof(double)); + if (!lons_) return GRIB_OUT_OF_MEMORY; iter = grib_iterator_new(h, GRIB_GEOITERATOR_NO_VALUES, &err); @@ -268,29 +185,29 @@ static int find_global(grib_nearest* nearest, grib_handle* h, return err; } while (grib_iterator_next(iter, &lat, &lon, NULL)) { - if (ilat < self->lats_count && olat != lat) { - self->lats[ilat++] = lat; + if (ilat < lats_count_ && olat != lat) { + lats_[ilat++] = lat; olat = lat; } while (lon > 360) lon -= 360; - if (!self->global) { /* ECC-756 */ - if (self->legacy == 0) /*TODO*/ + if (!global_) { /* ECC-756 */ + if (legacy_ == 0) /*TODO*/ if (lon > 180 && lon < 360) lon -= 360; } - DEBUG_ASSERT_ACCESS(self->lons, (long)ilon, (long)nearest->values_count); - self->lons[ilon++] = lon; + DEBUG_ASSERT_ACCESS(lons_, (long)ilon, (long)nearest->values_count); + lons_[ilon++] = lon; } - self->lats_count = ilat; + lats_count_ = ilat; grib_iterator_delete(iter); } - nearest->h = h; + h_ = h; /* Compute distances if it's the 1st time or different point or different grid. * This is for performance: if the grid and the input point have not changed * we only do this once and reuse for other messages */ - if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT) == 0 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { + if (!distances_ || (flags & GRIB_NEAREST_SAME_POINT) == 0 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { double* lons = NULL; int nlon = 0; size_t plsize = 0; @@ -298,40 +215,40 @@ static int find_global(grib_nearest* nearest, grib_handle* h, int nearest_lons_found = 0; long row_count, ilon_first, ilon_last; - if (self->global) { + if (global_) { inlon = normalise_longitude_in_degrees(inlon); } else { /* TODO: Experimental */ - if (self->legacy == 0) + if (legacy_ == 0) if (inlon > 180 && inlon < 360) inlon -= 360; } - ilat = self->lats_count; - if (self->lats[ilat - 1] > self->lats[0]) { - if (inlat < self->lats[0] || inlat > self->lats[ilat - 1]) + ilat = lats_count_; + if (lats_[ilat - 1] > lats_[0]) { + if (inlat < lats_[0] || inlat > lats_[ilat - 1]) return GRIB_OUT_OF_AREA; } else { - if (inlat > self->lats[0] || inlat < self->lats[ilat - 1]) + if (inlat > lats_[0] || inlat < lats_[ilat - 1]) return GRIB_OUT_OF_AREA; } - if (!self->distances) - self->distances = (double*)grib_context_malloc(h->context, NUM_NEIGHBOURS * sizeof(double)); - if (!self->distances) + if (!distances_) + distances_ = (double*)grib_context_malloc(h->context, NUM_NEIGHBOURS * sizeof(double)); + if (!distances_) return GRIB_OUT_OF_MEMORY; - grib_binary_search(self->lats, ilat - 1, inlat, &(self->j[0]), &(self->j[1])); + grib_binary_search(lats_, ilat - 1, inlat, &(j_[0]), &(j_[1])); - plsize = self->lats_count; - if ((err = grib_get_size(h, self->pl, &plsize)) != GRIB_SUCCESS) + plsize = lats_count_; + if ((err = grib_get_size(h, pl_, &plsize)) != GRIB_SUCCESS) return err; pla = (long*)grib_context_malloc(h->context, plsize * sizeof(long)); if (!pla) return GRIB_OUT_OF_MEMORY; - if ((err = grib_get_long_array(h, self->pl, pla, &plsize)) != GRIB_SUCCESS) + if ((err = grib_get_long_array(h, pl_, pla, &plsize)) != GRIB_SUCCESS) return err; pl = pla; @@ -340,27 +257,27 @@ static int find_global(grib_nearest* nearest, grib_handle* h, } nlon = 0; - if (self->global) { - for (jj = 0; jj < self->j[0]; jj++) + if (global_) { + for (jj = 0; jj < j_[0]; jj++) nlon += pl[jj]; - nplm1 = pl[self->j[0]] - 1; + nplm1 = pl[j_[0]] - 1; } else { nlon = 0; - for (jj = 0; jj < self->j[0]; jj++) { + for (jj = 0; jj < j_[0]; jj++) { row_count = 0; ilon_first = 0; ilon_last = 0; - get_reduced_row_func(pl[jj], self->lon_first, self->lon_last, &row_count, &ilon_first, &ilon_last); + get_reduced_row_func(pl[jj], lon_first_, lon_last_, &row_count, &ilon_first, &ilon_last); nlon += row_count; } row_count = 0; ilon_first = 0; ilon_last = 0; - get_reduced_row_func(pl[self->j[0]], self->lon_first, self->lon_last, &row_count, &ilon_first, &ilon_last); + get_reduced_row_func(pl[j_[0]], lon_first_, lon_last_, &row_count, &ilon_first, &ilon_last); nplm1 = row_count - 1; } - lons = self->lons + nlon; + lons = lons_ + nlon; nearest_lons_found = 0; /* ECC-756: The comparisons of longitudes here depends on the longitude values @@ -371,8 +288,8 @@ static int find_global(grib_nearest* nearest, grib_handle* h, if (lons[nplm1] > lons[0]) { if (inlon < lons[0] || inlon > lons[nplm1]) { if (lons[nplm1] - lons[0] - 360 <= lons[nplm1] - lons[nplm1 - 1]) { - self->k[0] = 0; - self->k[1] = nplm1; + k_[0] = 0; + k_[1] = nplm1; nearest_lons_found = 1; } else @@ -382,8 +299,8 @@ static int find_global(grib_nearest* nearest, grib_handle* h, else { if (inlon > lons[0] || inlon < lons[nplm1]) { if (lons[0] - lons[nplm1] - 360 <= lons[0] - lons[1]) { - self->k[0] = 0; - self->k[1] = nplm1; + k_[0] = 0; + k_[1] = nplm1; nearest_lons_found = 1; } else @@ -392,51 +309,51 @@ static int find_global(grib_nearest* nearest, grib_handle* h, } if (!nearest_lons_found) { - if (!self->global) { + if (!global_) { row_count = 0; ilon_first = 0; ilon_last = 0; - get_reduced_row_func(pl[self->j[0]], self->lon_first, self->lon_last, &row_count, &ilon_first, &ilon_last); + get_reduced_row_func(pl[j_[0]], lon_first_, lon_last_, &row_count, &ilon_first, &ilon_last); } else { - row_count = pl[self->j[0]]; + row_count = pl[j_[0]]; } grib_binary_search(lons, row_count - 1, inlon, - &(self->k[0]), &(self->k[1])); + &(k_[0]), &(k_[1])); } - self->k[0] += nlon; - self->k[1] += nlon; + k_[0] += nlon; + k_[1] += nlon; nlon = 0; - if (self->global) { - for (jj = 0; jj < self->j[1]; jj++) + if (global_) { + for (jj = 0; jj < j_[1]; jj++) nlon += pl[jj]; - nplm1 = pl[self->j[1]] - 1; + nplm1 = pl[j_[1]] - 1; } else { - for (jj = 0; jj < self->j[1]; jj++) { + for (jj = 0; jj < j_[1]; jj++) { row_count = 0; ilon_first = 0; ilon_last = 0; - get_reduced_row_func(pl[jj], self->lon_first, self->lon_last, &row_count, &ilon_first, &ilon_last); + get_reduced_row_func(pl[jj], lon_first_, lon_last_, &row_count, &ilon_first, &ilon_last); nlon += row_count; } row_count = 0; ilon_first = 0; ilon_last = 0; - get_reduced_row_func(pl[self->j[1]], self->lon_first, self->lon_last, &nplm1, &ilon_first, &ilon_last); + get_reduced_row_func(pl[j_[1]], lon_first_, lon_last_, &nplm1, &ilon_first, &ilon_last); nplm1--; } - lons = self->lons + nlon; + lons = lons_ + nlon; nearest_lons_found = 0; if (lons[nplm1] > lons[0]) { if (inlon < lons[0] || inlon > lons[nplm1]) { if (lons[nplm1] - lons[0] - 360 <= lons[nplm1] - lons[nplm1 - 1]) { - self->k[2] = 0; - self->k[3] = nplm1; + k_[2] = 0; + k_[3] = nplm1; nearest_lons_found = 1; } else @@ -447,8 +364,8 @@ static int find_global(grib_nearest* nearest, grib_handle* h, if (inlon > lons[0] || inlon < lons[nplm1]) { if (lons[0] - lons[nplm1] - 360 <= lons[0] - lons[1]) { - self->k[2] = 0; - self->k[3] = nplm1; + k_[2] = 0; + k_[3] = nplm1; nearest_lons_found = 1; } else @@ -457,28 +374,28 @@ static int find_global(grib_nearest* nearest, grib_handle* h, } if (!nearest_lons_found) { - if (!self->global) { + if (!global_) { row_count = 0; ilon_first = 0; ilon_last = 0; - get_reduced_row_func(pl[self->j[1]], self->lon_first, self->lon_last, &row_count, &ilon_first, &ilon_last); + get_reduced_row_func(pl[j_[1]], lon_first_, lon_last_, &row_count, &ilon_first, &ilon_last); } else { - row_count = pl[self->j[1]]; + row_count = pl[j_[1]]; } grib_binary_search(lons, row_count - 1, inlon, - &(self->k[2]), &(self->k[3])); + &(k_[2]), &(k_[3])); } - self->k[2] += nlon; - self->k[3] += nlon; + k_[2] += nlon; + k_[3] += nlon; kk = 0; for (jj = 0; jj < 2; jj++) { for (ii = 0; ii < 2; ii++) { - self->distances[kk] = geographic_distance_spherical(radiusInKm, inlon, inlat, - self->lons[self->k[kk]], self->lats[self->j[jj]]); + distances_[kk] = geographic_distance_spherical(radiusInKm, inlon, inlat, + lons_[k_[kk]], lats_[j_[jj]]); kk++; } } @@ -490,24 +407,24 @@ static int find_global(grib_nearest* nearest, grib_handle* h, if (values) { /* See ECC-1403 and ECC-499 */ /* Performance: Decode the field once and get all 4 values */ - err = grib_get_double_element_set(h, self->values_key, self->k, NUM_NEIGHBOURS, values); + err = grib_get_double_element_set(h, values_key_, k_, NUM_NEIGHBOURS, values); if (err != GRIB_SUCCESS) return err; } for (jj = 0; jj < 2; jj++) { for (ii = 0; ii < 2; ii++) { - distances[kk] = self->distances[kk]; - outlats[kk] = self->lats[self->j[jj]]; - outlons[kk] = self->lons[self->k[kk]]; + distances[kk] = distances_[kk]; + outlats[kk] = lats_[j_[jj]]; + outlons[kk] = lons_[k_[kk]]; /*if (values) { - * grib_get_double_element_internal(h, self->values_key, self->k[kk], &(values[kk])); + * grib_get_double_element_internal(h, values_key_, k_[kk], &(values[kk])); *} */ - if (self->k[kk] >= INT_MAX) { + if (k_[kk] >= INT_MAX) { /* Current interface uses an 'int' for 'indexes' which is 32bits! We should change this to a 64bit type */ grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_reduced: Unable to compute index. Value too large"); return GRIB_OUT_OF_RANGE; } else { - indexes[kk] = (int)self->k[kk]; + indexes[kk] = (int)k_[kk]; } kk++; } @@ -516,16 +433,17 @@ static int find_global(grib_nearest* nearest, grib_handle* h, return GRIB_SUCCESS; } -static int destroy(grib_nearest* nearest) +int Reduced::destroy() { - grib_nearest_reduced* self = (grib_nearest_reduced*)nearest; grib_context* c = grib_context_get_default(); - if (self->lats) grib_context_free(c, self->lats); - if (self->lons) grib_context_free(c, self->lons); - if (self->j) grib_context_free(c, self->j); - if (self->k) grib_context_free(c, self->k); - if (self->distances) grib_context_free(c, self->distances); + if (lats_) grib_context_free(c, lats_); + if (lons_) grib_context_free(c, lons_); + if (j_) grib_context_free(c, j_); + if (k_) grib_context_free(c, k_); + if (distances_) grib_context_free(c, distances_); return GRIB_SUCCESS; } + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_reduced.h b/src/geo_nearest/grib_nearest_class_reduced.h new file mode 100644 index 000000000..ccf499cf5 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_reduced.h @@ -0,0 +1,48 @@ +/* + * (C) Copyright 2005- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * + * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by + * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. + */ + +#pragma once + +#include "grib_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +class Reduced : public Gen { +public: + Reduced() { + class_name_ = "reduced"; + } + Nearest* create() override { return new Reduced(); } + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +private: + double* lats_ = nullptr; + int lats_count_ = 0; + double* lons_ = nullptr; + double* distances_ = nullptr; + size_t* k_ = nullptr; + size_t* j_ = nullptr; + const char* Nj_ = nullptr; + const char* pl_ = nullptr; + long global_ = 0.0; + double lon_first_ = 0.0; + double lon_last_ = 0.0; + int legacy_ = 0; + int rotated_ = 0; + + int find_global(grib_handle*, + double, double, unsigned long, + double*, double*, double*, + double*, int*, size_t*); +}; + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_regular.cc b/src/geo_nearest/grib_nearest_class_regular.cc new file mode 100644 index 000000000..5f04c1839 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_regular.cc @@ -0,0 +1,302 @@ +/* + * (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_nearest_class_regular.h" + +eccodes::geo_nearest::Regular _grib_nearest_regular{}; +eccodes::geo_nearest::Regular* grib_nearest_regular = &_grib_nearest_regular; + +namespace eccodes::geo_nearest { + +#define NUM_NEIGHBOURS 4 + +int Regular::init(grib_handle* h, grib_arguments* args) +{ + int ret = GRIB_SUCCESS; + if ((ret = Gen::init(h, args) != GRIB_SUCCESS)) + return ret; + + Ni_ = grib_arguments_get_name(h, args, cargs_++); + Nj_ = grib_arguments_get_name(h, args, cargs_++); + i_ = (size_t*)grib_context_malloc(h->context, 2 * sizeof(size_t)); + j_ = (size_t*)grib_context_malloc(h->context, 2 * sizeof(size_t)); + return ret; +} + +static bool is_rotated_grid(grib_handle* h) +{ + long is_rotated = 0; + int err = grib_get_long(h, "isRotatedGrid", &is_rotated); + if (!err && is_rotated) + return true; + return false; +} + +// Old implementation in +// src/deprecated/grib_nearest_class_regular.cc +// +int Regular::find(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, + double* values, double* distances, int* indexes, size_t* len) +{ + int ret = 0, kk = 0, ii = 0, jj = 0; + size_t nvalues = 0; + double radiusInKm; + + grib_iterator* iter = NULL; + double lat = 0, lon = 0; + const bool is_rotated = is_rotated_grid(h); + double angleOfRotation = 0, southPoleLat = 0, southPoleLon = 0; + grib_context* c = h->context; + + while (inlon < 0) + inlon += 360; + while (inlon > 360) + inlon -= 360; + + if ((ret = grib_get_size(h, values_key_, &nvalues)) != GRIB_SUCCESS) + return ret; + values_count_ = nvalues; + + if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS) + return ret; + + /* Compute lat/lon info, create iterator etc if it's the 1st time or different grid. + * This is for performance: if the grid has not changed, we only do this once + * and reuse for other messages */ + if (!h_ || (flags & GRIB_NEAREST_SAME_GRID) == 0) { + double olat = 1.e10, olon = 1.e10; + int ilat = 0, ilon = 0; + long n = 0; + + if (grib_is_missing(h, Ni_, &ret)) { + grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Ni_); + return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; + } + + if (grib_is_missing(h, Nj_, &ret)) { + grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Nj_); + return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; + } + + /* ECC-600: Support for rotated grids + * First: rotate the input point + * Then: run the lat/lon iterator over the rotated grid (disableUnrotate) + * Finally: unrotate the resulting point + */ + if (is_rotated) { + double new_lat = 0, new_lon = 0; + ret = grib_get_double_internal(h, "angleOfRotation", &angleOfRotation); + if (ret) + return ret; + ret = grib_get_double_internal(h, "latitudeOfSouthernPoleInDegrees", &southPoleLat); + if (ret) + return ret; + ret = grib_get_double_internal(h, "longitudeOfSouthernPoleInDegrees", &southPoleLon); + if (ret) + return ret; + ret = grib_set_long(h, "iteratorDisableUnrotate", 1); + if (ret) + return ret; + /* Rotate the inlat, inlon */ + rotate(inlat, inlon, angleOfRotation, southPoleLat, southPoleLon, &new_lat, &new_lon); + inlat = new_lat; + inlon = new_lon; + /*if(h->context->debug) printf("nearest find: rotated grid: new point=(%g,%g)\n",new_lat,new_lon);*/ + } + + if ((ret = grib_get_long(h, Ni_, &n)) != GRIB_SUCCESS) + return ret; + lons_count_ = n; + + if ((ret = grib_get_long(h, Nj_, &n)) != GRIB_SUCCESS) + return ret; + lats_count_ = n; + + if (lats_) + grib_context_free(c, lats_); + lats_ = (double*)grib_context_malloc(c, lats_count_ * sizeof(double)); + if (!lats_) + return GRIB_OUT_OF_MEMORY; + + if (lons_) + grib_context_free(c, lons_); + lons_ = (double*)grib_context_malloc(c, lons_count_ * sizeof(double)); + if (!lons_) + return GRIB_OUT_OF_MEMORY; + + iter = grib_iterator_new(h, GRIB_GEOITERATOR_NO_VALUES, &ret); + if (ret != GRIB_SUCCESS) { + grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_regular: Unable to create lat/lon iterator"); + return ret; + } + while (grib_iterator_next(iter, &lat, &lon, NULL)) { + if (ilat < lats_count_ && olat != lat) { + lats_[ilat++] = lat; + olat = lat; + } + if (ilon < lons_count_ && olon != lon) { + lons_[ilon++] = lon; + olon = lon; + } + } + grib_iterator_delete(iter); + } + h_ = h; + + /* Compute distances if it's the 1st time or different point or different grid. + * This is for performance: if the grid and the input point have not changed + * we only do this once and reuse for other messages */ + if (!distances_ || (flags & GRIB_NEAREST_SAME_POINT) == 0 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { + int nearest_lons_found = 0; + + if (lats_[lats_count_ - 1] > lats_[0]) { + if (inlat < lats_[0] || inlat > lats_[lats_count_ - 1]) + return GRIB_OUT_OF_AREA; + } + else { + if (inlat > lats_[0] || inlat < lats_[lats_count_ - 1]) + return GRIB_OUT_OF_AREA; + } + + if (lons_[lons_count_ - 1] > lons_[0]) { + if (inlon < lons_[0] || inlon > lons_[lons_count_ - 1]) { + /* try to scale*/ + if (inlon > 0) + inlon -= 360; + else + inlon += 360; + + if (inlon < lons_[0] || inlon > lons_[lons_count_ - 1]) { + if (lons_[0] + 360 - lons_[lons_count_ - 1] <= + lons_[1] - lons_[0]) { + /*it's a global field in longitude*/ + i_[0] = 0; + i_[1] = lons_count_ - 1; + nearest_lons_found = 1; + } + else + return GRIB_OUT_OF_AREA; + } + } + } + else { + if (inlon > lons_[0] || inlon < lons_[lons_count_ - 1]) { + /* try to scale*/ + if (inlon > 0) + inlon -= 360; + else + inlon += 360; + if (lons_[0] - lons_[lons_count_ - 1] - 360 <= + lons_[0] - lons_[1]) { + /*it's a global field in longitude*/ + i_[0] = 0; + i_[1] = lons_count_ - 1; + nearest_lons_found = 1; + } + else if (inlon > lons_[0] || inlon < lons_[lons_count_ - 1]) + return GRIB_OUT_OF_AREA; + } + } + + grib_binary_search(lats_, lats_count_ - 1, inlat, + &(j_[0]), &(j_[1])); + + if (!nearest_lons_found) + grib_binary_search(lons_, lons_count_ - 1, inlon, + &(i_[0]), &(i_[1])); + + if (!distances_) + distances_ = (double*)grib_context_malloc(c, NUM_NEIGHBOURS * sizeof(double)); + if (!k_) + k_ = (size_t*)grib_context_malloc(c, NUM_NEIGHBOURS * sizeof(size_t)); + kk = 0; + for (jj = 0; jj < 2; jj++) { + for (ii = 0; ii < 2; ii++) { + k_[kk] = i_[ii] + lons_count_ * j_[jj]; + distances_[kk] = geographic_distance_spherical(radiusInKm, inlon, inlat, + lons_[i_[ii]], lats_[j_[jj]]); + kk++; + } + } + } + + kk = 0; + + /* + * Brute force algorithm: + * First unpack all the values into an array. Then when we need the 4 points + * we just index into this array so no need to call grib_get_double_element_internal + * + * if (nearest->values) grib_context_free(c,nearest->values); + * nearest->values = grib_context_malloc(h->context,nvalues*sizeof(double)); + * if (!nearest->values) return GRIB_OUT_OF_MEMORY; + * ret = grib_get_double_array(h, values_key_, nearest->values ,&nvalues); + * if (ret) return ret; + */ + + if (values) { + /* See ECC-1403 and ECC-499 */ + /* Performance: Decode the field once and get all 4 values */ + if ((ret = grib_get_double_element_set(h, values_key_, k_, NUM_NEIGHBOURS, values)) != GRIB_SUCCESS) + return ret; + } + + for (jj = 0; jj < 2; jj++) { + for (ii = 0; ii < 2; ii++) { + distances[kk] = distances_[kk]; + outlats[kk] = lats_[j_[jj]]; + outlons[kk] = lons_[i_[ii]]; + if (is_rotated) { + /* Unrotate resulting lat/lon */ + double new_lat = 0, new_lon = 0; + unrotate(outlats[kk], outlons[kk], angleOfRotation, southPoleLat, southPoleLon, &new_lat, &new_lon); + outlats[kk] = new_lat; + outlons[kk] = new_lon; + } + /* See ECC-1403 and ECC-499 + * if (values) { + * grib_get_double_element_internal(h, values_key_, k_[kk], &(values[kk])); + *} + */ + /* Using the brute force approach described above */ + /* Assert(k_[kk] < nvalues); */ + /* values[kk]=nearest->values[k_[kk]]; */ + + if (k_[kk] >= INT_MAX) { + /* Current interface uses an 'int' for 'indexes' which is 32bits! We should change this to a 64bit type */ + grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_regular: Unable to compute index. Value too large"); + return GRIB_OUT_OF_RANGE; + } else { + indexes[kk] = (int)k_[kk]; + } + kk++; + } + } + + return GRIB_SUCCESS; +} + +int Regular::destroy() +{ + grib_context* c = grib_context_get_default(); + + if (lats_) grib_context_free(c, lats_); + if (lons_) grib_context_free(c, lons_); + if (i_) grib_context_free(c, i_); + if (j_) grib_context_free(c, j_); + if (k_) grib_context_free(c, k_); + if (distances_) grib_context_free(c, distances_); + + return Gen::destroy(); +} + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_regular.h b/src/geo_nearest/grib_nearest_class_regular.h new file mode 100644 index 000000000..0113100da --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_regular.h @@ -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_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +class Regular : public Gen { +public: + Regular() { + class_name_ = "regular"; + } + Nearest* create() override { return new Regular(); } + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +private: + double* lats_ = nullptr; + int lats_count_ = 0; + double* lons_ = nullptr; + int lons_count_ = 0; + double* distances_ = nullptr; + size_t* k_ = nullptr; + size_t* i_ = nullptr; + size_t* j_ = nullptr; + const char* Ni_ = nullptr; + const char* Nj_ = nullptr; +}; + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_space_view.cc b/src/geo_nearest/grib_nearest_class_space_view.cc new file mode 100644 index 000000000..59d4edec5 --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_space_view.cc @@ -0,0 +1,65 @@ +/* + * (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_nearest_class_space_view.h" + +eccodes::geo_nearest::SpaceView _grib_nearest_space_view{}; +eccodes::geo_nearest::SpaceView* grib_nearest_space_view = &_grib_nearest_space_view; + +namespace eccodes::geo_nearest { + +int SpaceView::init(grib_handle* h, grib_arguments* args) +{ + int ret = GRIB_SUCCESS; + if ((ret = Gen::init(h, args) != GRIB_SUCCESS)) + return ret; + + Ni_ = grib_arguments_get_name(h, args, cargs_++); + Nj_ = grib_arguments_get_name(h, args, cargs_++); + i_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + j_ = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); + + return ret; +} + +int SpaceView::find(grib_handle* h, + double inlat, double inlon, unsigned long flags, + double* outlats, double* outlons, + double* values, double* distances, int* indexes, size_t* len) +{ + return grib_nearest_find_generic( + h, inlat, inlon, flags, /* inputs */ + + values_key_, /* outputs to set the 'self' object */ + &(lats_), + &(lats_count_), + &(lons_), + &(lons_count_), + &(distances_), + + outlats, outlons, /* outputs of the find function */ + values, distances, indexes, len); +} + +int SpaceView::destroy() +{ + grib_context* c = grib_context_get_default(); + + if (lats_) grib_context_free(c, lats_); + if (lons_) grib_context_free(c, lons_); + if (i_) grib_context_free(c, i_); + if (j_) grib_context_free(c, j_); + if (k_) grib_context_free(c, k_); + if (distances_) grib_context_free(c, distances_); + + return Gen::destroy(); +} + +} // namespace eccodes::geo_nearest diff --git a/src/geo_nearest/grib_nearest_class_space_view.h b/src/geo_nearest/grib_nearest_class_space_view.h new file mode 100644 index 000000000..161fb583a --- /dev/null +++ b/src/geo_nearest/grib_nearest_class_space_view.h @@ -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_nearest_class_gen.h" + +namespace eccodes::geo_nearest { + +class SpaceView : public Gen { +public: + SpaceView() { + class_name_ = "space_view"; + } + Nearest* create() override { return new SpaceView(); } + int init(grib_handle*, grib_arguments*) override; + int find(grib_handle*, double, double, unsigned long, double*, double*, double*, double*, int*, size_t*) override; + int destroy() override; + +private: + double* lats_ = nullptr; + int lats_count_ = 0; + double* lons_ = nullptr; + int lons_count_ = 0; + double* distances_ = nullptr; + int* k_ = nullptr; + int* i_ = nullptr; + int* j_ = nullptr; + const char* Ni_ = nullptr; + const char* Nj_ = nullptr; +}; + +} // namespace eccodes::geo_nearest diff --git a/src/grib_api_internal.h b/src/grib_api_internal.h index b093c2e6c..5a1f05f03 100644 --- a/src/grib_api_internal.h +++ b/src/grib_api_internal.h @@ -256,7 +256,6 @@ typedef struct grib_iterator { eccodes::geo_iterator::Iterator* iterator; } grib_iterator; -typedef struct grib_nearest_class grib_nearest_class; typedef struct grib_dumper grib_dumper; typedef struct grib_dumper_class grib_dumper_class; typedef struct grib_dependency grib_dependency; @@ -264,16 +263,6 @@ typedef struct grib_dependency grib_dependency; typedef struct codes_condition codes_condition; /* typedef void (*dynamic_key_proc) (const char*, void*) */ -typedef void (*nearest_init_class_proc)(grib_nearest_class*); -typedef int (*nearest_init_proc)(grib_nearest* i, grib_handle*, grib_arguments*); - -typedef int (*nearest_find_proc)(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, - unsigned long flags, double* outlats, - double* outlons, double* values, - double* distances, int* indexes, size_t* len); -typedef int (*nearest_destroy_proc)(grib_nearest* nearest); - 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); @@ -458,17 +447,13 @@ struct grib_section size_t padding; }; -struct grib_nearest_class -{ - grib_nearest_class** super; - const char* name; - size_t size; - int inited; - nearest_init_class_proc init_class; - nearest_init_proc init; - nearest_destroy_proc destroy; - nearest_find_proc find; -}; +namespace eccodes::geo_nearest { +class Nearest; +} + +typedef struct grib_nearest { + eccodes::geo_nearest::Nearest* nearest; +} grib_nearest; /* --------------- */ typedef int (*dumper_init_proc)(grib_dumper*); @@ -513,15 +498,6 @@ struct grib_dumper_class dumper_footer_proc footer; }; -struct grib_nearest -{ - grib_handle* h; - double* values; - size_t values_count; - grib_nearest_class* cclass; - unsigned long flags; -}; - struct grib_dependency { grib_dependency* next; @@ -1243,13 +1219,14 @@ typedef struct j2k_encode_helper } j2k_encode_helper; -#include "eccodes_prototypes.h" +#include "eccodes_prototypes.h" #ifdef __cplusplus } #include "accessor/grib_accessor.h" #include "accessor/grib_accessors_list.h" #include "geo_iterator/grib_iterator.h" +#include "geo_nearest/grib_nearest.h" #endif #endif diff --git a/src/grib_iterator_class.h b/src/grib_iterator_class.h deleted file mode 100644 index 31aabe985..000000000 --- a/src/grib_iterator_class.h +++ /dev/null @@ -1,13 +0,0 @@ -/* This file is automatically generated by ./make_class.pl, do not edit */ -extern eccodes::geo_iterator::Iterator* grib_iterator_gaussian; -extern eccodes::geo_iterator::Iterator* grib_iterator_gaussian_reduced; -//extern eccodes::geo_iterator::Iterator* grib_iterator_gen; -extern eccodes::geo_iterator::Iterator* grib_iterator_healpix; -extern eccodes::geo_iterator::Iterator* grib_iterator_lambert_azimuthal_equal_area; -extern eccodes::geo_iterator::Iterator* grib_iterator_lambert_conformal; -extern eccodes::geo_iterator::Iterator* grib_iterator_latlon; -extern eccodes::geo_iterator::Iterator* grib_iterator_latlon_reduced; -extern eccodes::geo_iterator::Iterator* grib_iterator_mercator; -extern eccodes::geo_iterator::Iterator* grib_iterator_polar_stereographic; -extern eccodes::geo_iterator::Iterator* grib_iterator_regular; -extern eccodes::geo_iterator::Iterator* grib_iterator_space_view; diff --git a/src/grib_iterator_class.cc b/src/grib_iterator_factory.cc similarity index 83% rename from src/grib_iterator_class.cc rename to src/grib_iterator_factory.cc index fbf208818..578926889 100644 --- a/src/grib_iterator_class.cc +++ b/src/grib_iterator_factory.cc @@ -8,18 +8,10 @@ * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. */ -/*************************************************************************** - * Jean Baptiste Filippi - 01.11.2005 * - * Enrico Fucile * - ***************************************************************************/ - #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 */ -#include "grib_iterator_class.h" - #if GRIB_PTHREADS static pthread_once_t once = PTHREAD_ONCE_INIT; static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER; @@ -54,11 +46,31 @@ struct table_entry eccodes::geo_iterator::Iterator** iterator; }; -static const struct table_entry table[] = { - /* This file is generated by ./make_class.pl */ - #include "grib_iterator_factory.h" -}; +extern eccodes::geo_iterator::Iterator* grib_iterator_gaussian; +extern eccodes::geo_iterator::Iterator* grib_iterator_gaussian_reduced; +extern eccodes::geo_iterator::Iterator* grib_iterator_healpix; +extern eccodes::geo_iterator::Iterator* grib_iterator_lambert_azimuthal_equal_area; +extern eccodes::geo_iterator::Iterator* grib_iterator_lambert_conformal; +extern eccodes::geo_iterator::Iterator* grib_iterator_latlon; +extern eccodes::geo_iterator::Iterator* grib_iterator_latlon_reduced; +extern eccodes::geo_iterator::Iterator* grib_iterator_mercator; +extern eccodes::geo_iterator::Iterator* grib_iterator_polar_stereographic; +extern eccodes::geo_iterator::Iterator* grib_iterator_regular; +extern eccodes::geo_iterator::Iterator* grib_iterator_space_view; +static const struct table_entry table[] = { + { "gaussian", &grib_iterator_gaussian, }, + { "gaussian_reduced", &grib_iterator_gaussian_reduced, }, + { "healpix", &grib_iterator_healpix, }, + { "lambert_azimuthal_equal_area", &grib_iterator_lambert_azimuthal_equal_area, }, + { "lambert_conformal", &grib_iterator_lambert_conformal, }, + { "latlon", &grib_iterator_latlon, }, + { "latlon_reduced", &grib_iterator_latlon_reduced, }, + { "mercator", &grib_iterator_mercator, }, + { "polar_stereographic", &grib_iterator_polar_stereographic, }, + { "regular", &grib_iterator_regular, }, + { "space_view", &grib_iterator_space_view, }, +}; eccodes::geo_iterator::Iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsigned long flags, int* error) { diff --git a/src/grib_iterator_factory.h b/src/grib_iterator_factory.h index ed247bebd..094df5d58 100644 --- a/src/grib_iterator_factory.h +++ b/src/grib_iterator_factory.h @@ -1,13 +1,3 @@ -/* This file is automatically generated by ./make.pl, do not edit */ -{ "gaussian", &grib_iterator_gaussian, }, -{ "gaussian_reduced", &grib_iterator_gaussian_reduced, }, -//{ "gen", &grib_iterator_gen, }, -{ "healpix", &grib_iterator_healpix, }, -{ "lambert_azimuthal_equal_area", &grib_iterator_lambert_azimuthal_equal_area, }, -{ "lambert_conformal", &grib_iterator_lambert_conformal, }, -{ "latlon", &grib_iterator_latlon, }, -{ "latlon_reduced", &grib_iterator_latlon_reduced, }, -{ "mercator", &grib_iterator_mercator, }, -{ "polar_stereographic", &grib_iterator_polar_stereographic, }, -{ "regular", &grib_iterator_regular, }, -{ "space_view", &grib_iterator_space_view, }, +#pragma once + +eccodes::geo_iterator::Iterator* grib_iterator_factory(grib_handle* h, grib_arguments* args, unsigned long flags, int* error); diff --git a/src/grib_nearest_class.cc b/src/grib_nearest_class.cc deleted file mode 100644 index 8128e844e..000000000 --- a/src/grib_nearest_class.cc +++ /dev/null @@ -1,52 +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 file is generated by ./make_class.pl */ -#include "grib_nearest_class.h" - -struct table_entry -{ - const char* type; - grib_nearest_class** cclass; -}; - -static const struct table_entry table[] = { -/* This file is generated by ./make_class.pl */ -#include "grib_nearest_factory.h" -}; - -grib_nearest* grib_nearest_factory(grib_handle* h, grib_arguments* args, int* error) -{ - size_t i = 0, num_table_entries = 0; - *error = GRIB_NOT_IMPLEMENTED; - char* type = (char*)grib_arguments_get_name(h, args, 0); - - num_table_entries = sizeof(table) / sizeof(table[0]); - for (i = 0; i < num_table_entries; i++) { - if (strcmp(type, table[i].type) == 0) { - grib_nearest_class* c = *(table[i].cclass); - grib_nearest* it = (grib_nearest*)grib_context_malloc_clear(h->context, c->size); - it->cclass = c; - *error = grib_nearest_init(it, h, args); - if (*error == GRIB_SUCCESS) - return it; - grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_factory: Error instantiating nearest %s (%s)", - table[i].type, grib_get_error_message(*error)); - grib_nearest_delete(it); - return NULL; - } - } - - grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_factory: Unknown type: %s", type); - - return NULL; -} diff --git a/src/grib_nearest_class.h b/src/grib_nearest_class.h deleted file mode 100644 index a1341806a..000000000 --- a/src/grib_nearest_class.h +++ /dev/null @@ -1,11 +0,0 @@ -/* This file is automatically generated by ./make_class.pl, do not edit */ -extern grib_nearest_class* grib_nearest_class_gen; -extern grib_nearest_class* grib_nearest_class_healpix; -extern grib_nearest_class* grib_nearest_class_lambert_azimuthal_equal_area; -extern grib_nearest_class* grib_nearest_class_lambert_conformal; -extern grib_nearest_class* grib_nearest_class_latlon_reduced; -extern grib_nearest_class* grib_nearest_class_mercator; -extern grib_nearest_class* grib_nearest_class_polar_stereographic; -extern grib_nearest_class* grib_nearest_class_reduced; -extern grib_nearest_class* grib_nearest_class_regular; -extern grib_nearest_class* grib_nearest_class_space_view; diff --git a/src/grib_nearest_class_gen.cc b/src/grib_nearest_class_gen.cc deleted file mode 100644 index c97e6b8b9..000000000 --- a/src/grib_nearest_class_gen.cc +++ /dev/null @@ -1,102 +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 = nearest - IMPLEMENTS = destroy - IMPLEMENTS = find - IMPLEMENTS = init - MEMBERS = const char* values_key - MEMBERS = const char* radius - MEMBERS = int cargs - 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_gen{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; -} grib_nearest_gen; - - -static grib_nearest_class _grib_nearest_class_gen = { - 0, /* super */ - "gen", /* name */ - sizeof(grib_nearest_gen), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_gen = &_grib_nearest_class_gen; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ - -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) -{ - grib_nearest_gen* self = (grib_nearest_gen*)nearest; - int ret = GRIB_SUCCESS; - - self->cargs = 1; - - self->values_key = grib_arguments_get_name(h, args, self->cargs++); - self->radius = grib_arguments_get_name(h, args, self->cargs++); - nearest->values = NULL; - - return ret; -} - -static int destroy(grib_nearest* nearest) -{ - grib_context* c = grib_context_get_default(); - if (nearest->values) - grib_context_free(c, nearest->values); - grib_context_free(c, nearest); - return GRIB_SUCCESS; -} - -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, double* values, - double* distances, int* indexes, size_t* len) -{ - return GRIB_NOT_IMPLEMENTED; -} diff --git a/src/grib_nearest_class_healpix.cc b/src/grib_nearest_class_healpix.cc deleted file mode 100644 index d8618ba34..000000000 --- a/src/grib_nearest_class_healpix.cc +++ /dev/null @@ -1,136 +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 = nearest - SUPER = grib_nearest_class_gen - IMPLEMENTS = init - IMPLEMENTS = find;destroy - MEMBERS = double* lats - MEMBERS = int lats_count - MEMBERS = double* lons - MEMBERS = int lons_count - MEMBERS = double* distances - MEMBERS = int* k - MEMBERS = int* i - MEMBERS = int* j - MEMBERS = const char* Ni - MEMBERS = const char* 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_healpix{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; - /* Members defined in healpix */ - double* lats; - int lats_count; - double* lons; - int lons_count; - double* distances; - int* k; - int* i; - int* j; - const char* Ni; - const char* Nj; -} grib_nearest_healpix; - -extern grib_nearest_class* grib_nearest_class_gen; - -static grib_nearest_class _grib_nearest_class_healpix = { - &grib_nearest_class_gen, /* super */ - "healpix", /* name */ - sizeof(grib_nearest_healpix), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_healpix = &_grib_nearest_class_healpix; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ - -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) -{ - grib_nearest_healpix* self = (grib_nearest_healpix*)nearest; - self->Ni = grib_arguments_get_name(h, args, self->cargs++); - self->Nj = grib_arguments_get_name(h, args, self->cargs++); - self->lats = self->lons = self->distances = NULL; - self->lats_count = self->lons_count = 0; - self->i = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - self->j = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - return GRIB_SUCCESS; -} - -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, - double* values, double* distances, int* indexes, size_t* len) -{ - grib_nearest_healpix* self = (grib_nearest_healpix*)nearest; - return grib_nearest_find_generic( - nearest, h, inlat, inlon, flags, /* inputs */ - - self->values_key, /* outputs to set the 'self' object */ - &(self->lats), - &(self->lats_count), - &(self->lons), - &(self->lons_count), - &(self->distances), - - outlats, outlons, /* outputs of the find function */ - values, distances, indexes, len); -} - -static int destroy(grib_nearest* nearest) -{ - grib_nearest_healpix* self = (grib_nearest_healpix*)nearest; - grib_context* c = grib_context_get_default(); - - if (self->lats) grib_context_free(c, self->lats); - if (self->lons) grib_context_free(c, self->lons); - if (self->i) grib_context_free(c, self->i); - if (self->j) grib_context_free(c, self->j); - if (self->k) grib_context_free(c, self->k); - if (self->distances) grib_context_free(c, self->distances); - return GRIB_SUCCESS; -} diff --git a/src/grib_nearest_class_lambert_azimuthal_equal_area.cc b/src/grib_nearest_class_lambert_azimuthal_equal_area.cc deleted file mode 100644 index 951f3e571..000000000 --- a/src/grib_nearest_class_lambert_azimuthal_equal_area.cc +++ /dev/null @@ -1,134 +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 = nearest - SUPER = grib_nearest_class_gen - IMPLEMENTS = init - IMPLEMENTS = find;destroy - MEMBERS = double* lats - MEMBERS = int lats_count - MEMBERS = double* lons - MEMBERS = int lons_count - MEMBERS = double* distances - MEMBERS = int* k - MEMBERS = int* i - MEMBERS = int* j - MEMBERS = const char* Ni - MEMBERS = const char* 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_lambert_azimuthal_equal_area{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; - /* Members defined in lambert_azimuthal_equal_area */ - double* lats; - int lats_count; - double* lons; - int lons_count; - double* distances; - int* k; - int* i; - int* j; - const char* Ni; - const char* Nj; -} grib_nearest_lambert_azimuthal_equal_area; - -extern grib_nearest_class* grib_nearest_class_gen; - -static grib_nearest_class _grib_nearest_class_lambert_azimuthal_equal_area = { - &grib_nearest_class_gen, /* super */ - "lambert_azimuthal_equal_area", /* name */ - sizeof(grib_nearest_lambert_azimuthal_equal_area), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_lambert_azimuthal_equal_area = &_grib_nearest_class_lambert_azimuthal_equal_area; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ - -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) -{ - grib_nearest_lambert_azimuthal_equal_area* self = (grib_nearest_lambert_azimuthal_equal_area*)nearest; - self->Ni = grib_arguments_get_name(h, args, self->cargs++); - self->Nj = grib_arguments_get_name(h, args, self->cargs++); - self->i = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - self->j = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - return 0; -} - -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, - double* values, double* distances, int* indexes, size_t* len) -{ - grib_nearest_lambert_azimuthal_equal_area* self = (grib_nearest_lambert_azimuthal_equal_area*)nearest; - return grib_nearest_find_generic( - nearest, h, inlat, inlon, flags, /* inputs */ - - self->values_key, /* outputs to set the 'self' object */ - &(self->lats), - &(self->lats_count), - &(self->lons), - &(self->lons_count), - &(self->distances), - - outlats, outlons, /* outputs of the find function */ - values, distances, indexes, len); -} - -static int destroy(grib_nearest* nearest) -{ - grib_nearest_lambert_azimuthal_equal_area* self = (grib_nearest_lambert_azimuthal_equal_area*)nearest; - grib_context* c = grib_context_get_default(); - - if (self->lats) grib_context_free(c, self->lats); - if (self->lons) grib_context_free(c, self->lons); - if (self->i) grib_context_free(c, self->i); - if (self->j) grib_context_free(c, self->j); - if (self->k) grib_context_free(c, self->k); - if (self->distances) grib_context_free(c, self->distances); - return GRIB_SUCCESS; -} diff --git a/src/grib_nearest_class_lambert_conformal.cc b/src/grib_nearest_class_lambert_conformal.cc deleted file mode 100644 index 9f96c6489..000000000 --- a/src/grib_nearest_class_lambert_conformal.cc +++ /dev/null @@ -1,134 +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 = nearest - SUPER = grib_nearest_class_gen - IMPLEMENTS = init - IMPLEMENTS = find;destroy - MEMBERS = double* lats - MEMBERS = int lats_count - MEMBERS = double* lons - MEMBERS = int lons_count - MEMBERS = double* distances - MEMBERS = int* k - MEMBERS = int* i - MEMBERS = int* j - MEMBERS = const char* Ni - MEMBERS = const char* 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_lambert_conformal{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; - /* Members defined in lambert_conformal */ - double* lats; - int lats_count; - double* lons; - int lons_count; - double* distances; - int* k; - int* i; - int* j; - const char* Ni; - const char* Nj; -} grib_nearest_lambert_conformal; - -extern grib_nearest_class* grib_nearest_class_gen; - -static grib_nearest_class _grib_nearest_class_lambert_conformal = { - &grib_nearest_class_gen, /* super */ - "lambert_conformal", /* name */ - sizeof(grib_nearest_lambert_conformal), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_lambert_conformal = &_grib_nearest_class_lambert_conformal; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ - -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) -{ - grib_nearest_lambert_conformal* self = (grib_nearest_lambert_conformal*)nearest; - self->Ni = grib_arguments_get_name(h, args, self->cargs++); - self->Nj = grib_arguments_get_name(h, args, self->cargs++); - self->i = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - self->j = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - return 0; -} - -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, - double* values, double* distances, int* indexes, size_t* len) -{ - grib_nearest_lambert_conformal* self = (grib_nearest_lambert_conformal*)nearest; - return grib_nearest_find_generic( - nearest, h, inlat, inlon, flags, /* inputs */ - - self->values_key, /* outputs to set the 'self' object */ - &(self->lats), - &(self->lats_count), - &(self->lons), - &(self->lons_count), - &(self->distances), - - outlats, outlons, /* outputs of the find function */ - values, distances, indexes, len); -} - -static int destroy(grib_nearest* nearest) -{ - grib_nearest_lambert_conformal* self = (grib_nearest_lambert_conformal*)nearest; - grib_context* c = grib_context_get_default(); - - if (self->lats) grib_context_free(c, self->lats); - if (self->lons) grib_context_free(c, self->lons); - if (self->i) grib_context_free(c, self->i); - if (self->j) grib_context_free(c, self->j); - if (self->k) grib_context_free(c, self->k); - if (self->distances) grib_context_free(c, self->distances); - return GRIB_SUCCESS; -} diff --git a/src/grib_nearest_class_latlon_reduced.cc b/src/grib_nearest_class_latlon_reduced.cc deleted file mode 100644 index ad07f3816..000000000 --- a/src/grib_nearest_class_latlon_reduced.cc +++ /dev/null @@ -1,436 +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 = nearest - SUPER = grib_nearest_class_gen - IMPLEMENTS = init;destroy;find - MEMBERS = double* lats - MEMBERS = int lats_count - MEMBERS = double* lons - MEMBERS = double* distances - MEMBERS = size_t* k - MEMBERS = size_t* j - MEMBERS = const char* Nj - MEMBERS = const char* pl - MEMBERS = const char* lonFirst - MEMBERS = const char* lonLast - 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_latlon_reduced{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; - /* Members defined in latlon_reduced */ - double* lats; - int lats_count; - double* lons; - double* distances; - size_t* k; - size_t* j; - const char* Nj; - const char* pl; - const char* lonFirst; - const char* lonLast; -} grib_nearest_latlon_reduced; - -extern grib_nearest_class* grib_nearest_class_gen; - -static grib_nearest_class _grib_nearest_class_latlon_reduced = { - &grib_nearest_class_gen, /* super */ - "latlon_reduced", /* name */ - sizeof(grib_nearest_latlon_reduced), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_latlon_reduced = &_grib_nearest_class_latlon_reduced; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ - -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) -{ - grib_nearest_latlon_reduced* self = (grib_nearest_latlon_reduced*)nearest; - self->Nj = grib_arguments_get_name(h, args, self->cargs++); - self->pl = grib_arguments_get_name(h, args, self->cargs++); - self->lonFirst = grib_arguments_get_name(h, args, self->cargs++); - self->lonLast = grib_arguments_get_name(h, args, self->cargs++); - self->j = (size_t*)grib_context_malloc(h->context, 2 * sizeof(size_t)); - if (!self->j) - return GRIB_OUT_OF_MEMORY; - self->k = (size_t*)grib_context_malloc(h->context, 4 * sizeof(size_t)); - if (!self->k) - return GRIB_OUT_OF_MEMORY; - - return 0; -} - -static int find_global(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, double* values, - double* distances, int* indexes, size_t* len); - -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, double* values, - double* distances, int* indexes, size_t* len) -{ - int err = 0; - grib_nearest_latlon_reduced* self = (grib_nearest_latlon_reduced*)nearest; - double lat1, lat2, lon1, lon2; - int is_global = 1; - - if (grib_get_double(h, "longitudeFirstInDegrees", &lon1) == GRIB_SUCCESS && - grib_get_double(h, "longitudeLastInDegrees", &lon2) == GRIB_SUCCESS && - grib_get_double(h, "latitudeFirstInDegrees", &lat1) == GRIB_SUCCESS && - grib_get_double(h, "latitudeLastInDegrees", &lat2) == GRIB_SUCCESS) - { - const double difflat = fabs(lat1-lat2); - if (difflat < 180 || lon1 != 0 || lon2 < 359) { - is_global = 0; /* subarea */ - } - } - - if (is_global) { - err = find_global(nearest, h, inlat, inlon, flags, - outlats, outlons, values, - distances, indexes, len); - } - else - { - int lons_count = 0; /*dummy*/ - err = grib_nearest_find_generic( - nearest, h, inlat, inlon, flags, - self->values_key, - &(self->lats), - &(self->lats_count), - &(self->lons), - &(lons_count), - &(self->distances), - outlats, outlons, - values, distances, indexes, len); - } - return err; -} - -static int find_global(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, double* values, - double* distances, int* indexes, size_t* len) -{ - grib_nearest_latlon_reduced* self = (grib_nearest_latlon_reduced*)nearest; - int ret = 0, kk = 0, ii = 0, jj = 0; - int j = 0; - long* pla = NULL; - long* pl = NULL; - size_t nvalues = 0; - grib_iterator* iter = NULL; - double lat = 0, lon = 0; - double radiusInKm; - int ilat = 0, ilon = 0; - - if ((ret = grib_get_size(h, self->values_key, &nvalues)) != GRIB_SUCCESS) - return ret; - nearest->values_count = nvalues; - - if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS) - return ret; - - /* Compute lat/lon info, create iterator etc if it's the 1st time or different grid. - * This is for performance: if the grid has not changed, we only do this once - * and reuse for other messages */ - if (!nearest->h || (flags & GRIB_NEAREST_SAME_GRID) == 0) { - double olat = 1.e10; - long n = 0; - - ilat = 0; - ilon = 0; - if (grib_is_missing(h, self->Nj, &ret)) { - grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->Nj); - return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; - } - - if ((ret = grib_get_long(h, self->Nj, &n)) != GRIB_SUCCESS) - return ret; - self->lats_count = n; - - if (self->lats) - grib_context_free(h->context, self->lats); - self->lats = (double*)grib_context_malloc(h->context, - self->lats_count * sizeof(double)); - if (!self->lats) - return GRIB_OUT_OF_MEMORY; - - if (self->lons) - grib_context_free(h->context, self->lons); - self->lons = (double*)grib_context_malloc(h->context, - nearest->values_count * sizeof(double)); - if (!self->lons) - return GRIB_OUT_OF_MEMORY; - - iter = grib_iterator_new(h, GRIB_GEOITERATOR_NO_VALUES, &ret); - if (ret) { - grib_context_log(h->context, GRIB_LOG_ERROR, "Unable to create iterator"); - return ret; - } - while (grib_iterator_next(iter, &lat, &lon, NULL)) { - if (ilat < self->lats_count && olat != lat) { - self->lats[ilat++] = lat; - olat = lat; - } - self->lons[ilon++] = lon; - } - self->lats_count = ilat; - grib_iterator_delete(iter); - } - nearest->h = h; - - /* Compute distances if it's the 1st time or different point or different grid. - * This is for performance: if the grid and the input point have not changed - * we only do this once and reuse for other messages */ - if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT) == 0 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { - double* lons = NULL; - int nlon = 0; - size_t plsize = 0; - long nplm1 = 0; - int nearest_lons_found = 0; - double lon_first, lon_last; - int islocal = 0; - long plmax; - double dimin; - - if ((ret = grib_get_double(h, self->lonFirst, &lon_first)) != GRIB_SUCCESS) { - grib_context_log(h->context, GRIB_LOG_ERROR, - "grib_nearest_latlon_reduced.find(): unable to get %s %s\n", self->lonFirst, - grib_get_error_message(ret)); - return ret; - } - if ((ret = grib_get_double(h, self->lonLast, &lon_last)) != GRIB_SUCCESS) { - grib_context_log(h->context, GRIB_LOG_ERROR, - "grib_nearest_latlon_reduced.find(): unable to get %s %s\n", self->lonLast, - grib_get_error_message(ret)); - return ret; - } - - plsize = self->lats_count; - if ((ret = grib_get_size(h, self->pl, &plsize)) != GRIB_SUCCESS) - return ret; - pla = (long*)grib_context_malloc(h->context, plsize * sizeof(long)); - if (!pla) - return GRIB_OUT_OF_MEMORY; - if ((ret = grib_get_long_array(h, self->pl, pla, &plsize)) != GRIB_SUCCESS) - return ret; - - pl = pla; - while ((*pl) == 0) { - pl++; - } - - plmax = pla[0]; - for (j = 0; j < plsize; j++) - if (plmax < pla[j]) - plmax = pla[j]; - dimin = 360.0 / plmax; - - if (360 - fabs(lon_last - lon_first) < 2 * dimin) { - islocal = 0; - } - else { - islocal = 1; - } - - if (islocal) - for (j = 0; j < plsize; j++) - pla[j]--; - - /* printf("XXXX islocal=%d\n",islocal); */ - while (inlon < 0) - inlon += 360; - while (inlon > 360) - inlon -= 360; - - ilat = self->lats_count; - if (self->lats[ilat - 1] > self->lats[0]) { - if (inlat < self->lats[0] || inlat > self->lats[ilat - 1]) - return GRIB_OUT_OF_AREA; - } - else { - if (inlat > self->lats[0] || inlat < self->lats[ilat - 1]) - return GRIB_OUT_OF_AREA; - } - - if (!self->distances) - self->distances = (double*)grib_context_malloc(h->context, 4 * sizeof(double)); - if (!self->distances) - return GRIB_OUT_OF_MEMORY; - - grib_binary_search(self->lats, ilat - 1, inlat, - &(self->j[0]), &(self->j[1])); - - nlon = 0; - for (jj = 0; jj < self->j[0]; jj++) - nlon += pl[jj]; - nplm1 = pl[self->j[0]] - 1; - - lons = self->lons + nlon; - - nearest_lons_found = 0; - if (lons[nplm1] > lons[0]) { - if (inlon < lons[0] || inlon > lons[nplm1]) { - if (lons[nplm1] - lons[0] - 360 <= - lons[nplm1] - lons[nplm1 - 1]) { - self->k[0] = 0; - self->k[1] = nplm1; - nearest_lons_found = 1; - } - else - return GRIB_OUT_OF_AREA; - } - } - else { - if (inlon > lons[0] || inlon < lons[nplm1]) { - if (lons[0] - lons[nplm1] - 360 <= - lons[0] - lons[1]) { - self->k[0] = 0; - self->k[1] = nplm1; - nearest_lons_found = 1; - } - else - return GRIB_OUT_OF_AREA; - } - } - - if (!nearest_lons_found) { - grib_binary_search(lons, pl[self->j[0]] - 1, inlon, - &(self->k[0]), &(self->k[1])); - } - self->k[0] += nlon; - self->k[1] += nlon; - - nlon = 0; - for (jj = 0; jj < self->j[1]; jj++) - nlon += pl[jj]; - nplm1 = pl[self->j[1]] - 1; - - lons = self->lons + nlon; - - nearest_lons_found = 0; - if (lons[nplm1] > lons[0]) { - if (inlon < lons[0] || inlon > lons[nplm1]) { - if (lons[nplm1] - lons[0] - 360 <= - lons[nplm1] - lons[nplm1 - 1]) { - self->k[2] = 0; - self->k[3] = nplm1; - nearest_lons_found = 1; - } - else - return GRIB_OUT_OF_AREA; - } - } - else { - if (inlon > lons[0] || inlon < lons[nplm1]) { - if (lons[0] - lons[nplm1] - 360 <= - lons[0] - lons[1]) { - self->k[2] = 0; - self->k[3] = nplm1; - nearest_lons_found = 1; - } - else - return GRIB_OUT_OF_AREA; - } - } - - if (!nearest_lons_found) { - grib_binary_search(lons, pl[self->j[1]] - 1, inlon, - &(self->k[2]), &(self->k[3])); - } - - self->k[2] += nlon; - self->k[3] += nlon; - - kk = 0; - for (jj = 0; jj < 2; jj++) { - for (ii = 0; ii < 2; ii++) { - self->distances[kk] = geographic_distance_spherical(radiusInKm, inlon, inlat, - self->lons[self->k[kk]], self->lats[self->j[jj]]); - kk++; - } - } - - grib_context_free(h->context, pla); - } - - kk = 0; - for (jj = 0; jj < 2; jj++) { - for (ii = 0; ii < 2; ii++) { - distances[kk] = self->distances[kk]; - outlats[kk] = self->lats[self->j[jj]]; - outlons[kk] = self->lons[self->k[kk]]; - if (values) { /* ECC-499 */ - grib_get_double_element_internal(h, self->values_key, self->k[kk], &(values[kk])); - } - indexes[kk] = (int)self->k[kk]; - kk++; - } - } - - return GRIB_SUCCESS; -} - -static int destroy(grib_nearest* nearest) -{ - grib_nearest_latlon_reduced* self = (grib_nearest_latlon_reduced*)nearest; - grib_context* c = grib_context_get_default(); - if (self->lats) - grib_context_free(c, self->lats); - if (self->lons) - grib_context_free(c, self->lons); - if (self->j) - grib_context_free(c, self->j); - if (self->k) - grib_context_free(c, self->k); - if (self->distances) - grib_context_free(c, self->distances); - - return GRIB_SUCCESS; -} diff --git a/src/grib_nearest_class_mercator.cc b/src/grib_nearest_class_mercator.cc deleted file mode 100644 index a160d9dfb..000000000 --- a/src/grib_nearest_class_mercator.cc +++ /dev/null @@ -1,136 +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 = nearest - SUPER = grib_nearest_class_gen - IMPLEMENTS = init - IMPLEMENTS = find;destroy - MEMBERS = double* lats - MEMBERS = int lats_count - MEMBERS = double* lons - MEMBERS = int lons_count - MEMBERS = double* distances - MEMBERS = int* k - MEMBERS = int* i - MEMBERS = int* j - MEMBERS = const char* Ni - MEMBERS = const char* 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_mercator{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; - /* Members defined in mercator */ - double* lats; - int lats_count; - double* lons; - int lons_count; - double* distances; - int* k; - int* i; - int* j; - const char* Ni; - const char* Nj; -} grib_nearest_mercator; - -extern grib_nearest_class* grib_nearest_class_gen; - -static grib_nearest_class _grib_nearest_class_mercator = { - &grib_nearest_class_gen, /* super */ - "mercator", /* name */ - sizeof(grib_nearest_mercator), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_mercator = &_grib_nearest_class_mercator; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ - -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) -{ - grib_nearest_mercator* self = (grib_nearest_mercator*)nearest; - self->Ni = grib_arguments_get_name(h, args, self->cargs++); - self->Nj = grib_arguments_get_name(h, args, self->cargs++); - self->lats = self->lons = self->distances = NULL; - self->lats_count = self->lons_count = 0; - self->i = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - self->j = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - return GRIB_SUCCESS; -} - -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, - double* values, double* distances, int* indexes, size_t* len) -{ - grib_nearest_mercator* self = (grib_nearest_mercator*)nearest; - return grib_nearest_find_generic( - nearest, h, inlat, inlon, flags, /* inputs */ - - self->values_key, /* outputs to set the 'self' object */ - &(self->lats), - &(self->lats_count), - &(self->lons), - &(self->lons_count), - &(self->distances), - - outlats, outlons, /* outputs of the find function */ - values, distances, indexes, len); -} - -static int destroy(grib_nearest* nearest) -{ - grib_nearest_mercator* self = (grib_nearest_mercator*)nearest; - grib_context* c = grib_context_get_default(); - - if (self->lats) grib_context_free(c, self->lats); - if (self->lons) grib_context_free(c, self->lons); - if (self->i) grib_context_free(c, self->i); - if (self->j) grib_context_free(c, self->j); - if (self->k) grib_context_free(c, self->k); - if (self->distances) grib_context_free(c, self->distances); - return GRIB_SUCCESS; -} diff --git a/src/grib_nearest_class_polar_stereographic.cc b/src/grib_nearest_class_polar_stereographic.cc deleted file mode 100644 index 7f04e4ef8..000000000 --- a/src/grib_nearest_class_polar_stereographic.cc +++ /dev/null @@ -1,134 +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 = nearest - SUPER = grib_nearest_class_gen - IMPLEMENTS = init - IMPLEMENTS = find;destroy - MEMBERS = double* lats - MEMBERS = int lats_count - MEMBERS = double* lons - MEMBERS = int lons_count - MEMBERS = double* distances - MEMBERS = int* k - MEMBERS = int* i - MEMBERS = int* j - MEMBERS = const char* Ni - MEMBERS = const char* 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_polar_stereographic{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; - /* Members defined in polar_stereographic */ - double* lats; - int lats_count; - double* lons; - int lons_count; - double* distances; - int* k; - int* i; - int* j; - const char* Ni; - const char* Nj; -} grib_nearest_polar_stereographic; - -extern grib_nearest_class* grib_nearest_class_gen; - -static grib_nearest_class _grib_nearest_class_polar_stereographic = { - &grib_nearest_class_gen, /* super */ - "polar_stereographic", /* name */ - sizeof(grib_nearest_polar_stereographic), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_polar_stereographic = &_grib_nearest_class_polar_stereographic; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ - -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) -{ - grib_nearest_polar_stereographic* self = (grib_nearest_polar_stereographic*)nearest; - self->Ni = grib_arguments_get_name(h, args, self->cargs++); - self->Nj = grib_arguments_get_name(h, args, self->cargs++); - self->i = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - self->j = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - return GRIB_SUCCESS; -} - -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, - double* values, double* distances, int* indexes, size_t* len) -{ - grib_nearest_polar_stereographic* self = (grib_nearest_polar_stereographic*)nearest; - return grib_nearest_find_generic( - nearest, h, inlat, inlon, flags, /* inputs */ - - self->values_key, /* outputs to set the 'self' object */ - &(self->lats), - &(self->lats_count), - &(self->lons), - &(self->lons_count), - &(self->distances), - - outlats, outlons, /* outputs of the find function */ - values, distances, indexes, len); -} - -static int destroy(grib_nearest* nearest) -{ - grib_nearest_polar_stereographic* self = (grib_nearest_polar_stereographic*)nearest; - grib_context* c = grib_context_get_default(); - - if (self->lats) grib_context_free(c, self->lats); - if (self->lons) grib_context_free(c, self->lons); - if (self->i) grib_context_free(c, self->i); - if (self->j) grib_context_free(c, self->j); - if (self->k) grib_context_free(c, self->k); - if (self->distances) grib_context_free(c, self->distances); - return GRIB_SUCCESS; -} diff --git a/src/grib_nearest_class_regular.cc b/src/grib_nearest_class_regular.cc deleted file mode 100644 index 5d591a6e3..000000000 --- a/src/grib_nearest_class_regular.cc +++ /dev/null @@ -1,373 +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 = nearest - SUPER = grib_nearest_class_gen - IMPLEMENTS = init - IMPLEMENTS = find;destroy - MEMBERS = double* lats - MEMBERS = int lats_count - MEMBERS = double* lons - MEMBERS = int lons_count - MEMBERS = double* distances - MEMBERS = size_t* k - MEMBERS = size_t* i - MEMBERS = size_t* j - MEMBERS = const char* Ni - MEMBERS = const char* 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_regular{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; - /* Members defined in regular */ - double* lats; - int lats_count; - double* lons; - int lons_count; - double* distances; - size_t* k; - size_t* i; - size_t* j; - const char* Ni; - const char* Nj; -} grib_nearest_regular; - -extern grib_nearest_class* grib_nearest_class_gen; - -static grib_nearest_class _grib_nearest_class_regular = { - &grib_nearest_class_gen, /* super */ - "regular", /* name */ - sizeof(grib_nearest_regular), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_regular = &_grib_nearest_class_regular; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ - - -#define NUM_NEIGHBOURS 4 - -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) -{ - grib_nearest_regular* self = (grib_nearest_regular*)nearest; - self->Ni = grib_arguments_get_name(h, args, self->cargs++); - self->Nj = grib_arguments_get_name(h, args, self->cargs++); - self->i = (size_t*)grib_context_malloc(h->context, 2 * sizeof(size_t)); - self->j = (size_t*)grib_context_malloc(h->context, 2 * sizeof(size_t)); - return 0; -} - -static bool is_rotated_grid(grib_handle* h) -{ - long is_rotated = 0; - int err = grib_get_long(h, "isRotatedGrid", &is_rotated); - if (!err && is_rotated) - return true; - return false; -} - -// Old implementation in -// src/deprecated/grib_nearest_class_regular.cc -// -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, - double* values, double* distances, int* indexes, size_t* len) -{ - grib_nearest_regular* self = (grib_nearest_regular*)nearest; - int ret = 0, kk = 0, ii = 0, jj = 0; - size_t nvalues = 0; - double radiusInKm; - - grib_iterator* iter = NULL; - double lat = 0, lon = 0; - const bool is_rotated = is_rotated_grid(h); - double angleOfRotation = 0, southPoleLat = 0, southPoleLon = 0; - grib_context* c = h->context; - - while (inlon < 0) - inlon += 360; - while (inlon > 360) - inlon -= 360; - - if ((ret = grib_get_size(h, self->values_key, &nvalues)) != GRIB_SUCCESS) - return ret; - nearest->values_count = nvalues; - - if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS) - return ret; - - /* Compute lat/lon info, create iterator etc if it's the 1st time or different grid. - * This is for performance: if the grid has not changed, we only do this once - * and reuse for other messages */ - if (!nearest->h || (flags & GRIB_NEAREST_SAME_GRID) == 0) { - double olat = 1.e10, olon = 1.e10; - int ilat = 0, ilon = 0; - long n = 0; - - if (grib_is_missing(h, self->Ni, &ret)) { - grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->Ni); - return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; - } - - if (grib_is_missing(h, self->Nj, &ret)) { - grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->Nj); - return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; - } - - /* ECC-600: Support for rotated grids - * First: rotate the input point - * Then: run the lat/lon iterator over the rotated grid (disableUnrotate) - * Finally: unrotate the resulting point - */ - if (is_rotated) { - double new_lat = 0, new_lon = 0; - ret = grib_get_double_internal(h, "angleOfRotation", &angleOfRotation); - if (ret) - return ret; - ret = grib_get_double_internal(h, "latitudeOfSouthernPoleInDegrees", &southPoleLat); - if (ret) - return ret; - ret = grib_get_double_internal(h, "longitudeOfSouthernPoleInDegrees", &southPoleLon); - if (ret) - return ret; - ret = grib_set_long(h, "iteratorDisableUnrotate", 1); - if (ret) - return ret; - /* Rotate the inlat, inlon */ - rotate(inlat, inlon, angleOfRotation, southPoleLat, southPoleLon, &new_lat, &new_lon); - inlat = new_lat; - inlon = new_lon; - /*if(h->context->debug) printf("nearest find: rotated grid: new point=(%g,%g)\n",new_lat,new_lon);*/ - } - - if ((ret = grib_get_long(h, self->Ni, &n)) != GRIB_SUCCESS) - return ret; - self->lons_count = n; - - if ((ret = grib_get_long(h, self->Nj, &n)) != GRIB_SUCCESS) - return ret; - self->lats_count = n; - - if (self->lats) - grib_context_free(c, self->lats); - self->lats = (double*)grib_context_malloc(c, self->lats_count * sizeof(double)); - if (!self->lats) - return GRIB_OUT_OF_MEMORY; - - if (self->lons) - grib_context_free(c, self->lons); - self->lons = (double*)grib_context_malloc(c, self->lons_count * sizeof(double)); - if (!self->lons) - return GRIB_OUT_OF_MEMORY; - - iter = grib_iterator_new(h, GRIB_GEOITERATOR_NO_VALUES, &ret); - if (ret != GRIB_SUCCESS) { - grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_regular: Unable to create lat/lon iterator"); - return ret; - } - while (grib_iterator_next(iter, &lat, &lon, NULL)) { - if (ilat < self->lats_count && olat != lat) { - self->lats[ilat++] = lat; - olat = lat; - } - if (ilon < self->lons_count && olon != lon) { - self->lons[ilon++] = lon; - olon = lon; - } - } - grib_iterator_delete(iter); - } - nearest->h = h; - - /* Compute distances if it's the 1st time or different point or different grid. - * This is for performance: if the grid and the input point have not changed - * we only do this once and reuse for other messages */ - if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT) == 0 || (flags & GRIB_NEAREST_SAME_GRID) == 0) { - int nearest_lons_found = 0; - - if (self->lats[self->lats_count - 1] > self->lats[0]) { - if (inlat < self->lats[0] || inlat > self->lats[self->lats_count - 1]) - return GRIB_OUT_OF_AREA; - } - else { - if (inlat > self->lats[0] || inlat < self->lats[self->lats_count - 1]) - return GRIB_OUT_OF_AREA; - } - - if (self->lons[self->lons_count - 1] > self->lons[0]) { - if (inlon < self->lons[0] || inlon > self->lons[self->lons_count - 1]) { - /* try to scale*/ - if (inlon > 0) - inlon -= 360; - else - inlon += 360; - - if (inlon < self->lons[0] || inlon > self->lons[self->lons_count - 1]) { - if (self->lons[0] + 360 - self->lons[self->lons_count - 1] <= - self->lons[1] - self->lons[0]) { - /*it's a global field in longitude*/ - self->i[0] = 0; - self->i[1] = self->lons_count - 1; - nearest_lons_found = 1; - } - else - return GRIB_OUT_OF_AREA; - } - } - } - else { - if (inlon > self->lons[0] || inlon < self->lons[self->lons_count - 1]) { - /* try to scale*/ - if (inlon > 0) - inlon -= 360; - else - inlon += 360; - if (self->lons[0] - self->lons[self->lons_count - 1] - 360 <= - self->lons[0] - self->lons[1]) { - /*it's a global field in longitude*/ - self->i[0] = 0; - self->i[1] = self->lons_count - 1; - nearest_lons_found = 1; - } - else if (inlon > self->lons[0] || inlon < self->lons[self->lons_count - 1]) - return GRIB_OUT_OF_AREA; - } - } - - grib_binary_search(self->lats, self->lats_count - 1, inlat, - &(self->j[0]), &(self->j[1])); - - if (!nearest_lons_found) - grib_binary_search(self->lons, self->lons_count - 1, inlon, - &(self->i[0]), &(self->i[1])); - - if (!self->distances) - self->distances = (double*)grib_context_malloc(c, NUM_NEIGHBOURS * sizeof(double)); - if (!self->k) - self->k = (size_t*)grib_context_malloc(c, NUM_NEIGHBOURS * sizeof(size_t)); - kk = 0; - for (jj = 0; jj < 2; jj++) { - for (ii = 0; ii < 2; ii++) { - self->k[kk] = self->i[ii] + self->lons_count * self->j[jj]; - self->distances[kk] = geographic_distance_spherical(radiusInKm, inlon, inlat, - self->lons[self->i[ii]], self->lats[self->j[jj]]); - kk++; - } - } - } - - kk = 0; - - /* - * Brute force algorithm: - * First unpack all the values into an array. Then when we need the 4 points - * we just index into this array so no need to call grib_get_double_element_internal - * - * if (nearest->values) grib_context_free(c,nearest->values); - * nearest->values = grib_context_malloc(h->context,nvalues*sizeof(double)); - * if (!nearest->values) return GRIB_OUT_OF_MEMORY; - * ret = grib_get_double_array(h, self->values_key, nearest->values ,&nvalues); - * if (ret) return ret; - */ - - if (values) { - /* See ECC-1403 and ECC-499 */ - /* Performance: Decode the field once and get all 4 values */ - if ((ret = grib_get_double_element_set(h, self->values_key, self->k, NUM_NEIGHBOURS, values)) != GRIB_SUCCESS) - return ret; - } - - for (jj = 0; jj < 2; jj++) { - for (ii = 0; ii < 2; ii++) { - distances[kk] = self->distances[kk]; - outlats[kk] = self->lats[self->j[jj]]; - outlons[kk] = self->lons[self->i[ii]]; - if (is_rotated) { - /* Unrotate resulting lat/lon */ - double new_lat = 0, new_lon = 0; - unrotate(outlats[kk], outlons[kk], angleOfRotation, southPoleLat, southPoleLon, &new_lat, &new_lon); - outlats[kk] = new_lat; - outlons[kk] = new_lon; - } - /* See ECC-1403 and ECC-499 - * if (values) { - * grib_get_double_element_internal(h, self->values_key, self->k[kk], &(values[kk])); - *} - */ - /* Using the brute force approach described above */ - /* Assert(self->k[kk] < nvalues); */ - /* values[kk]=nearest->values[self->k[kk]]; */ - - if (self->k[kk] >= INT_MAX) { - /* Current interface uses an 'int' for 'indexes' which is 32bits! We should change this to a 64bit type */ - grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_regular: Unable to compute index. Value too large"); - return GRIB_OUT_OF_RANGE; - } else { - indexes[kk] = (int)self->k[kk]; - } - kk++; - } - } - - return GRIB_SUCCESS; -} - -static int destroy(grib_nearest* nearest) -{ - grib_nearest_regular* self = (grib_nearest_regular*)nearest; - grib_context* c = grib_context_get_default(); - - if (self->lats) grib_context_free(c, self->lats); - if (self->lons) grib_context_free(c, self->lons); - if (self->i) grib_context_free(c, self->i); - if (self->j) grib_context_free(c, self->j); - if (self->k) grib_context_free(c, self->k); - if (self->distances) grib_context_free(c, self->distances); - return GRIB_SUCCESS; -} diff --git a/src/grib_nearest_class_space_view.cc b/src/grib_nearest_class_space_view.cc deleted file mode 100644 index ccba897b9..000000000 --- a/src/grib_nearest_class_space_view.cc +++ /dev/null @@ -1,134 +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 = nearest - SUPER = grib_nearest_class_gen - IMPLEMENTS = init - IMPLEMENTS = find;destroy - MEMBERS = double* lats - MEMBERS = int lats_count - MEMBERS = double* lons - MEMBERS = int lons_count - MEMBERS = double* distances - MEMBERS = int* k - MEMBERS = int* i - MEMBERS = int* j - MEMBERS = const char* Ni - MEMBERS = const char* 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 "nearest.class" and rerun ./make_class.pl - -*/ - - -static void init_class (grib_nearest_class*); - -static int init (grib_nearest* nearest,grib_handle* h,grib_arguments* args); -static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len); -static int destroy (grib_nearest* nearest); - -typedef struct grib_nearest_space_view{ - grib_nearest nearest; - /* Members defined in gen */ - const char* values_key; - const char* radius; - int cargs; - /* Members defined in space_view */ - double* lats; - int lats_count; - double* lons; - int lons_count; - double* distances; - int* k; - int* i; - int* j; - const char* Ni; - const char* Nj; -} grib_nearest_space_view; - -extern grib_nearest_class* grib_nearest_class_gen; - -static grib_nearest_class _grib_nearest_class_space_view = { - &grib_nearest_class_gen, /* super */ - "space_view", /* name */ - sizeof(grib_nearest_space_view), /* size of instance */ - 0, /* inited */ - &init_class, /* init_class */ - &init, /* constructor */ - &destroy, /* destructor */ - &find, /* find nearest */ -}; - -grib_nearest_class* grib_nearest_class_space_view = &_grib_nearest_class_space_view; - - -static void init_class(grib_nearest_class* c) -{ -} -/* END_CLASS_IMP */ - -static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args) -{ - grib_nearest_space_view* self = (grib_nearest_space_view*)nearest; - self->Ni = grib_arguments_get_name(h, args, self->cargs++); - self->Nj = grib_arguments_get_name(h, args, self->cargs++); - self->i = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - self->j = (int*)grib_context_malloc(h->context, 2 * sizeof(int)); - return 0; -} - -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon, unsigned long flags, - double* outlats, double* outlons, - double* values, double* distances, int* indexes, size_t* len) -{ - grib_nearest_space_view* self = (grib_nearest_space_view*)nearest; - return grib_nearest_find_generic( - nearest, h, inlat, inlon, flags, /* inputs */ - - self->values_key, /* outputs to set the 'self' object */ - &(self->lats), - &(self->lats_count), - &(self->lons), - &(self->lons_count), - &(self->distances), - - outlats, outlons, /* outputs of the find function */ - values, distances, indexes, len); -} - -static int destroy(grib_nearest* nearest) -{ - grib_nearest_space_view* self = (grib_nearest_space_view*)nearest; - grib_context* c = grib_context_get_default(); - - if (self->lats) grib_context_free(c, self->lats); - if (self->lons) grib_context_free(c, self->lons); - if (self->i) grib_context_free(c, self->i); - if (self->j) grib_context_free(c, self->j); - if (self->k) grib_context_free(c, self->k); - if (self->distances) grib_context_free(c, self->distances); - return GRIB_SUCCESS; -} diff --git a/src/grib_nearest_factory.cc b/src/grib_nearest_factory.cc new file mode 100644 index 000000000..bfe6c9ecc --- /dev/null +++ b/src/grib_nearest_factory.cc @@ -0,0 +1,68 @@ +/* + * (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_nearest_factory.h" +#include "accessor/grib_accessor_class_nearest.h" + +extern eccodes::geo_nearest::Nearest* grib_nearest_healpix; +extern eccodes::geo_nearest::Nearest* grib_nearest_lambert_azimuthal_equal_area; +extern eccodes::geo_nearest::Nearest* grib_nearest_lambert_conformal; +extern eccodes::geo_nearest::Nearest* grib_nearest_latlon_reduced; +extern eccodes::geo_nearest::Nearest* grib_nearest_mercator; +extern eccodes::geo_nearest::Nearest* grib_nearest_polar_stereographic; +extern eccodes::geo_nearest::Nearest* grib_nearest_reduced; +extern eccodes::geo_nearest::Nearest* grib_nearest_regular; +extern eccodes::geo_nearest::Nearest* grib_nearest_space_view; + +struct table_entry +{ + const char* type; + eccodes::geo_nearest::Nearest** nearest; +}; + +static const struct table_entry table[] = { + { "healpix", &grib_nearest_healpix, }, + { "lambert_azimuthal_equal_area", &grib_nearest_lambert_azimuthal_equal_area, }, + { "lambert_conformal", &grib_nearest_lambert_conformal, }, + { "latlon_reduced", &grib_nearest_latlon_reduced, }, + { "mercator", &grib_nearest_mercator, }, + { "polar_stereographic", &grib_nearest_polar_stereographic, }, + { "reduced", &grib_nearest_reduced, }, + { "regular", &grib_nearest_regular, }, + { "space_view", &grib_nearest_space_view, }, +}; + +eccodes::geo_nearest::Nearest* grib_nearest_factory(grib_handle* h, grib_arguments* args, int* error) +{ + size_t i = 0, num_table_entries = 0; + *error = GRIB_NOT_IMPLEMENTED; + char* type = (char*)grib_arguments_get_name(h, args, 0); + + num_table_entries = sizeof(table) / sizeof(table[0]); + for (i = 0; i < num_table_entries; i++) { + if (strcmp(type, table[i].type) == 0) { + eccodes::geo_nearest::Nearest* creator = *(table[i].nearest); + eccodes::geo_nearest::Nearest* it = creator->create(); + + *error = it->init(h, args); + + if (*error == GRIB_SUCCESS) + return it; + grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_factory: Error instantiating nearest %s (%s)", + table[i].type, grib_get_error_message(*error)); + gribNearestDelete(it); + return NULL; + } + } + + grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_factory: Unknown type: %s", type); + + return NULL; +} diff --git a/src/grib_nearest_factory.h b/src/grib_nearest_factory.h index 837524914..b203f7590 100644 --- a/src/grib_nearest_factory.h +++ b/src/grib_nearest_factory.h @@ -1,11 +1,16 @@ -/* This file is automatically generated by ./make_class.pl, do not edit */ -{ "gen", &grib_nearest_class_gen, }, -{ "healpix", &grib_nearest_class_healpix, }, -{ "lambert_azimuthal_equal_area", &grib_nearest_class_lambert_azimuthal_equal_area, }, -{ "lambert_conformal", &grib_nearest_class_lambert_conformal, }, -{ "latlon_reduced", &grib_nearest_class_latlon_reduced, }, -{ "mercator", &grib_nearest_class_mercator, }, -{ "polar_stereographic", &grib_nearest_class_polar_stereographic, }, -{ "reduced", &grib_nearest_class_reduced, }, -{ "regular", &grib_nearest_class_regular, }, -{ "space_view", &grib_nearest_class_space_view, }, +/* + * (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" +#include "geo_nearest/grib_nearest.h" + +eccodes::geo_nearest::Nearest* grib_nearest_factory(grib_handle* h, grib_arguments* args, int* error);