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

ECC-1055: Shape of the earth not taken into account for Lambert Conformal Conic (part 5)

This commit is contained in:
Shahram Najm 2020-03-28 17:52:19 +00:00
parent 0f898f4336
commit 6cb27aaf37
1 changed files with 57 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

99
src/grib_iterator_class_lambert_conformal.c
View File

@ -82,20 +82,7 @@ static void init_class(grib_iterator_class* c)
}
/* END_CLASS_IMP */
static int next(grib_iterator* i, double* lat, double* lon, double* val)
{
grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)i;
if ((long)i->e >= (long)(i->nv - 1))
return 0;
i->e++;
*lat = self->lats[i->e];
*lon = self->lons[i->e];
*val = i->data[i->e];
return 1;
}
#define EPSILON 1.0e-10
#ifndef M_PI
#define M_PI 3.14159265358979323846 /* Whole pie */
@ -112,12 +99,12 @@ static int next(grib_iterator* i, double* lat, double* lon, double* val)
#define RAD2DEG 57.29577951308232087684 /* 180 over pi */
#define DEG2RAD 0.01745329251994329576 /* pi over 180 */
/* Adjust longitude to range -180 to 180 */
static double adjust_lon(double x)
/* Adjust longitude (in radians) to range -180 to 180 */
static double adjust_lon_radians(double lon)
{
if (x > M_PI) x -= 2 * M_PI;
if (x < -M_PI) x += 2 * M_PI;
return x;
if (lon > M_PI) lon -= 2 * M_PI;
if (lon < -M_PI) lon += 2 * M_PI;
return lon;
}
/* Function to compute the latitude angle, phi2, for the inverse */
@ -130,14 +117,14 @@ double phi2z(
int i, MAX_ITER = 15;
eccnth = 0.5 * eccent;
phi = M_PI_2 - 2 * atan(ts);
phi = M_PI_2 - 2 * atan(ts);
for (i = 0; i <= MAX_ITER; i++) {
sinpi = sin(phi);
con = eccent * sinpi;
dphi = M_PI_2 - 2 * atan(ts *(pow(((1.0 - con)/(1.0 + con)),eccnth))) - phi;
con = eccent * sinpi;
dphi = M_PI_2 - 2 * atan(ts * (pow(((1.0 - con) / (1.0 + con)), eccnth))) - phi;
phi += dphi;
if (fabs(dphi) <= .0000000001)
return(phi);
return (phi);
}
/*Assert(!"Convergence error");*/
*error = GRIB_INTERNAL_ERROR;
@ -245,10 +232,8 @@ static int init_sphere(grib_handle* h,
rho = -rho;
lonDeg = LoVInDegrees + (angle / n) * RAD2DEG;
latDeg = (2.0 * atan(pow(radius * f / rho, 1.0 / n)) - M_PI_2) * RAD2DEG;
while (lonDeg >= 360.0)
lonDeg -= 360.0;
while (lonDeg < 0.0)
lonDeg += 360.0;
while (lonDeg >= 360.0) lonDeg -= 360.0;
while (lonDeg < 0.0) lonDeg += 360.0;
self->lons[index] = lonDeg;
self->lats[index] = latDeg;
/*printf("DBK: llat[%d] = %g \t llon[%d] = %g\n", index,lats[index], index,lons[index]);*/
@ -258,8 +243,6 @@ static int init_sphere(grib_handle* h,
return GRIB_SUCCESS;
}
#define EPSILON 1.0e-10
/* Oblate spheroid */
static int init_oblate(grib_handle* h,
grib_iterator_lambert_conformal* self,
@ -292,7 +275,8 @@ static int init_oblate(grib_handle* h,
cos_po = cos(Latin1InRadians);
con = sin_po;
ms1 = msfnz(e, sin_po, cos_po);
ts1 = tsfnz(e, Latin2InRadians, sin_po);
ts1 = tsfnz(e, Latin1InRadians, sin_po);
sin_po = sin(Latin2InRadians);
cos_po = cos(Latin2InRadians);
ms2 = msfnz(e, sin_po, cos_po);
@ -302,7 +286,8 @@ static int init_oblate(grib_handle* h,
if (fabs(Latin1InRadians - Latin2InRadians) > EPSILON) {
ns = log(ms1 / ms2) / log(ts1 / ts2);
} else {
}
else {
ns = con;
}
f0 = ms1 / (ns * pow(ts1, ns));
@ -314,7 +299,8 @@ static int init_oblate(grib_handle* h,
sinphi = sin(latFirstInRadians);
ts = tsfnz(e, latFirstInRadians, sinphi);
rh1 = earthMajorAxisInMetres * f0 * pow(ts, ns);
} else {
}
else {
con = latFirstInRadians * ns;
if (con <= 0) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Point cannot be projected");
@ -322,9 +308,11 @@ static int init_oblate(grib_handle* h,
}
rh1 = 0;
}
theta = ns * adjust_lon(lonFirstInRadians - LoVInRadians);
theta = ns * adjust_lon_radians(lonFirstInRadians - LoVInRadians);
x0 = rh1 * sin(theta);
y0 = rh - rh1 * cos(theta);
x0 = -x0;
y0 = -y0;
/* Allocate latitude and longitude arrays */
self->lats = (double*)grib_context_malloc(h->context, nv * sizeof(double));
@ -343,13 +331,14 @@ static int init_oblate(grib_handle* h,
y = y0 + j * Dy;
for (i = 0; i < nx; i++) {
const int index = i + j * nx;
x = x0 + i * Dx;
x = x0 + i * Dx;
/* from x,y to lat,lon */
y = rh - y;
if (ns > 0) {
rh1 = sqrt(x * x + y * y);
con = 1.0;
} else {
}
else {
rh1 = -sqrt(x * x + y * y);
con = -1.0;
}
@ -357,22 +346,26 @@ static int init_oblate(grib_handle* h,
if (rh1 != 0)
theta = atan2((con * x), (con * y));
if ((rh1 != 0) || (ns > 0.0)) {
con = 1.0 / ns;
ts = pow((rh1 / (earthMajorAxisInMetres * f0)), con);
con = 1.0 / ns;
ts = pow((rh1 / (earthMajorAxisInMetres * f0)), con);
latDeg = phi2z(e, ts, &err);
if (err) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Failed to compute the latitude angle, phi2, for the inverse");
return err;
}
} else {
}
else {
latDeg = -M_PI_2;
}
lonDeg = adjust_lon(theta / ns + LoVInRadians);
lonDeg = adjust_lon_radians(theta / ns + LoVInRadians);
// Convert to degrees
latDeg *= RAD2DEG;
lonDeg *= RAD2DEG;
while (lonDeg >= 360.0) lonDeg -= 360.0;
while (lonDeg < 0.0) lonDeg += 360.0;
self->lons[index] = lonDeg;
self->lats[index] = latDeg;
/*printf("DBK: llat[%d] = %g \t llon[%d] = %g\n", index,lats[index], index,lons[index]);*/
//printf("DBK: llat[%d] = %g \t llon[%d] = %g\n", index, self->lats[index], index, self->lons[index]);
}
}
return GRIB_SUCCESS;
@ -417,14 +410,15 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
if (is_oblate) {
if ((err = grib_get_double_internal(h, "earthMajorAxisInMetres", &earthMajorAxisInMetres)) != GRIB_SUCCESS) return err;
if ((err = grib_get_double_internal(h, "earthMinorAxisInMetres", &earthMinorAxisInMetres)) != GRIB_SUCCESS) return err;
} else {
if ((err = grib_get_double_internal(h, sradius, &radius)) != 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, sradius, &radius)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(h, sLoVInDegrees, &LoVInDegrees)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(h, sLaDInDegrees, &LaDInDegrees)) != GRIB_SUCCESS)
@ -450,6 +444,12 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
if ((err = grib_get_long_internal(h, salternativeRowScanning, &alternativeRowScanning)) != GRIB_SUCCESS)
return err;
/* Standard Parallels cannot be equal and on opposite sides of the equator */
if (fabs(Latin1InDegrees + Latin2InDegrees) < EPSILON) {
grib_context_log(h->context, GRIB_LOG_ERROR, "Cannot have equal latitudes for standard parallels on opposite sides of equator");
return GRIB_WRONG_GRID;
}
/*
* See Wolfram MathWorld: http://mathworld.wolfram.com/LambertConformalConicProjection.html
*/
@ -487,6 +487,21 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
return err;
}
static int next(grib_iterator* i, double* lat, double* lon, double* val)
{
grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)i;
if ((long)i->e >= (long)(i->nv - 1))
return 0;
i->e++;
*lat = self->lats[i->e];
*lon = self->lons[i->e];
*val = i->data[i->e];
return 1;
}
static int destroy(grib_iterator* i)
{
grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)i;