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

Merge branch 'develop' of ssh://git.ecmwf.int:7999/eccodes/eccodes into develop

This commit is contained in:
Sebastien Villaume 2021-11-04 16:39:07 +00:00
parent 3f3a504b45 040ab8d66d
commit fee4f7fe3b
13 changed files with 1491 additions and 1222 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
definitions/grib2/localConcepts/kwbc/name.def
View File

@ -1625,30 +1625,6 @@
parameterCategory = 19 ; parameterCategory = 19 ;
parameterNumber = 234 ; parameterNumber = 234 ;
} }
#Ventilation Rate
'Ventilation Rate' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 224 ;
}
#Convective inhibition
'Convective inhibition' = {
discipline = 0 ;
parameterCategory = 7 ;
parameterNumber = 7 ;
}
#Total Cloud Cover
'Total Cloud Cover' = {
discipline = 0 ;
parameterCategory = 6 ;
parameterNumber = 1 ;
}
#Total Precipitation
'Total Precipitation' = {
discipline = 0 ;
parameterCategory = 1 ;
parameterNumber = 8 ;
}
#Simulated Brightness Temperature for AMSRE on Aqua, Channel 10 #Simulated Brightness Temperature for AMSRE on Aqua, Channel 10
'Simulated Brightness Temperature for AMSRE on Aqua, Channel 10' = { 'Simulated Brightness Temperature for AMSRE on Aqua, Channel 10' = {
discipline = 3 ; discipline = 3 ;
@ -1925,3 +1901,27 @@
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 35 ; parameterNumber = 35 ;
} }
#Ventilation Rate
'Ventilation Rate' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 224 ;
}
#Convective inhibition
'Convective inhibition' = {
discipline = 0 ;
parameterCategory = 7 ;
parameterNumber = 7 ;
}
#Total Cloud Cover
'Total Cloud Cover' = {
discipline = 0 ;
parameterCategory = 6 ;
parameterNumber = 1 ;
}
#Total Precipitation
'Total Precipitation' = {
discipline = 0 ;
parameterCategory = 1 ;
parameterNumber = 8 ;
}

48
definitions/grib2/localConcepts/kwbc/paramId.def
View File

@ -1625,30 +1625,6 @@
parameterCategory = 19 ; parameterCategory = 19 ;
parameterNumber = 234 ; parameterNumber = 234 ;
} }
#Ventilation Rate
'7001353' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 224 ;
}
#Convective inhibition
'228001' = {
discipline = 0 ;
parameterCategory = 7 ;
parameterNumber = 7 ;
}
#Total Cloud Cover
'228164' = {
discipline = 0 ;
parameterCategory = 6 ;
parameterNumber = 1 ;
}
#Total Precipitation
'228228' = {
discipline = 0 ;
parameterCategory = 1 ;
parameterNumber = 8 ;
}
#Simulated Brightness Temperature for AMSRE on Aqua, Channel 10 #Simulated Brightness Temperature for AMSRE on Aqua, Channel 10
'7001294' = { '7001294' = {
discipline = 3 ; discipline = 3 ;
@ -1925,3 +1901,27 @@
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 35 ; parameterNumber = 35 ;
} }
#Ventilation Rate
'7001353' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 224 ;
}
#Convective inhibition
'228001' = {
discipline = 0 ;
parameterCategory = 7 ;
parameterNumber = 7 ;
}
#Total Cloud Cover
'228164' = {
discipline = 0 ;
parameterCategory = 6 ;
parameterNumber = 1 ;
}
#Total Precipitation
'228228' = {
discipline = 0 ;
parameterCategory = 1 ;
parameterNumber = 8 ;
}

48
definitions/grib2/localConcepts/kwbc/shortName.def
View File

@ -1625,30 +1625,6 @@
parameterCategory = 19 ; parameterCategory = 19 ;
parameterNumber = 234 ; parameterNumber = 234 ;
} }
#Ventilation Rate
'VRATE' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 224 ;
}
#Convective inhibition
'cin' = {
discipline = 0 ;
parameterCategory = 7 ;
parameterNumber = 7 ;
}
#Total Cloud Cover
'tcc' = {
discipline = 0 ;
parameterCategory = 6 ;
parameterNumber = 1 ;
}
#Total Precipitation
'tp' = {
discipline = 0 ;
parameterCategory = 1 ;
parameterNumber = 8 ;
}
#Simulated Brightness Temperature for AMSRE on Aqua, Channel 10 #Simulated Brightness Temperature for AMSRE on Aqua, Channel 10
'AMSRE10' = { 'AMSRE10' = {
discipline = 3 ; discipline = 3 ;
@ -1925,3 +1901,27 @@
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 35 ; parameterNumber = 35 ;
} }
#Ventilation Rate
'VRATE' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 224 ;
}
#Convective inhibition
'cin' = {
discipline = 0 ;
parameterCategory = 7 ;
parameterNumber = 7 ;
}
#Total Cloud Cover
'tcc' = {
discipline = 0 ;
parameterCategory = 6 ;
parameterNumber = 1 ;
}
#Total Precipitation
'tp' = {
discipline = 0 ;
parameterCategory = 1 ;
parameterNumber = 8 ;
}

