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

Merge branch 'develop' into feature/LAEA_Oblate

This commit is contained in:
Shahram Najm 2020-06-24 21:47:20 +01:00
parent 77b632fedc d291c25cd6
commit fbdee2773a
41 changed files with 906 additions and 488 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.appveyor.yml
View File

@ -2,7 +2,7 @@
# general configuration #
#---------------------------------#
version: 2.18.0-{build}-{branch}
version: 2.19.0-{build}-{branch}
branches:
only:

4
CMakeLists.txt
View File

@ -18,7 +18,7 @@
cmake_minimum_required( VERSION 3.6 FATAL_ERROR )
project( eccodes VERSION 2.18.0 LANGUAGES C )
project( eccodes VERSION 2.19.0 LANGUAGES C )
set( CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/../ecbuild/cmake" )
@ -119,7 +119,7 @@ ecbuild_add_option( FEATURE AEC
ecbuild_add_option( FEATURE PYTHON
DESCRIPTION "Build the ecCodes Python2 interface (deprecated)"
DEFAULT ON
DEFAULT OFF
REQUIRED_PACKAGES "Python VERSION 2.6 NO_LIBS" NumPy )
# For Python2 we build our own bindings (using SWIG) in the build directory
# but for Python3 one has to add the eccodes from pip3 AFTER the install

2
definitions/grib1/grid_definition_1.def
View File

@ -64,6 +64,8 @@ iterator mercator(numberOfPoints,missingValue,values,
jPointsAreConsecutive,
alternativeRowScanning);
nearest mercator(values,radius,Nx,Ny);
meta latLonValues latlonvalues(values);
alias latitudeLongitudeValues=latLonValues;
meta latitudes latitudes(values,0);

2
definitions/grib2/template.3.10.def
View File

@ -75,6 +75,8 @@ iterator mercator(numberOfPoints,missingValue,values,
iScansNegatively, jScansPositively,
jPointsAreConsecutive, alternativeRowScanning);
nearest mercator(values,radius,Nx,Ny);
meta latLonValues latlonvalues(values);
alias latitudeLongitudeValues=latLonValues;
meta latitudes latitudes(values,0);

2
definitions/grib2/template.3.140.def
View File

@ -60,6 +60,8 @@ iterator lambert_azimuthal_equal_area(numberOfPoints,missingValue,values,
iScansNegatively, jScansPositively,
jPointsAreConsecutive, alternativeRowScanning);
nearest lambert_azimuthal_equal_area(values,radius,Nx,Ny);
meta latLonValues latlonvalues(values);
alias latitudeLongitudeValues=latLonValues;
meta latitudes latitudes(values,0);

2
definitions/grib2/template.3.90.def
View File

@ -80,6 +80,8 @@ iterator space_view(numberOfPoints, missingValue, values, radius,
iScansNegatively, jScansPositively,
jPointsAreConsecutive, alternativeRowScanning);
nearest space_view(values,radius,Nx,Ny);
meta latLonValues latlonvalues(values);
alias latitudeLongitudeValues=latLonValues;
meta latitudes latitudes(values,0);

2
definitions/mars/grib.elda.eme.def Normal file
View File

@ -0,0 +1,2 @@
alias mars.anoffset=offsetToEndOf4DvarWindow;
alias mars.number=componentIndex;

1
definitions/mars/grib.enda.eme.def Normal file
View File

@ -0,0 +1 @@
alias mars.number=componentIndex;

3
src/CMakeLists.txt
View File

@ -346,7 +346,10 @@ list( APPEND grib_api_srcs
grib_nearest_class_latlon_reduced.c
grib_nearest_class_sh.c
grib_nearest_class_lambert_conformal.c
grib_nearest_class_lambert_azimuthal_equal_area.c
grib_nearest_class_mercator.c
grib_nearest_class_polar_stereographic.c
grib_nearest_class_space_view.c
grib_iterator_class_polar_stereographic.c
grib_iterator_class_lambert_azimuthal_equal_area.c
grib_iterator_class_lambert_conformal.c

3
src/Makefile.am
View File

@ -361,7 +361,10 @@ libeccodes_la_prototypes= \
grib_nearest_class_latlon_reduced.c \
grib_nearest_class_sh.c \
grib_nearest_class_lambert_conformal.c \
grib_nearest_class_lambert_azimuthal_equal_area.c \
grib_nearest_class_mercator.c \
grib_nearest_class_polar_stereographic.c \
grib_nearest_class_space_view.c \
grib_iterator_class_polar_stereographic.c \
grib_iterator_class_lambert_azimuthal_equal_area.c \
grib_iterator_class_lambert_conformal.c \

12
src/grib_api_prototypes.h
View File

@ -1048,6 +1048,8 @@ int grib_get_gaussian_latitudes(long trunc, double* lats);
int is_gaussian_global(double lat1, double lat2, double lon1, double lon2, long num_points_equator, const double* latitudes, double angular_precision);
void rotate(const double inlat, const double inlon, const double angleOfRot, const double southPoleLat, const double southPoleLon, double* outlat, double* outlon);
void unrotate(const double inlat, const double inlon, const double angleOfRot, const double southPoleLat, const double southPoleLon, double* outlat, double* outlon);
double geographic_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2);
double geographic_distance_ellipsoid(double major, double minor, double lon1, double lat1, double lon2, double lat2);
/* grib_handle.c */
grib_section* grib_section_create(grib_handle* h, grib_accessor* owner);
@ -1376,8 +1378,16 @@ int grib_nearest_find(grib_nearest* nearest, const grib_handle* h, double inlat,
int grib_nearest_init(grib_nearest* i, grib_handle* h, grib_arguments* args);
int grib_nearest_delete(grib_nearest* i);
void grib_binary_search(double xx[], const unsigned long n, double x, int* ju, int* jl);
double grib_nearest_distance(double radius, double lon1, double lat1, double lon2, double lat2);
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, const char* radius_keyname,
const char* Ni_keyname, const char* Nj_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.c */
grib_nearest* grib_nearest_factory(grib_handle* h, grib_arguments* args);

2
src/grib_bits_any_endian.c
View File

@ -420,7 +420,7 @@ int grib_encode_size_tb(unsigned char* p, size_t val, long* bitp, long nb)
#if OMP_PACKING
#include "grib_bits_any_endian_omp.c"
#elif VECTOR
#include "grib_bits_any_endian_vector.c"
#include "grib_bits_any_endian_vector.c" /* Experimental */
#else
#include "grib_bits_any_endian_simple.c"
#endif

12
src/grib_bits_any_endian_vector.c
View File

@ -10,7 +10,7 @@
/***************************************************************************
* Enrico Fucile - 19.06.2007 *
* *
* EXPERIMENTAL CODE - NOT FULLY TESTED *
***************************************************************************/
int grib_decode_long_array(const unsigned char* p, long* bitp, long bitsPerValue,
@ -19,7 +19,8 @@ int grib_decode_long_array(const unsigned char* p, long* bitp, long bitsPerValue
long i = 0;
unsigned long lvalue = 0;
if (bitsPerValue % 8) {
/* SUP-3196: Thanks to Daniel Tameling */
if (bitsPerValue % 8 || (*bitp & 7)) {
int j = 0;
for (i = 0; i < n_vals; i++) {
lvalue = 0;
@ -54,7 +55,6 @@ int grib_decode_long_array(const unsigned char* p, long* bitp, long bitsPerValue
return 0;
}
int grib_decode_double_array(const unsigned char* p, long* bitp, long bitsPerValue,
double reference_value, double s, double d,
size_t n_vals, double* val)
@ -96,12 +96,14 @@ int grib_decode_double_array(const unsigned char* p, long* bitp, long bitsPerVal
return 0;
}
int grib_decode_double_array_complex(const unsigned char* p, long* bitp, long nbits, double reference_value, double s, double* d, size_t size, double* val)
int grib_decode_double_array_complex(const unsigned char* p, long* bitp, long nbits, double reference_value,
double s, double* d, size_t size, double* val)
{
return GRIB_NOT_IMPLEMENTED;
}
int grib_encode_double_array(size_t n_vals, const double* val, long bits_per_value, double reference_value, double d, double divisor, unsigned char* p, long* off)
int grib_encode_double_array(size_t n_vals, const double* val, long bits_per_value, double reference_value,
double d, double divisor, unsigned char* p, long* off)
{
size_t i = 0;
unsigned long unsigned_val = 0;

2
src/grib_bits_fast_big_endian.c
View File

@ -186,7 +186,7 @@ int grib_encode_unsigned_longb(unsigned char* p, unsigned long val, long* bitp,
}
#if VECTOR
#include "grib_bits_fast_big_endian_vector.c"
#include "grib_bits_fast_big_endian_vector.c" /* Experimental */
#elif OMP
#include "grib_bits_fast_big_endian_omp.c"
#else

70
src/grib_geography.c
View File

@ -4140,3 +4140,73 @@ void unrotate(const double inlat, const double inlon,
*outlat = ret_lat;
*outlon = ret_lon;
}
#define RADIAN(x) ((x)*acos(0.0) / 90.0)
/* radius is in km, angles in degrees */
/* Spherical Law of Cosines */
double geographic_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2)
{
double rlat1 = RADIAN(lat1);
double rlat2 = RADIAN(lat2);
double rlon1 = lon1;
double rlon2 = lon2;
double a;
if (lat1 == lat2 && lon1 == lon2) {
return 0.0; /* the two points are identical */
}
if (rlon1 >= 360) rlon1 -= 360.0;
rlon1 = RADIAN(rlon1);
if (rlon2 >= 360) rlon2 -= 360.0;
rlon2 = RADIAN(rlon2);
a = sin(rlat1) * sin(rlat2) + cos(rlat1) * cos(rlat2) * cos(rlon2 - rlon1);
DebugAssert(a >= -1 && a <= 1);
return radius * acos(a);
}
/* major and minor axes in km, angles in degrees */
double geographic_distance_ellipsoid(double major, double minor, double lon1, double lat1, double lon2, double lat2)
{
/* Lambert's formula */
double rlat1 = RADIAN(lat1);
double rlat2 = RADIAN(lat2);
double rlon1 = RADIAN(lon1);
double rlon2 = RADIAN(lon2);
double deltaLat = rlat2 - rlat1;
double deltaLon = rlon2 - rlon1;
double sinDlat = sin(deltaLat/2.0);
double sinDlon = sin(deltaLon/2.0);
double sin2Dlat = sinDlat*sinDlat;
double sin2Dlon = sinDlon*sinDlon;
double a = sin2Dlat + cos(rlat1) * cos(rlat2) * sin2Dlon;
double c = 2 * atan2(sqrt(a), sqrt(1.0-a));
double f = (major - minor)/major; /*flattening*/
double latr1 = atan( (1.0-f)*tan(rlat1) ); /*Reduced latitude1*/
double latr2 = atan( (1.0-f)*tan(rlat2) ); /*Reduced latitude2*/
double P = (latr1+latr2)/2;
double Q = (latr2-latr1)/2;
double sinP = sin(P);
double sin2P = sinP*sinP;
double cosQ = cos(Q);
double cos2Q = cosQ*cosQ;
double cosc2 = cos(c/2);
double cos2c2 = cosc2*cosc2;
double sinQ = sin(Q);
double sin2Q = sinQ*sinQ;
double cosP = cos(P);
double cos2P = cosP*cosP;
double sinc2 = sin(c/2);
double sin2c2 = sinc2*sinc2;
double X = (c - sin(c))* sin2P * cos2Q / cos2c2;
double Y = (c + sin(c))*sin2Q*cos2P/sin2c2;
double dist = major * (c - f*(X+Y)/2);
return dist;
}

15
src/grib_io.c
View File

@ -90,6 +90,7 @@ static int read_the_rest(reader* r, size_t message_length, unsigned char* tmp, i
size_t buffer_size;
size_t rest;
unsigned char* buffer;
grib_context* c = grib_context_get_default();
if (message_length == 0)
return GRIB_BUFFER_TOO_SMALL;
@ -109,12 +110,20 @@ static int read_the_rest(reader* r, size_t message_length, unsigned char* tmp, i
if ((r->read(r->read_data, buffer + already_read, rest, &err) != rest) || err) {
/*fprintf(stderr, "read_the_rest: r->read failed: %s\n", grib_get_error_message(err));*/
if (c->debug)
fprintf(stderr, "ECCODES DEBUG: read_the_rest: Read failed (Coded length=%lu, Already read=%d)\n",
message_length, already_read);
return err;
}
if (check7777 && !r->headers_only && (buffer[message_length - 4] != '7' || buffer[message_length - 3] != '7' || buffer[message_length - 2] != '7' || buffer[message_length - 1] != '7')) {
grib_context* c = grib_context_get_default();
grib_context_log(c, GRIB_LOG_DEBUG, "read_the_rest: No final 7777 at expected location (Coded length=%lu)", message_length);
if (check7777 && !r->headers_only &&
(buffer[message_length - 4] != '7' ||
buffer[message_length - 3] != '7' ||
buffer[message_length - 2] != '7' ||
buffer[message_length - 1] != '7'))
{
if (c->debug)
fprintf(stderr, "ECCODES DEBUG: read_the_rest: No final 7777 at expected location (Coded length=%lu)\n", message_length);
return GRIB_WRONG_LENGTH;
}

6
src/grib_iterator_class_lambert_conformal.c
View File

@ -360,9 +360,11 @@ static int init_oblate(grib_handle* h,
latRad = -M_PI_2;
}
lonRad = adjust_lon_radians(theta / ns + LoVInRadians);
if (i == 0 && j == 0) {
DebugAssert(fabs(latFirstInRadians - latRad) <= EPSILON);
}
latDeg = latRad * RAD2DEG; /* Convert to degrees */
lonDeg = lonRad * RAD2DEG;
lonDeg = normalise_longitude_in_degrees(lonDeg);
lonDeg = normalise_longitude_in_degrees(lonRad * RAD2DEG);
self->lons[index] = lonDeg;
self->lats[index] = latDeg;
}

8
src/grib_iterator_class_mercator.c
View File

@ -215,10 +215,12 @@ static int init_mercator(grib_handle* h,
grib_context_free(h->context, self->lons);
return err;
}
lonRad = adjust_lon_radians(orientationInRadians + _x / (earthMajorAxisInMetres * m1));
lonRad = adjust_lon_radians(orientationInRadians + _x / (earthMajorAxisInMetres * m1));
if (i == 0 && j == 0) {
DebugAssert(fabs(latFirstInRadians - latRad) <= EPSILON);
}
latDeg = latRad * RAD2DEG; /* Convert to degrees */
lonDeg = lonRad * RAD2DEG;
lonDeg = normalise_longitude_in_degrees(lonDeg);
lonDeg = normalise_longitude_in_degrees(lonRad * RAD2DEG);
self->lons[index] = lonDeg;
self->lats[index] = latDeg;
}

247
src/grib_nearest.c
View File

@ -114,35 +114,6 @@ void grib_binary_search(double xx[], const unsigned long n, double x,
}
}
#define RADIAN(x) ((x)*acos(0.0) / 90.0)
double grib_nearest_distance(double radius, double lon1, double lat1, double lon2, double lat2)
{
double rlat1 = RADIAN(lat1);
double rlat2 = RADIAN(lat2);
double rlon1 = lon1;
double rlon2 = lon2;
double a;
if (lat1 == lat2 && lon1 == lon2) {
return 0.0; /* the two points are identical */
}
if (rlon1 >= 360)
rlon1 -= 360.0;
rlon1 = RADIAN(rlon1);
if (rlon2 >= 360)
rlon2 -= 360.0;
rlon2 = RADIAN(rlon2);
a = sin(rlat1) * sin(rlat2) + cos(rlat1) * cos(rlat2) * cos(rlon2 - rlon1);
if (a > 1 || a < -1)
a = (int)a;
return radius * acos(a);
}
int grib_nearest_find_multiple(
const grib_handle* h, int is_lsm,
const double* inlats, const double* inlons, long npoints,
@ -241,3 +212,221 @@ 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,
const char* radius_keyname,
const char* Ni_keyname,
const char* Nj_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, i = 0;
size_t nvalues = 0, nneighbours = 0;
double radiusInMetres, 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;
nearest->values_count = nvalues;
/* We need the radius to calculate the nearest distance. For an oblate earth
approximate this using the average of the semimajor and semiminor axes */
if ((ret = grib_get_double(h, radius_keyname, &radiusInMetres)) == GRIB_SUCCESS &&
!grib_is_missing(h, radius_keyname, &ret)) {
radiusInKm = radiusInMetres / 1000.0;
}
else {
double minor = 0, major = 0;
if ((ret = grib_get_double_internal(h, "earthMinorAxisInMetres", &minor)) != GRIB_SUCCESS) return ret;
if ((ret = grib_get_double_internal(h, "earthMajorAxisInMetres", &major)) != GRIB_SUCCESS) return ret;
if (grib_is_missing(h, "earthMinorAxisInMetres", &ret)) return GRIB_GEOCALCULUS_PROBLEM;
if (grib_is_missing(h, "earthMajorAxisInMetres", &ret)) return GRIB_GEOCALCULUS_PROBLEM;
radiusInMetres = (major + minor) / 2;
radiusInKm = radiusInMetres / 1000.0;
}
neighbours = (PointStore*)grib_context_malloc(nearest->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;
int idx_upper = 0, idx_lower = 0;
double lat1 = 0, lat2 = 0; /* inlat will be between these */
const double LAT_DELTA = 10.0; /* in degrees */
if (grib_is_missing(h, Ni_keyname, &ret)) {
grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Ni_keyname);
return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
}
if (grib_is_missing(h, Nj_keyname, &ret)) {
grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Nj_keyname);
return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
}
*out_lons_count = nvalues; /* Maybe overestimate but safe */
*out_lats_count = nvalues;
if (*out_lats)
grib_context_free(nearest->context, *out_lats);
*out_lats = (double*)grib_context_malloc(nearest->context, nvalues * sizeof(double));
if (!*out_lats)
return GRIB_OUT_OF_MEMORY;
if (*out_lons)
grib_context_free(nearest->context, *out_lons);
*out_lons = (double*)grib_context_malloc(nearest->context, nvalues * sizeof(double));
if (!*out_lons)
return GRIB_OUT_OF_MEMORY;
iter = grib_iterator_new(h, 0, &ret);
if (ret)
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 = 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);
}
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(nearest->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;
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;
}

3
src/grib_nearest_class.h
View File

@ -1,8 +1,11 @@
/* 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_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_sh;
extern grib_nearest_class* grib_nearest_class_space_view;

136
src/grib_nearest_class_lambert_azimuthal_equal_area.c Normal file
View File

@ -0,0 +1,136 @@
/*
* (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->radius,
self->Ni,
self->Nj,
&(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;
if (self->lats) grib_context_free(nearest->context, self->lats);
if (self->lons) grib_context_free(nearest->context, self->lons);
if (self->i) grib_context_free(nearest->context, self->i);
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;
}

225
src/grib_nearest_class_lambert_conformal.c
View File

@ -100,224 +100,37 @@ static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args)
return 0;
}
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;
}
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;
int ret = 0, i = 0;
size_t nvalues = 0;
long iradius;
double radius;
grib_iterator* iter = NULL;
double lat = 0, lon = 0;
return grib_nearest_find_generic(
nearest, h, inlat, inlon, flags, /* inputs */
/* array of candidates for nearest neighbours */
PointStore* neighbours = NULL;
self->values_key, /* outputs to set the 'self' object */
self->radius,
self->Ni,
self->Nj,
&(self->lats),
&(self->lats_count),
&(self->lons),
&(self->lons_count),
&(self->distances),
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 (grib_is_missing(h, self->radius, &ret)) {
grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->radius);
return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
}
if ((ret = grib_get_long(h, self->radius, &iradius)) != GRIB_SUCCESS)
return ret;
radius = ((double)iradius) / 1000.0;
neighbours = (PointStore*)grib_context_malloc(nearest->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;
int idx_upper = 0, idx_lower = 0;
double lat1 = 0, lat2 = 0; /* inlat will be between these */
double dist = 0;
const double LAT_DELTA = 10.0; /* in degrees */
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;
}
self->lons_count = nvalues; /* Maybe overestimate but safe */
self->lats_count = nvalues;
if (self->lats)
grib_context_free(nearest->context, self->lats);
self->lats = (double*)grib_context_malloc(nearest->context, nvalues * 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, nvalues * sizeof(double));
if (!self->lons)
return GRIB_OUT_OF_MEMORY;
iter = grib_iterator_new(h, 0, &ret);
if (ret)
return ret;
/* First pass: collect all latitudes and longitudes */
while (grib_iterator_next(iter, &lat, &lon, &the_value)) {
++the_index;
Assert(ilat < self->lats_count);
Assert(ilon < self->lons_count);
self->lats[ilat++] = lat;
self->lons[ilon++] = lon;
}
/* See between which 2 latitudes our point lies */
qsort(self->lats, nvalues, sizeof(double), &compare_doubles_ascending);
grib_binary_search(self->lats, self->lats_count - 1, inlat, &idx_upper, &idx_lower);
lat2 = self->lats[idx_upper];
lat1 = self->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 {
dist = grib_nearest_distance(radius, 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 = the_index;
neighbours[i].m_lat = lat;
neighbours[i].m_lon = lon;
neighbours[i].m_value = the_value;
i++;
}
++the_index;
}
/* Sort the candidate neighbours in ascending order of distance */
/* The first 4 entries will now be the closest 4 neighbours */
qsort(neighbours, nvalues, sizeof(PointStore), &compare_points);
grib_iterator_delete(iter);
}
nearest->h = h;
/* Sanity check for sorting */
#ifdef DEBUG
for (i = 0; i < nvalues - 1; ++i) {
Assert(neighbours[i].m_dist <= neighbours[i + 1].m_dist);
}
#endif
/* GRIB_NEAREST_SAME_XXX not yet implemented */
if (!self->distances) {
self->distances = (double*)grib_context_malloc(nearest->context, 4 * sizeof(double));
}
self->distances[0] = neighbours[0].m_dist;
self->distances[1] = neighbours[1].m_dist;
self->distances[2] = neighbours[2].m_dist;
self->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;
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;
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;
if (self->lats)
grib_context_free(nearest->context, self->lats);
if (self->lons)
grib_context_free(nearest->context, self->lons);
if (self->i)
grib_context_free(nearest->context, self->i);
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);
if (self->lats) grib_context_free(nearest->context, self->lats);
if (self->lons) grib_context_free(nearest->context, self->lons);
if (self->i) grib_context_free(nearest->context, self->i);
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;
}

2
src/grib_nearest_class_latlon_reduced.c
View File

@ -346,7 +346,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
kk = 0;
for (jj = 0; jj < 2; jj++) {
for (ii = 0; ii < 2; ii++) {
self->distances[kk] = grib_nearest_distance(radius, inlon, inlat,
self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
self->lons[self->k[kk]], self->lats[self->j[jj]]);
kk++;
}

138
src/grib_nearest_class_mercator.c Normal file
View File

@ -0,0 +1,138 @@
/*
* (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->radius,
self->Ni,
self->Nj,
&(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;
if (self->lats) grib_context_free(nearest->context, self->lats);
if (self->lons) grib_context_free(nearest->context, self->lons);
if (self->i) grib_context_free(nearest->context, self->i);
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;
}

224
src/grib_nearest_class_polar_stereographic.c
View File

@ -97,53 +97,7 @@ static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args)
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 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;
return GRIB_SUCCESS;
}
static int find(grib_nearest* nearest, grib_handle* h,
@ -152,169 +106,31 @@ static int find(grib_nearest* nearest, grib_handle* h,
double* values, double* distances, int* indexes, size_t* len)
{
grib_nearest_polar_stereographic* self = (grib_nearest_polar_stereographic*)nearest;
int ret = 0, i = 0;
size_t nvalues = 0;
long iradius;
double radius;
grib_iterator* iter = NULL;
double lat = 0, lon = 0;
return grib_nearest_find_generic(
nearest, h, inlat, inlon, flags, /* inputs */
/* array of candidates for nearest neighbours */
PointStore* neighbours = NULL;
self->values_key, /* outputs to set the 'self' object */
self->radius,
self->Ni,
self->Nj,
&(self->lats),
&(self->lats_count),
&(self->lons),
&(self->lons_count),
&(self->distances),
inlon = normalise_longitude_in_degrees(inlon);
if ((ret = grib_get_size(h, self->values_key, &nvalues)) != GRIB_SUCCESS)
return ret;
nearest->values_count = nvalues;
if (grib_is_missing(h, self->radius, &ret)) {
grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->radius);
return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
}
if ((ret = grib_get_long(h, self->radius, &iradius)) != GRIB_SUCCESS)
return ret;
radius = ((double)iradius) / 1000.0;
neighbours = (PointStore*)grib_context_malloc(nearest->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;
int idx_upper = 0, idx_lower = 0;
double lat1 = 0, lat2 = 0; /* inlat will be between these */
double dist = 0;
const double LAT_DELTA = 10.0; /* in degrees */
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;
}
self->lons_count = nvalues; /* Maybe overestimate but safe */
self->lats_count = nvalues;
if (self->lats)
grib_context_free(nearest->context, self->lats);
self->lats = (double*)grib_context_malloc(nearest->context, nvalues * 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, nvalues * sizeof(double));
if (!self->lons)
return GRIB_OUT_OF_MEMORY;
iter = grib_iterator_new(h, 0, &ret);
if (ret)
return ret;
/* First pass: collect all latitudes and longitudes */
while (grib_iterator_next(iter, &lat, &lon, &the_value)) {
++the_index;
Assert(ilat < self->lats_count);
Assert(ilon < self->lons_count);
self->lats[ilat++] = lat;
self->lons[ilon++] = lon;
}
/* See between which 2 latitudes our point lies */
qsort(self->lats, nvalues, sizeof(double), &compare_doubles_ascending);
grib_binary_search(self->lats, self->lats_count - 1, inlat, &idx_upper, &idx_lower);
lat2 = self->lats[idx_upper];
lat1 = self->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 {
dist = grib_nearest_distance(radius, 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 = the_index;
neighbours[i].m_lat = lat;
neighbours[i].m_lon = lon;
neighbours[i].m_value = the_value;
i++;
}
++the_index;
}
/* Sort the candidate neighbours in ascending order of distance */
/* The first 4 entries will now be the closest 4 neighbours */
qsort(neighbours, nvalues, sizeof(PointStore), &compare_points);
grib_iterator_delete(iter);
}
nearest->h = h;
/* Sanity check for sorting */
#ifdef DEBUG
for (i = 0; i < nvalues - 1; ++i) {
Assert(neighbours[i].m_dist <= neighbours[i + 1].m_dist);
}
#endif
/* GRIB_NEAREST_SAME_XXX not yet implemented */
if (!self->distances) {
self->distances = (double*)grib_context_malloc(nearest->context, 4 * sizeof(double));
}
self->distances[0] = neighbours[0].m_dist;
self->distances[1] = neighbours[1].m_dist;
self->distances[2] = neighbours[2].m_dist;
self->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;
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;
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;
if (self->lats)
grib_context_free(nearest->context, self->lats);
if (self->lons)
grib_context_free(nearest->context, self->lons);
if (self->i)
grib_context_free(nearest->context, self->i);
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);
if (self->lats) grib_context_free(nearest->context, self->lats);
if (self->lons) grib_context_free(nearest->context, self->lons);
if (self->i) grib_context_free(nearest->context, self->i);
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;
}

