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

Modernisation: Nearest hierarchy

This commit is contained in:
Eugen Betke 2024-11-03 21:01:15 +01:00
parent 1c45f774de
commit 2a330d3b2a
43 changed files with 2137 additions and 2404 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/CMakeLists.txt
View File

@ -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

31
src/accessor/grib_accessor_class_nearest.cc
View File

@ -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

6
src/accessor/grib_accessor_class_nearest.h
View File

@ -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);

30
src/eccodes_prototypes.h
View File

@ -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);

553
src/grib_nearest.cc → src/geo_nearest/grib_nearest.cc
View File

@ -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

48
src/geo_nearest/grib_nearest.h Normal file
View File

@ -0,0 +1,48 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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);

47
src/geo_nearest/grib_nearest_class_gen.cc Normal file
View File

@ -0,0 +1,47 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#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

46
src/geo_nearest/grib_nearest_class_gen.h Normal file
View File

@ -0,0 +1,46 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

65
src/geo_nearest/grib_nearest_class_healpix.cc Normal file
View File

@ -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

40
src/geo_nearest/grib_nearest_class_healpix.h Normal file
View File

@ -0,0 +1,40 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

64
src/geo_nearest/grib_nearest_class_lambert_azimuthal_equal_area.cc Normal file
View File

@ -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

40
src/geo_nearest/grib_nearest_class_lambert_azimuthal_equal_area.h Normal file
View File

@ -0,0 +1,40 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

65
src/geo_nearest/grib_nearest_class_lambert_conformal.cc Normal file
View File

@ -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

40
src/geo_nearest/grib_nearest_class_lambert_conformal.h Normal file
View File

@ -0,0 +1,40 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

360
src/geo_nearest/grib_nearest_class_latlon_reduced.cc Normal file
View File

@ -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

46
src/geo_nearest/grib_nearest_class_latlon_reduced.h Normal file
View File

@ -0,0 +1,46 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

65
src/geo_nearest/grib_nearest_class_mercator.cc Normal file
View File

@ -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

40
src/geo_nearest/grib_nearest_class_mercator.h Normal file
View File

@ -0,0 +1,40 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

63
src/geo_nearest/grib_nearest_class_polar_stereographic.cc Normal file
View File

@ -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

40
src/geo_nearest/grib_nearest_class_polar_stereographic.h Normal file
View File

@ -0,0 +1,40 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

332
src/grib_nearest_class_reduced.cc → src/geo_nearest/grib_nearest_class_reduced.cc
View File

@ -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

48
src/geo_nearest/grib_nearest_class_reduced.h Normal file
View File

@ -0,0 +1,48 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

302
src/geo_nearest/grib_nearest_class_regular.cc Normal file
View File

@ -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

40
src/geo_nearest/grib_nearest_class_regular.h Normal file
View File

@ -0,0 +1,40 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

65
src/geo_nearest/grib_nearest_class_space_view.cc Normal file
View File

@ -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

40
src/geo_nearest/grib_nearest_class_space_view.h Normal file
View File

@ -0,0 +1,40 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_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

41
src/grib_api_internal.h
View File

@ -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

13
src/grib_iterator_class.h
View File

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

36
src/grib_iterator_class.cc → src/grib_iterator_factory.cc
View File

@ -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)
{

16
src/grib_iterator_factory.h
View File

@ -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);

52
src/grib_nearest_class.cc
View File

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

11
src/grib_nearest_class.h
View File

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

102
src/grib_nearest_class_gen.cc
View File

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

136
src/grib_nearest_class_healpix.cc
View File

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

134
src/grib_nearest_class_lambert_azimuthal_equal_area.cc
View File

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

134
src/grib_nearest_class_lambert_conformal.cc
View File

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

436
src/grib_nearest_class_latlon_reduced.cc
View File

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

136
src/grib_nearest_class_mercator.cc
View File

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

134
src/grib_nearest_class_polar_stereographic.cc
View File

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

373
src/grib_nearest_class_regular.cc
View File

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

134
src/grib_nearest_class_space_view.cc
View File

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

68
src/grib_nearest_factory.cc Normal file
View File

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

27
src/grib_nearest_factory.h
View File

@ -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);