72
definitions/grib2/localConcepts/kwbc/units.def
View File

@ -1625,30 +1625,6 @@
parameterCategory = 19 ; parameterCategory = 19 ;
parameterNumber = 234 ; parameterNumber = 234 ;
} }
#Ventilation Rate
'm**2 s**-1' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 224 ;
}
#Convective inhibition
'J kg**-1' = {
discipline = 0 ;
parameterCategory = 7 ;
parameterNumber = 7 ;
}
#Total Cloud Cover
'%' = {
discipline = 0 ;
parameterCategory = 6 ;
parameterNumber = 1 ;
}
#Total Precipitation
'kg m**-2' = {
discipline = 0 ;
parameterCategory = 1 ;
parameterNumber = 8 ;
}
#Simulated Brightness Temperature for AMSRE on Aqua, Channel 10 #Simulated Brightness Temperature for AMSRE on Aqua, Channel 10
'K' = { 'K' = {
discipline = 3 ; discipline = 3 ;
@ -1854,74 +1830,98 @@
parameterNumber = 38 ; parameterNumber = 38 ;
} }
#Simulated Reflectance Factor for ABI GOES-16, Band-1 #Simulated Reflectance Factor for ABI GOES-16, Band-1
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 14 ; parameterNumber = 14 ;
} }
#Simulated Reflectance Factor for ABI GOES-16, Band-2 #Simulated Reflectance Factor for ABI GOES-16, Band-2
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 15 ; parameterNumber = 15 ;
} }
#Simulated Reflectance Factor for ABI GOES-16, Band-3 #Simulated Reflectance Factor for ABI GOES-16, Band-3
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 16 ; parameterNumber = 16 ;
} }
#Simulated Reflectance Factor for ABI GOES-16, Band-4 #Simulated Reflectance Factor for ABI GOES-16, Band-4
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 17 ; parameterNumber = 17 ;
} }
#Simulated Reflectance Factor for ABI GOES-16, Band-5 #Simulated Reflectance Factor for ABI GOES-16, Band-5
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 18 ; parameterNumber = 18 ;
} }
#Simulated Reflectance Factor for ABI GOES-16, Band-6 #Simulated Reflectance Factor for ABI GOES-16, Band-6
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 19 ; parameterNumber = 19 ;
} }
#Simulated Reflectance Factor for ABI GOES-17, Band-1 #Simulated Reflectance Factor for ABI GOES-17, Band-1
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 30 ; parameterNumber = 30 ;
} }
#Simulated Reflectance Factor for ABI GOES-17, Band-2 #Simulated Reflectance Factor for ABI GOES-17, Band-2
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 31 ; parameterNumber = 31 ;
} }
#Simulated Reflectance Factor for ABI GOES-17, Band-3 #Simulated Reflectance Factor for ABI GOES-17, Band-3
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 32 ; parameterNumber = 32 ;
} }
#Simulated Reflectance Factor for ABI GOES-17, Band-4 #Simulated Reflectance Factor for ABI GOES-17, Band-4
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 33 ; parameterNumber = 33 ;
} }
#Simulated Reflectance Factor for ABI GOES-17, Band-5 #Simulated Reflectance Factor for ABI GOES-17, Band-5
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 34 ; parameterNumber = 34 ;
} }
#Simulated Reflectance Factor for ABI GOES-17, Band-6 #Simulated Reflectance Factor for ABI GOES-17, Band-6
'-' = { '~' = {
discipline = 3 ; discipline = 3 ;
parameterCategory = 192 ; parameterCategory = 192 ;
parameterNumber = 35 ; parameterNumber = 35 ;
} }
#Ventilation Rate
'm**2 s**-1' = {
discipline = 0 ;
parameterCategory = 2 ;
parameterNumber = 224 ;
}
#Convective inhibition
'J kg**-1' = {
discipline = 0 ;
parameterCategory = 7 ;
parameterNumber = 7 ;
}
#Total Cloud Cover
'%' = {
discipline = 0 ;
parameterCategory = 6 ;
parameterNumber = 1 ;
}
#Total Precipitation
'kg m**-2' = {
discipline = 0 ;
parameterCategory = 1 ;
parameterNumber = 8 ;
}

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

