mirror of https://github.com/ecmwf/eccodes.git
520 lines
18 KiB
C++
520 lines
18 KiB
C++
/*
|
|
* (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 = long global
|
|
MEMBERS = double lon_first
|
|
MEMBERS = double lon_last
|
|
MEMBERS = int legacy
|
|
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;
|
|
} 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 */
|
|
|
|
#define NUM_NEIGHBOURS 4
|
|
|
|
static int init(grib_nearest* nearest, 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;
|
|
if (!self->j)
|
|
return GRIB_OUT_OF_MEMORY;
|
|
self->k = (size_t*)grib_context_malloc(nearest->context, NUM_NEIGHBOURS * sizeof(size_t));
|
|
if (!self->k)
|
|
return GRIB_OUT_OF_MEMORY;
|
|
grib_get_long(h, "global", &self->global);
|
|
if (!self->global) {
|
|
int err;
|
|
/*TODO longitudeOfFirstGridPointInDegrees from the def file*/
|
|
if ((err = grib_get_double(h, "longitudeOfFirstGridPointInDegrees", &self->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) {
|
|
grib_context_log(h->context, GRIB_LOG_ERROR,
|
|
"grib_nearest_reduced: Unable to get longitudeOfLastGridPointInDegrees %s\n",
|
|
grib_get_error_message(err));
|
|
return err;
|
|
}
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
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;
|
|
long lLegacy = 0;
|
|
err = grib_get_long(h, "legacyGaussSubarea", &lLegacy);
|
|
if (err) return err;
|
|
*legacy = (int)lLegacy;
|
|
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)
|
|
{
|
|
int err = 0;
|
|
grib_nearest_reduced* self = (grib_nearest_reduced*)nearest;
|
|
|
|
if (self->global) {
|
|
err = find_global(nearest, h,
|
|
inlat, inlon, flags,
|
|
outlats, outlons, values,
|
|
distances, indexes, len);
|
|
}
|
|
else
|
|
{
|
|
/* ECC-762, ECC-1432: Use brute force generic algorithm
|
|
* for reduced grid subareas. Review in the future
|
|
*/
|
|
int lons_count = 0; /*dummy*/
|
|
|
|
err = grib_nearest_find_generic(
|
|
nearest, h, inlat, inlon, flags,
|
|
self->values_key,
|
|
"Ni",
|
|
self->Nj,
|
|
&(self->lats),
|
|
&(self->lats_count),
|
|
&(self->lons),
|
|
&(lons_count),
|
|
&(self->distances),
|
|
outlats, outlons,
|
|
values, distances, indexes, len);
|
|
}
|
|
return err;
|
|
}
|
|
|
|
/* Old implementation in src/deprecated/grib_nearest_class_reduced.old */
|
|
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_reduced* self = (grib_nearest_reduced*)nearest;
|
|
int err = 0, kk = 0, ii = 0;
|
|
size_t jj = 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;
|
|
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 (err) return err;
|
|
}
|
|
if (self->legacy == 1) {
|
|
get_reduced_row_func = &grib_get_reduced_row_legacy;
|
|
}
|
|
|
|
if ((err = grib_get_size(h, self->values_key, &nvalues)) != GRIB_SUCCESS)
|
|
return err;
|
|
nearest->values_count = nvalues;
|
|
|
|
if ((err = grib_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS)
|
|
return err;
|
|
|
|
/* 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, &err)) {
|
|
grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->Nj);
|
|
return err ? err : GRIB_GEOCALCULUS_PROBLEM;
|
|
}
|
|
|
|
if ((err = grib_get_long(h, self->Nj, &n)) != GRIB_SUCCESS)
|
|
return err;
|
|
self->lats_count = n;
|
|
|
|
if (self->lats)
|
|
grib_context_free(nearest->context, self->lats);
|
|
self->lats = (double*)grib_context_malloc(nearest->context,
|
|
self->lats_count * sizeof(double));
|
|
if (!self->lats)
|
|
return GRIB_OUT_OF_MEMORY;
|
|
|
|
if (self->lons)
|
|
grib_context_free(nearest->context, self->lons);
|
|
self->lons = (double*)grib_context_malloc(nearest->context,
|
|
nearest->values_count * sizeof(double));
|
|
if (!self->lons)
|
|
return GRIB_OUT_OF_MEMORY;
|
|
|
|
iter = grib_iterator_new(h, GRIB_GEOITERATOR_NO_VALUES, &err);
|
|
if (err != GRIB_SUCCESS) {
|
|
grib_context_log(h->context, GRIB_LOG_ERROR, "grib_nearest_reduced: Unable to create lat/lon iterator");
|
|
return err;
|
|
}
|
|
while (grib_iterator_next(iter, &lat, &lon, NULL)) {
|
|
if (olat != lat) {
|
|
self->lats[ilat++] = lat;
|
|
olat = lat;
|
|
}
|
|
while (lon > 360)
|
|
lon -= 360;
|
|
if (!self->global) { /* ECC-756 */
|
|
if (self->legacy == 0) /*TODO*/
|
|
if (lon > 180 && lon < 360)
|
|
lon -= 360;
|
|
}
|
|
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;
|
|
long row_count, ilon_first, ilon_last;
|
|
|
|
if (self->global) {
|
|
inlon = normalise_longitude_in_degrees(inlon);
|
|
}
|
|
else {
|
|
/* TODO: Experimental */
|
|
if (self->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])
|
|
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(nearest->context, NUM_NEIGHBOURS * sizeof(double));
|
|
if (!self->distances)
|
|
return GRIB_OUT_OF_MEMORY;
|
|
|
|
grib_binary_search(self->lats, ilat - 1, inlat,
|
|
&(self->j[0]), &(self->j[1]));
|
|
|
|
plsize = self->lats_count;
|
|
if ((err = grib_get_size(h, self->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)
|
|
return err;
|
|
|
|
pl = pla;
|
|
while ((*pl) == 0) {
|
|
pl++;
|
|
}
|
|
|
|
nlon = 0;
|
|
if (self->global) {
|
|
for (jj = 0; jj < self->j[0]; jj++)
|
|
nlon += pl[jj];
|
|
nplm1 = pl[self->j[0]] - 1;
|
|
}
|
|
else {
|
|
nlon = 0;
|
|
for (jj = 0; jj < self->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);
|
|
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);
|
|
nplm1 = row_count - 1;
|
|
}
|
|
lons = self->lons + nlon;
|
|
|
|
nearest_lons_found = 0;
|
|
/* ECC-756: The comparisons of longitudes here depends on the longitude values
|
|
* from the point iterator. The old values could be -ve but the new algorithm
|
|
* generates +ve values which break this test:
|
|
* lons[nplm1]>lons[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) {
|
|
if (!self->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);
|
|
}
|
|
else {
|
|
row_count = pl[self->j[0]];
|
|
}
|
|
|
|
grib_binary_search(lons, row_count - 1, inlon,
|
|
&(self->k[0]), &(self->k[1]));
|
|
}
|
|
self->k[0] += nlon;
|
|
self->k[1] += nlon;
|
|
|
|
nlon = 0;
|
|
if (self->global) {
|
|
for (jj = 0; jj < self->j[1]; jj++)
|
|
nlon += pl[jj];
|
|
nplm1 = pl[self->j[1]] - 1;
|
|
}
|
|
else {
|
|
for (jj = 0; jj < self->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);
|
|
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);
|
|
nplm1--;
|
|
}
|
|
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) {
|
|
if (!self->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);
|
|
}
|
|
else {
|
|
row_count = pl[self->j[1]];
|
|
}
|
|
|
|
grib_binary_search(lons, row_count - 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;
|
|
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);
|
|
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]];
|
|
/*if (values) {
|
|
* grib_get_double_element_internal(h, self->values_key, self->k[kk], &(values[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_reduced: Unable to compute index. Value too large");
|
|
return GRIB_OUT_OF_RANGE;
|
|
} else {
|
|
indexes[kk] = self->k[kk];
|
|
}
|
|
kk++;
|
|
}
|
|
}
|
|
|
|
return GRIB_SUCCESS;
|
|
}
|
|
|
|
static int destroy(grib_nearest* nearest)
|
|
{
|
|
grib_nearest_reduced* self = (grib_nearest_reduced*)nearest;
|
|
if (self->lats)
|
|
grib_context_free(nearest->context, self->lats);
|
|
if (self->lons)
|
|
grib_context_free(nearest->context, self->lons);
|
|
if (self->j)
|
|
grib_context_free(nearest->context, self->j);
|
|
if (self->k)
|
|
grib_context_free(nearest->context, self->k);
|
|
if (self->distances)
|
|
grib_context_free(nearest->context, self->distances);
|
|
|
|
return GRIB_SUCCESS;
|
|
}
|