4
src/grib_nearest_class_reduced.c
View File

@ -409,7 +409,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
kk = 0;
for (jj = 0; jj < 2; jj++) {
for (ii = 0; ii < 2; ii++) {
self->distances[kk] = grib_nearest_distance(radius, inlon, inlat,
self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
self->lons[self->k[kk]], self->lats[self->j[jj]]);
kk++;
}
@ -548,7 +548,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
kk = 0;
for (ii = 0; ii < 2; ii++) {
for (jj = 0; jj < 2; jj++) {
self->distances[kk] = grib_nearest_distance(radius, inlon, inlat,
self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
self->lons[self->k[kk]], self->lats[self->j[jj]]);
kk++;
}

4
src/grib_nearest_class_regular.c
View File

@ -193,7 +193,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
for (ii=0;ii<2;ii++) {
for (jj=0;jj<2;jj++) {
self->k[kk]=self->i[ii]+self->lons_count*self->j[jj]-1;
self->distances[kk]=grib_nearest_distance(radius,inlon,inlat,
self->distances[kk]=geographic_distance_spherical(radius,inlon,inlat,
self->lons[self->i[ii]],self->lats[self->j[jj]]);
kk++;
}
@ -409,7 +409,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
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] = grib_nearest_distance(radius, inlon, inlat,
self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
self->lons[self->i[ii]], self->lats[self->j[jj]]);
kk++;
}

136
src/grib_nearest_class_space_view.c Normal file
View File

@ -0,0 +1,136 @@
/*
* (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->radius,
self->Ni,
self->Nj,
&(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;
if (self->lats) grib_context_free(nearest->context, self->lats);
if (self->lons) grib_context_free(nearest->context, self->lons);
if (self->i) grib_context_free(nearest->context, self->i);
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;
}

3
src/grib_nearest_factory.h
View File

@ -1,8 +1,11 @@
/* This file is automatically generated by ./make_class.pl, do not edit */
{ "gen", &grib_nearest_class_gen, },
{ "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, },
{ "sh", &grib_nearest_class_sh, },
{ "space_view", &grib_nearest_class_space_view, },

1
tests/CMakeLists.txt
View File

@ -69,6 +69,7 @@ list(APPEND tests_no_data_reqd
grib_sh_imag
pseudo_diag
grib_lambert_conformal
grib_space_view
grib_g1fcperiod)
# These tests do require data downloads

1
tests/Makefile.am
View File

@ -57,6 +57,7 @@ TESTS = definitions.sh \
grib_padding.sh \
grib_lamb_az_eq_area.sh \
grib_lambert_conformal.sh \
grib_space_view.sh \
grib_mercator.sh \
grib_to_netcdf.sh \
grib_dump_debug.sh \

4
tests/grib_lamb_az_eq_area.sh
View File

@ -55,5 +55,9 @@ diff $DATA_OUTFILE $REF_FILE
grib_check_key_equals $GRIB_OUTFILE standardParallelInDegrees,centralLongitudeInDegrees '48 9'
grib_check_key_equals $GRIB_OUTFILE xDirectionGridLengthInMetres,yDirectionGridLengthInMetres '5000 5000'
# Nearest
${tools_dir}/grib_ls -l 67,-33,1 $GRIB_OUTFILE
# Clean up
rm -f $FILTER_FILE $GRIB_OUTFILE $DATA_OUTFILE

8
tests/grib_lambert_conformal.sh
View File

@ -52,7 +52,6 @@ ${tools_dir}/grib_ls -l 50,0 $tempGrib
cat > $tempFilter <<EOF
set gridType="lambert";
set numberOfDataPoints=294000;
set shapeOfTheEarth=6;
set Nx=588;
set Ny=500;
set latitudeOfFirstGridPoint=40442000;
@ -63,7 +62,7 @@ cat > $tempFilter <<EOF
set Dy=2499000;
set Latin1=46401000;
set Latin2=46401000;
set shapeOfTheEarth=2;
set shapeOfTheEarth=2; # oblate earth
set numberOfValues=294000;
write;
EOF
@ -76,10 +75,11 @@ if [ ! -f "$tempGrib" ]; then
fi
grib_check_key_equals $tempGrib 'earthIsOblate,earthMinorAxisInMetres,earthMajorAxisInMetres' '1 6356775 6378160'
# Invoke Geoiterator on the newly created GRIB file
# Invoke Geoiterator on the oblate lambert GRIB
${tools_dir}/grib_get_data $tempGrib > $tempOut
# Nearest neighbour on the oblate lambert GRIB
${tools_dir}/grib_ls -l 40.44,353.56 $tempGrib
# Clean up

9
tests/grib_local.sh
View File

@ -132,6 +132,15 @@ ${tools_dir}/grib_set -s localDefinitionNumber=49,type=35 $sample_g1 $temp
grib_check_key_equals $temp 'perturbationNumber,numberOfForecastsInEnsemble' '0 0'
# Local Definition 39 and type 'eme' for GRIB2
# --------------------------------------------
${tools_dir}/grib_set -s \
setLocalDefinition=1,localDefinitionNumber=39,type=eme,stream=elda,componentIndex=11,offsetToEndOf4DvarWindow=12 \
$sample_g2 $temp
grib_check_key_equals $temp 'mars.number,mars.anoffset' '11 12'
${tools_dir}/grib_set -s setLocalDefinition=1,localDefinitionNumber=39,type=eme,stream=enda,componentIndex=11 $sample_g2 $temp
grib_check_key_equals $temp 'mars.number' '11'
# Local Definition 18 (list of ascii keys)
# ----------------------------------------

4
tests/grib_mercator.sh
View File

@ -30,5 +30,9 @@ EOF
${tools_dir}/grib_filter $tempFilter $input
# Nearest function
${tools_dir}/grib_ls -l 19,-97,1 $input > $tempOut
grep -q "Point chosen #1 index=618" $tempOut
# Clean up
rm -f $tempFilter $tempOut

48
tests/grib_space_view.sh Executable file
View File

@ -0,0 +1,48 @@
#!/bin/sh
# (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.sh
# Define a common label for all the tmp files
label="grib_space_view_test"
tempFilter="temp.${label}.filt"
tempGrib="temp.${label}.grib"
tempOut="temp.${label}.out"
input=$ECCODES_SAMPLES_PATH/GRIB2.tmpl
# Create a filter
cat > $tempFilter <<EOF
set gridType="space_view";
set Nx=1900;
set Ny=900;
set dx=3622;
set dy=3610;
set Xp=764000;
set Yp=1774000;
set Nr=6610700;
set numberOfDataPoints=1710000; # 1900 x 900
set numberOfValues=1710000;
write;
EOF
# Use filter on input to create a new GRIB
${tools_dir}/grib_filter -o $tempGrib $tempFilter $input
if [ ! -f "$tempGrib" ]; then
echo 'Failed to create output GRIB from filter' >&2
exit 1
fi
# Invoke Geoiterator on the newly created GRIB file
${tools_dir}/grib_get_data $tempGrib > $tempOut
${tools_dir}/grib_ls -l 50,0 $tempGrib
# Clean up
rm -f $tempFilter $tempGrib $tempOut

27
tools/grib_ls.c
View File

@ -156,14 +156,14 @@ int grib_tool_init(grib_runtime_options* options)
options->latlon_idx = -1;
max = options->distances[0];
for (i = 0; i < 4; i++)
for (i = 0; i < LATLON_SIZE; i++)
if (max < options->distances[i]) {
max = options->distances[i];
}
min = max;
min_overall = max;
/* See GRIB-213 */
for (i = 0; i < 4; i++) {
for (i = 0; i < LATLON_SIZE; i++) {
if (min_overall >= options->distances[i]) { /* find overall min and index ignoring mask */
min_overall = options->distances[i];
idx_overall = i;
@ -184,7 +184,7 @@ int grib_tool_init(grib_runtime_options* options)
if (options->latlon_idx < 0) {
min = 0;
options->latlon_idx = 0;
for (i = 1; i < 4; i++)
for (i = 1; i < LATLON_SIZE; i++)
if (min > options->distances[i]) {
min = options->distances[i];
options->latlon_idx = i;
@ -307,7 +307,7 @@ int grib_tool_new_handle_action(grib_runtime_options* options, grib_handle* h)
if (!options->latlon_mask) {
min = options->distances[0];
options->latlon_idx = 0;
for (i = 1; i < 4; i++) {
for (i = 1; i < LATLON_SIZE; i++) {
if (min > options->distances[i]) {
min = options->distances[i];
options->latlon_idx = i;
@ -330,7 +330,7 @@ int grib_tool_new_handle_action(grib_runtime_options* options, grib_handle* h)
else
printf("\"nearest\"");
printf("\n, \"neighbours\" : ");
for (i = 0; i < 4; i++) {
for (i = 0; i < LATLON_SIZE; i++) {
printf("%s", s);
len = MAX_STRING_LEN;
printf(
@ -398,15 +398,16 @@ int grib_tool_finalise_action(grib_runtime_options* options)
int i = 0;
if (options->latlon && options->verbose) {
printf("Input Point: latitude=%.2f longitude=%.2f\n", lat, lon);
printf("Grid Point chosen #%d index=%d latitude=%.2f longitude=%.2f distance=%.2f (Km)\n",
options->latlon_idx + 1, (int)options->indexes[options->latlon_idx],
options->lats[options->latlon_idx],
options->lons[options->latlon_idx],
options->distances[options->latlon_idx]);
if (options->latlon_idx >= 0 && options->latlon_idx < LATLON_SIZE) {
printf("Grid Point chosen #%d index=%d latitude=%.2f longitude=%.2f distance=%.2f (Km)\n",
options->latlon_idx + 1, (int)options->indexes[options->latlon_idx],
options->lats[options->latlon_idx],
options->lons[options->latlon_idx],
options->distances[options->latlon_idx]);
}
if (options->latlon_mask) {
printf("Mask values:\n");
for (i = 0; i < 4; i++) {
for (i = 0; i < LATLON_SIZE; i++) {
printf("- %d - index=%d latitude=%.2f longitude=%.2f distance=%.2f (Km) value=%.2f\n",
i + 1, (int)options->indexes[i], options->lats[i], options->lons[i],
options->distances[i], options->mask_values[i]);
@ -414,7 +415,7 @@ int grib_tool_finalise_action(grib_runtime_options* options)
}
else {
printf("Other grid Points\n");
for (i = 0; i < 4; i++) {
for (i = 0; i < LATLON_SIZE; i++) {
printf("- %d - index=%d latitude=%.2f longitude=%.2f distance=%.2f (Km)\n",
i + 1, (int)options->indexes[i], options->lats[i], options->lons[i],
options->distances[i]);

2
tools/grib_tools.c
View File

@ -1170,7 +1170,7 @@ void grib_print_key_values(grib_runtime_options* options, grib_handle* h)
if (options->latlon) {
if (options->latlon_mode == 4) {
int ii = 0;
for (ii = 0; ii < 4; ii++) {
for (ii = 0; ii < LATLON_SIZE; ii++) {
fprintf(dump_file, options->format, options->values[ii]);
fprintf(dump_file, " ");
}

13
tools/grib_tools.h
View File

@ -37,6 +37,7 @@
#define MAX_STRING_LEN 512
#define MAX_FAILED 1024
#define MAX_CONSTRAINT_VALUES 500
#define LATLON_SIZE 4 /* nearest */
#define MODE_GRIB 0
#define MODE_GTS 1
@ -140,15 +141,15 @@ typedef struct grib_runtime_options
int gts;
char* orderby;
char* latlon;
double lats[4];
double lons[4];
double values[4];
double distances[4];
int indexes[4];
double lats[LATLON_SIZE];
double lons[LATLON_SIZE];
double values[LATLON_SIZE];
double distances[LATLON_SIZE];
int indexes[LATLON_SIZE];
int latlon_mode;
char* latlon_mask;
int latlon_idx;
double mask_values[4];
double mask_values[LATLON_SIZE];
int index;
int index_on;
double constant;

2
version.sh
View File

@ -3,7 +3,7 @@ PACKAGE_NAME='eccodes'
# Package version
ECCODES_MAJOR_VERSION=2
ECCODES_MINOR_VERSION=18
ECCODES_MINOR_VERSION=19
ECCODES_REVISION_VERSION=0
ECCODES_CURRENT=1

3
windows/msvc/grib_api_lib/grib_api_lib.vcxproj
View File

@ -511,12 +511,15 @@
<ClCompile Include="..\..\..\src\grib_nearest.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_gen.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_lambert_azimuthal_equal_area.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_lambert_conformal.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_latlon_reduced.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_mercator.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_polar_stereographic.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_reduced.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_regular.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_sh.c" />
<ClCompile Include="..\..\..\src\grib_nearest_class_space_view.c" />
<ClCompile Include="..\..\..\src\grib_oarray.c" />
<ClCompile Include="..\..\..\src\grib_openjpeg_encoding.c" />
<ClCompile Include="..\..\..\src\grib_optimize_decimal_factor.c" />