@ -55,7 +55,7 @@ include "grib2/template.3.scanning_mode.def";
iterator lambert_azimuthal_equal_area(numberOfPoints,missingValue,values, iterator lambert_azimuthal_equal_area(numberOfPoints,missingValue,values,
radius,Nx,Ny, radius,Nx,Ny,
latitudeOfFirstGridPointInDegrees,longitudeOfFirstGridPointInDegrees, latitudeOfFirstGridPointInDegrees,longitudeOfFirstGridPointInDegrees,
standardParallel,centralLongitude, standardParallelInDegrees,centralLongitudeInDegrees,
Dx,Dy, Dx,Dy,
iScansNegatively, jScansPositively, iScansNegatively, jScansPositively,
jPointsAreConsecutive, alternativeRowScanning); jPointsAreConsecutive, alternativeRowScanning);

4
src/grib_accessor_class_bufr_data_array.c
View File

@ -766,7 +766,7 @@ static int descriptor_get_min_max(bufr_descriptor* bd, long width, long referenc
double* minAllowed, double* maxAllowed) double* minAllowed, double* maxAllowed)
{ {
/* Maximum value is allowed to be the largest number (all bits 1) which means it's MISSING */ /* Maximum value is allowed to be the largest number (all bits 1) which means it's MISSING */
unsigned long max1 = (1UL << width) - 1; /* Highest value for number with 'width' bits */ size_t max1 = (1ULL << width) - 1; /* Highest value for number with 'width' bits */
DebugAssert(width > 0 && width < 64); DebugAssert(width > 0 && width < 64);
*maxAllowed = (max1 + reference) * factor; *maxAllowed = (max1 + reference) * factor;
@ -1100,7 +1100,7 @@ static double decode_double_value(grib_context* c, unsigned char* data, long* po
dval = GRIB_MISSING_DOUBLE; dval = GRIB_MISSING_DOUBLE;
} }
else { else {
dval = ((long)lval + modifiedReference) * modifiedFactor; dval = ((int64_t)lval + modifiedReference) * modifiedFactor;
} }
return dval; return dval;
} }

5
src/grib_api_prototypes.h
View File

@ -1,3 +1,6 @@
#ifdef ECCODES_ON_WINDOWS
#include <stdint.h>
#endif
/* action.c */ /* action.c */
void grib_dump(grib_action* a, FILE* f, int l); void grib_dump(grib_action* a, FILE* f, int l);
@ -1550,7 +1553,7 @@ int grib_optimize_decimal_factor(grib_accessor* a, const char* reference_value,
const char* grib_get_git_sha1(void); const char* grib_get_git_sha1(void);
/* grib_bits_any_endian.c */ /* grib_bits_any_endian.c */
int grib_is_all_bits_one(long val, long nbits); int grib_is_all_bits_one(int64_t val, long nbits);
int grib_encode_string(unsigned char* bitStream, long* bitOffset, size_t numberOfCharacters, const char* string); int grib_encode_string(unsigned char* bitStream, long* bitOffset, size_t numberOfCharacters, const char* string);
char* grib_decode_string(const unsigned char* bitStream, long* bitOffset, size_t numberOfCharacters, char* string); char* grib_decode_string(const unsigned char* bitStream, long* bitOffset, size_t numberOfCharacters, char* string);
unsigned long grib_decode_unsigned_long(const unsigned char* p, long* bitp, long nbits); unsigned long grib_decode_unsigned_long(const unsigned char* p, long* bitp, long nbits);

16
src/grib_bits_any_endian.c
View File

@ -13,6 +13,10 @@
* * * *
***************************************************************************/ ***************************************************************************/
#ifdef ECCODES_ON_WINDOWS
#include <stdint.h>
#endif
#if GRIB_PTHREADS #if GRIB_PTHREADS
static pthread_once_t once = PTHREAD_ONCE_INIT; static pthread_once_t once = PTHREAD_ONCE_INIT;
static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER; static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
@ -45,16 +49,16 @@ typedef struct bits_all_one_t
{ {
int inited; int inited;
int size; int size;
long v[128]; int64_t v[128];
} bits_all_one_t; } bits_all_one_t;
static bits_all_one_t bits_all_one = { 0, 0, {0,} }; static bits_all_one_t bits_all_one = { 0, 0, {0,} };
static void init_bits_all_one() static void init_bits_all_one()
{ {
int size = sizeof(long) * 8; int size = sizeof(int64_t) * 8;
long* v = 0; int64_t* v = 0;
unsigned long cmask = -1; uint64_t cmask = -1;
DebugAssert(!bits_all_one.inited); DebugAssert(!bits_all_one.inited);
bits_all_one.size = size; bits_all_one.size = size;
@ -78,7 +82,7 @@ static void init_bits_all_one_if_needed()
init_bits_all_one(); init_bits_all_one();
GRIB_MUTEX_UNLOCK(&mutex); GRIB_MUTEX_UNLOCK(&mutex);
} }
int grib_is_all_bits_one(long val, long nbits) int grib_is_all_bits_one(int64_t val, long nbits)
{ {
/*if (!bits_all_one.inited) init_bits_all_one();*/ /*if (!bits_all_one.inited) init_bits_all_one();*/
init_bits_all_one_if_needed(); init_bits_all_one_if_needed();
@ -301,7 +305,7 @@ int grib_encode_unsigned_long(unsigned char* p, unsigned long val, long* bitp, l
* Note: On x64 Micrsoft Windows a "long" is 32 bits but "size_t" is 64 bits * Note: On x64 Micrsoft Windows a "long" is 32 bits but "size_t" is 64 bits
*/ */
#define BIT_MASK_SIZE_T(x) \ #define BIT_MASK_SIZE_T(x) \
(((x) == max_nbits_size_t) ? (size_t)-1UL : (1UL << (x)) - 1) (((x) == max_nbits_size_t) ? (size_t)-1ULL : (1ULL << (x)) - 1)
size_t grib_decode_size_t(const unsigned char* p, long* bitp, long nbits) size_t grib_decode_size_t(const unsigned char* p, long* bitp, long nbits)
{ {

2
src/grib_bits_fast_big_endian.c
View File

@ -41,7 +41,7 @@ static void init_bits_all_one()
*(--v) = ~(cmask << --size); *(--v) = ~(cmask << --size);
} }
int grib_is_all_bits_one(long val, long nbits) int grib_is_all_bits_one(int64_t val, long nbits)
{ {
if (!bits_all_one.inited) if (!bits_all_one.inited)
init_bits_all_one(); init_bits_all_one();

384
src/grib_iterator_class_lambert_azimuthal_equal_area.c
View File

@ -8,11 +8,6 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/ */
/**************************************
* Enrico Fucile
**************************************/
#include "grib_api_internal.h" #include "grib_api_internal.h"
#include <math.h> #include <math.h>
@ -102,111 +97,235 @@ static int next(grib_iterator* i, double* lat, double* lon, double* val)
return 1; return 1;
} }
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args) #ifndef M_PI
#define M_PI 3.14159265358979323846 /* Whole pie */
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923 /* Half a pie */
#endif
#ifndef M_PI_4
#define M_PI_4 0.78539816339744830962 /* Quarter of a pie */
#endif
#define RAD2DEG 57.29577951308232087684 /* 180 over pi */
#define DEG2RAD 0.01745329251994329576 /* pi over 180 */
#define P00 .33333333333333333333 /* 1 / 3 */
#define P01 .17222222222222222222 /* 31 / 180 */
#define P02 .10257936507936507937 /* 517 / 5040 */
#define P10 .06388888888888888888 /* 23 / 360 */
#define P11 .06640211640211640212 /* 251 / 3780 */
#define P20 .01677689594356261023 /* 761 / 45360 */
static void pj_authset(double es, double* APA)
{
double t;
APA[0] = es * P00;
t = es * es;
APA[0] += t * P01;
APA[1] = t * P10;
t *= es;
APA[0] += t * P02;
APA[1] += t * P11;
APA[2] = t * P20;
}
static double pj_authlat(double beta, double* APA)
{
double t = beta + beta;
return (beta + APA[0] * sin(t) + APA[1] * sin(t + t) + APA[2] * sin(t + t + t));
}
static double pj_qsfn(double sinphi, double e, double one_es)
{
double con, div1, div2;
const double EPSILON = 1.0e-7;
if (e >= EPSILON) {
con = e * sinphi;
div1 = 1.0 - con * con;
div2 = 1.0 + con;
/* avoid zero division, fail gracefully */
if (div1 == 0.0 || div2 == 0.0)
return HUGE_VAL;
return (one_es * (sinphi / div1 - (.5 / e) * log((1. - con) / div2)));
}
else
return (sinphi + sinphi);
}
#define EPS10 1.e-10
static int init_oblate(grib_handle* h,
grib_iterator_lambert_azimuthal_equal_area* self,
size_t nv, long nx, long ny,
double Dx, double Dy, double earthMinorAxisInMetres, double earthMajorAxisInMetres,
double latFirstInRadians, double lonFirstInRadians,
double centralLongitudeInRadians, double standardParallelInRadians,
long iScansNegatively, long jScansPositively, long jPointsAreConsecutive)
{ {
int ret = 0;
double *lats, *lons; double *lats, *lons;
double lonFirstInDegrees, latFirstInDegrees, lonFirst, latFirst, radius = 0; long i, j;
long nx, ny, standardParallel, centralLongitude; double x0, y0, x, y;
double phi1, lambda0, xFirst, yFirst, x, y, Dx, Dy; double coslam, sinlam, sinphi, sinphi_, q, sinb = 0.0, cosb = 0.0, b = 0.0, cosb2;
double Q__qp = 0, Q__rq = 0, Q__cosb1, Q__sinb1, Q__dd, Q__xmf, Q__ymf, t;
/* double Q__mmf = 0; */
double e, es, temp, one_es;
double APA[3] = {0,};
double xFirst, yFirst;
Dx = iScansNegatively == 0 ? Dx / 1000 : -Dx / 1000;
Dy = jScansPositively == 1 ? Dy / 1000 : -Dy / 1000;
temp = (earthMajorAxisInMetres - earthMinorAxisInMetres) / earthMajorAxisInMetres;
es = 2 * temp - temp * temp;
one_es = 1.0 - es;
e = sqrt(es);
coslam = cos(lonFirstInRadians - centralLongitudeInRadians); /* cos(lp.lam) */
sinlam = sin(lonFirstInRadians - centralLongitudeInRadians);
sinphi = sin(latFirstInRadians); /* sin(lp.phi) */
q = pj_qsfn(sinphi, e, one_es);
t = fabs(standardParallelInRadians);
if (t > M_PI_2 + EPS10) {
return GRIB_GEOCALCULUS_PROBLEM;
}
/* if (fabs(t - M_HALFPI) < EPS10)
Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
else if (fabs(t) < EPS10)
Q->mode = EQUIT;
else
Q->mode = OBLIQ;
*/
Q__qp = pj_qsfn(1.0, e, one_es);
/* Q__mmf = 0.5 / one_es; ---- TODO(masn): do I need this? */
pj_authset(es, APA); /* sets up APA array */
Q__rq = sqrt(0.5 * Q__qp);
sinphi_ = sin(standardParallelInRadians); /* (P->phi0); */
Q__sinb1 = pj_qsfn(sinphi_, e, one_es) / Q__qp;
Q__cosb1 = sqrt(1.0 - Q__sinb1 * Q__sinb1);
Q__dd = cos(standardParallelInRadians) / (sqrt(1. - es * sinphi_ * sinphi_) * Q__rq * Q__cosb1);
Q__ymf = (Q__xmf = Q__rq) / Q__dd;
Q__xmf *= Q__dd;
sinb = q / Q__qp;
cosb2 = 1.0 - sinb * sinb;
cosb = cosb2 > 0 ? sqrt(cosb2) : 0;
b = 1. + Q__sinb1 * sinb + Q__cosb1 * cosb * coslam;
if (fabs(b) < EPS10) {
return GRIB_GEOCALCULUS_PROBLEM;
}
b = sqrt(2.0 / b);
/* OBLIQUE */
y0 = Q__ymf * b * (Q__cosb1 * sinb - Q__sinb1 * cosb * coslam);
x0 = Q__xmf * b * cosb * sinlam;
/* Allocate latitude and longitude arrays */
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", nv * sizeof(double));
return GRIB_OUT_OF_MEMORY;
}
lats = self->lats;
lons = self->lons;
/* Populate the lat and lon arrays */
{
xFirst = x0;
yFirst = y0;
y = yFirst;
for (j = 0; j < ny; j++) {
x = xFirst;
for (i = 0; i < nx; i++) {
double cCe, sCe, rho, ab = 0.0, lp__lam, lp__phi, xy_x = x, xy_y = y;
xy_x /= Q__dd;
xy_y *= Q__dd;
rho = hypot(xy_x, xy_y);
Assert(rho >= EPS10); /* TODO(masn): check */
sCe = 2. * asin(.5 * rho / Q__rq);
cCe = cos(sCe);
sCe = sin(sCe);
xy_x *= sCe;
/* if oblique */
ab = cCe * Q__sinb1 + xy_y * sCe * Q__cosb1 / rho;
xy_y = rho * Q__cosb1 * cCe - xy_y * Q__sinb1 * sCe;
/* else
ab = xy.y * sCe / rho;
xy.y = rho * cCe;
*/
lp__lam = atan2(xy_x, xy_y); /* longitude */
lp__phi = pj_authlat(asin(ab), APA); /* latitude */
*lats = lp__phi * RAD2DEG;
*lons = (lp__lam + centralLongitudeInRadians) * RAD2DEG;
lons++;
lats++;
x += Dx / earthMajorAxisInMetres;
}
y += Dy / earthMajorAxisInMetres;
}
}
return GRIB_SUCCESS;
}
static int init_sphere(grib_handle* h,
grib_iterator_lambert_azimuthal_equal_area* self,
size_t nv, long nx, long ny,
double Dx, double Dy, double radius,
double latFirstInRadians, double lonFirstInRadians,
double centralLongitudeInRadians, double standardParallelInRadians,
long iScansNegatively, long jScansPositively, long jPointsAreConsecutive)
{
double *lats, *lons;
double phi1, lambda0, xFirst, yFirst, x, y;
double kp, sinphi1, cosphi1; double kp, sinphi1, cosphi1;
long alternativeRowScanning, iScansNegatively;
long jScansPositively, jPointsAreConsecutive;
double sinphi, cosphi, cosdlambda, sindlambda; double sinphi, cosphi, cosdlambda, sindlambda;
double cosc, sinc; double cosc, sinc;
long i, j; long i, j;
grib_iterator_lambert_azimuthal_equal_area* self =
(grib_iterator_lambert_azimuthal_equal_area*)iter;
const char* sradius = grib_arguments_get_name(h, args, self->carg++);
const char* snx = grib_arguments_get_name(h, args, self->carg++);
const char* sny = grib_arguments_get_name(h, args, self->carg++);
const char* slatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* slonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sstandardParallel = grib_arguments_get_name(h, args, self->carg++);
const char* scentralLongitude = grib_arguments_get_name(h, args, self->carg++);
const char* sDx = grib_arguments_get_name(h, args, self->carg++);
const char* sDy = grib_arguments_get_name(h, args, self->carg++);
const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++);
const char* sjScansPositively = grib_arguments_get_name(h, args, self->carg++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
const char* salternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
double c, rho; double c, rho;
double epsilon = 1.0e-20; const double epsilon = 1.0e-20;
double d2r = acos(0.0) / 90.0; const double d2r = acos(0.0) / 90.0;
if ((ret = grib_get_double_internal(h, sradius, &radius)) != GRIB_SUCCESS) { lambda0 = centralLongitudeInRadians;
/* Check if it's an oblate spheroid */ phi1 = standardParallelInRadians;
if (grib_is_earth_oblate(h))
grib_context_log(h->context, GRIB_LOG_ERROR, "Lambert Azimuthal Equal Area only supported for spherical earth.");
return ret;
}
if ((ret = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS)
return ret;
if (iter->nv != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR,
"Wrong number of points (%ld!=%ldx%ld)",
iter->nv, nx, ny);
return GRIB_WRONG_GRID;
}
if ((ret = grib_get_double_internal(h, slatFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_double_internal(h, slonFirstInDegrees, &lonFirstInDegrees)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_long_internal(h, sstandardParallel, &standardParallel)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_long_internal(h, scentralLongitude, &centralLongitude)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_double_internal(h, sDx, &Dx)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_double_internal(h, sDy, &Dy)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_long_internal(h,
sjPointsAreConsecutive, &jPointsAreConsecutive)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_long_internal(h, sjScansPositively, &jScansPositively)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_long_internal(h, siScansNegatively, &iScansNegatively)) != GRIB_SUCCESS)
return ret;
if ((ret = grib_get_long_internal(h,
salternativeRowScanning, &alternativeRowScanning)) != GRIB_SUCCESS)
return ret;
lambda0 = d2r * centralLongitude / 1000000;
phi1 = d2r * standardParallel / 1000000;
latFirst = latFirstInDegrees * d2r;
lonFirst = lonFirstInDegrees * d2r;
cosphi1 = cos(phi1); cosphi1 = cos(phi1);
sinphi1 = sin(phi1); sinphi1 = sin(phi1);
Dx = iScansNegatively == 0 ? Dx / 1000 : -Dx / 1000; Dx = iScansNegatively == 0 ? Dx / 1000 : -Dx / 1000;
Dy = jScansPositively == 1 ? Dy / 1000 : -Dy / 1000; Dy = jScansPositively == 1 ? Dy / 1000 : -Dy / 1000;
self->lats = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double)); self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!self->lats) {
grib_context_log(h->context, GRIB_LOG_ERROR, grib_context_log(h->context, GRIB_LOG_ERROR,
"Error allocating %ld bytes", iter->nv * sizeof(double)); "Error allocating %ld bytes", nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
self->lons = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double)); self->lons = (double*)grib_context_malloc(h->context, nv * sizeof(double));
if (!self->lats) { if (!self->lats) {
grib_context_log(h->context, GRIB_LOG_ERROR, grib_context_log(h->context, GRIB_LOG_ERROR,
"Error allocating %ld bytes", iter->nv * sizeof(double)); "Error allocating %ld bytes", nv * sizeof(double));
return GRIB_OUT_OF_MEMORY; return GRIB_OUT_OF_MEMORY;
} }
lats = self->lats; lats = self->lats;
lons = self->lons; lons = self->lons;
/* compute xFirst,yFirst in metres */ /* compute xFirst,yFirst in metres */
sinphi = sin(latFirst); sinphi = sin(latFirstInRadians);
cosphi = cos(latFirst); cosphi = cos(latFirstInRadians);
cosdlambda = cos(lonFirst - lambda0); cosdlambda = cos(lonFirstInRadians - lambda0);
sindlambda = sin(lonFirst - lambda0); sindlambda = sin(lonFirstInRadians - lambda0);
kp = radius * sqrt(2.0 / (1 + sinphi1 * sinphi + cosphi1 * cosphi * cosdlambda)); kp = radius * sqrt(2.0 / (1 + sinphi1 * sinphi + cosphi1 * cosphi * cosdlambda));
xFirst = kp * cosphi * sindlambda; xFirst = kp * cosphi * sindlambda;
yFirst = kp * (cosphi1 * sinphi - sinphi1 * cosphi * cosdlambda); yFirst = kp * (cosphi1 * sinphi - sinphi1 * cosphi * cosdlambda);
@ -267,10 +386,105 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
y += Dy; y += Dy;
} }
} }
return GRIB_SUCCESS;
}
static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
{
int err = 0;
int is_oblate = 0;
double lonFirstInDegrees, latFirstInDegrees, lonFirstInRadians, latFirstInRadians, radius = 0;
long nx, ny;
double standardParallelInDegrees, centralLongitudeInDegrees;
double standardParallelInRadians, centralLongitudeInRadians;
double Dx, Dy;
long alternativeRowScanning, iScansNegatively;
long jScansPositively, jPointsAreConsecutive;
double earthMajorAxisInMetres = 0, earthMinorAxisInMetres = 0;
grib_iterator_lambert_azimuthal_equal_area* self =
(grib_iterator_lambert_azimuthal_equal_area*)iter;
const char* sradius = grib_arguments_get_name(h, args, self->carg++);
const char* snx = grib_arguments_get_name(h, args, self->carg++);
const char* sny = grib_arguments_get_name(h, args, self->carg++);
const char* slatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* slonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
const char* sstandardParallel = grib_arguments_get_name(h, args, self->carg++);
const char* scentralLongitude = grib_arguments_get_name(h, args, self->carg++);
const char* sDx = grib_arguments_get_name(h, args, self->carg++);
const char* sDy = grib_arguments_get_name(h, args, self->carg++);
const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++);
const char* sjScansPositively = grib_arguments_get_name(h, args, self->carg++);
const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
const char* salternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
const double d2r = acos(0.0) / 90.0;
is_oblate = grib_is_earth_oblate(h);
if (is_oblate) {
if ((err = grib_get_double_internal(h, "earthMinorAxisInMetres", &earthMinorAxisInMetres)) != GRIB_SUCCESS) return err;
if ((err = grib_get_double_internal(h, "earthMajorAxisInMetres", &earthMajorAxisInMetres)) != GRIB_SUCCESS) return err;
}
else {
if ((err = grib_get_double_internal(h, sradius, &radius)) != GRIB_SUCCESS) return err;
}
if ((err = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS)
return err;
if (iter->nv != nx * ny) {
grib_context_log(h->context, GRIB_LOG_ERROR,
"Wrong number of points (%ld!=%ldx%ld)",
iter->nv, nx, ny);
return GRIB_WRONG_GRID;
}
if ((err = grib_get_double_internal(h, slatFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(h, slonFirstInDegrees, &lonFirstInDegrees)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(h, sstandardParallel, &standardParallelInDegrees)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(h, scentralLongitude, &centralLongitudeInDegrees)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(h, sDx, &Dx)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(h, sDy, &Dy)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(h, sjPointsAreConsecutive, &jPointsAreConsecutive)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(h, sjScansPositively, &jScansPositively)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(h, siScansNegatively, &iScansNegatively)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(h, salternativeRowScanning, &alternativeRowScanning)) != GRIB_SUCCESS)
return err;
latFirstInRadians = latFirstInDegrees * d2r;
lonFirstInRadians = lonFirstInDegrees * d2r;
centralLongitudeInRadians = centralLongitudeInDegrees * d2r;
standardParallelInRadians = standardParallelInDegrees * d2r;
if (is_oblate) {
err = init_oblate(h, self, iter->nv, nx, ny,
Dx, Dy, earthMinorAxisInMetres, earthMajorAxisInMetres,
latFirstInRadians, lonFirstInRadians,
centralLongitudeInRadians, standardParallelInRadians,
iScansNegatively, jScansPositively, jPointsAreConsecutive);
}
else {
err = init_sphere(h, self, iter->nv, nx, ny,
Dx, Dy, radius,
latFirstInRadians, lonFirstInRadians,
centralLongitudeInRadians, standardParallelInRadians,
iScansNegatively, jScansPositively, jPointsAreConsecutive);
}
if (err) return err;
iter->e = -1; iter->e = -1;
return ret; return err;
} }
static int destroy(grib_iterator* i) static int destroy(grib_iterator* i)

22
src/grib_nearest_class_regular.c
View File

@ -232,7 +232,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
int ret = 0, kk = 0, ii = 0, jj = 0; int ret = 0, kk = 0, ii = 0, jj = 0;
size_t nvalues = 0; size_t nvalues = 0;
long iradius; long iradius;
double radius; double radius; /* radius in km */
grib_iterator* iter = NULL; grib_iterator* iter = NULL;
double lat = 0, lon = 0; double lat = 0, lon = 0;
@ -248,14 +248,24 @@ static int find(grib_nearest* nearest, grib_handle* h,
return ret; return ret;
nearest->values_count = nvalues; nearest->values_count = nvalues;
if (grib_is_missing(h, self->radius, &ret)) { /* 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_long(h, self->radius, &iradius)) == GRIB_SUCCESS) {
if (grib_is_missing(h, self->radius, &ret) || iradius == GRIB_MISSING_LONG) {
grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->radius); grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->radius);
return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; return GRIB_GEOCALCULUS_PROBLEM;
} }
if ((ret = grib_get_long(h, self->radius, &iradius)) != GRIB_SUCCESS)
return ret;
radius = ((double)iradius) / 1000.0; radius = ((double)iradius) / 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;
radius = (major + minor) / 2.0;
radius = radius / 1000.0;
}
/* Compute lat/lon info, create iterator etc if it's the 1st time or different grid. /* 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 * This is for performance: if the grid has not changed, we only do this once

42
tests/grib_grid_lamb_az_eq_area.sh
View File

@ -8,7 +8,7 @@
# virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. # virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
. ./include.sh . ./include.sh
#set -x
GRIB_INFILE=${data_dir}/regular_gaussian_pressure_level_constant.grib2 GRIB_INFILE=${data_dir}/regular_gaussian_pressure_level_constant.grib2
REF_FILE=grib_lamb_az_eq_area.ref REF_FILE=grib_lamb_az_eq_area.ref
@ -19,6 +19,9 @@ GRIB_OUTFILE=lamb_az_eq_area.grib2
DATA_OUTFILE=lamb_data.txt DATA_OUTFILE=lamb_data.txt
rm -f $FILTER_FILE $GRIB_OUTFILE $DATA_OUTFILE rm -f $FILTER_FILE $GRIB_OUTFILE $DATA_OUTFILE
# Spherical Earth
# ----------------
# Create a filter # Create a filter
cat > $FILTER_FILE<<EOF cat > $FILTER_FILE<<EOF
set edition = 2; set edition = 2;
@ -45,11 +48,8 @@ EOF
# Use this filter and the input GRIB to create a new GRIB # Use this filter and the input GRIB to create a new GRIB
rm -f "$GRIB_OUTFILE" rm -f "$GRIB_OUTFILE"
${tools_dir}/grib_filter -o $GRIB_OUTFILE $FILTER_FILE $GRIB_INFILE ${tools_dir}/grib_filter -o $GRIB_OUTFILE $FILTER_FILE $GRIB_INFILE
if [ ! -f "$GRIB_OUTFILE" ]; then
echo Failed to create output GRIB from filter >&2 # Now run the Geoiterator on the newly created GRIB file
exit 1
fi
# Now get the data from the newly created GRIB file
${tools_dir}/grib_get_data $GRIB_OUTFILE > $DATA_OUTFILE ${tools_dir}/grib_get_data $GRIB_OUTFILE > $DATA_OUTFILE
# Compare output with reference. If the diff fails, script will immediately exit with status 1 # Compare output with reference. If the diff fails, script will immediately exit with status 1
@ -62,6 +62,36 @@ grib_check_key_equals $GRIB_OUTFILE xDirectionGridLengthInMetres,yDirectionGridL
${tools_dir}/grib_ls -l 67,-33,1 $GRIB_OUTFILE ${tools_dir}/grib_ls -l 67,-33,1 $GRIB_OUTFILE
# Oblate spheroid
# --------------------
cat > $FILTER_FILE<<EOF
set edition = 2;
set gridType = "lambert_azimuthal_equal_area";
set Nx = 10;
set Ny = 10;
set values = {2};
set numberOfDataPoints = 100;
set shapeOfTheEarth = 4; # Earth assumed oblate spheroid
set numberOfValues = 100;
set latitudeOfFirstGridPointInDegrees = 67.575;
set longitudeOfFirstGridPointInDegrees = 326.5056;
set Dx = 5000000;
set Dy = 5000000;
set standardParallel = 48000000;
set centralLongitude = 9000000;
write;
EOF
# Use this filter and the input GRIB to create a new GRIB
rm -f "$GRIB_OUTFILE"
${tools_dir}/grib_filter -o $GRIB_OUTFILE $FILTER_FILE $GRIB_INFILE
${tools_dir}/grib_get_data $GRIB_OUTFILE
# Clean up # Clean up
rm -f $FILTER_FILE $DATA_OUTFILE rm -f $FILTER_FILE $DATA_OUTFILE
rm -f $GRIB_OUTFILE rm -f $GRIB_OUTFILE

8
tests/grib_nearest_test.sh
View File

@ -51,5 +51,13 @@ Idx lat lon dist
EOF EOF
diff $tempRef $temp diff $tempRef $temp
# ECC-1295: regular lat/lon on ellipsoid
# ----------------------------------------
sample2=$ECCODES_SAMPLES_PATH/GRIB2.tmpl
${tools_dir}/grib_set -s shapeOfTheEarth=5 $sample2 $temp
grib_check_key_equals $sample2 earthIsOblate 0
grib_check_key_equals $temp earthIsOblate 1
${tools_dir}/grib_ls -l 0,0 $temp
# Clean up # Clean up
rm -f $temp $tempRef rm -f $temp $tempRef