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

Merge pull request #195 from ecmwf/bugfix/ECC-1364-LambertSouthern 

Bugfix/ECC-1364 lambert southern
This commit is contained in:
shahramn 2024-02-14 13:22:53 +00:00 committed by GitHub
parent 951ff99b50 404778c7c0
commit 08d09ead58
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
3 changed files with 109 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
src/grib_iterator_class_lambert_conformal.cc
View File

@ -162,7 +162,32 @@ static double calculate_eccentricity(double minor, double major)
return sqrt(1.0 - temp * temp);
}
static int init_sphere(grib_handle* h,
static void xy2latlon(double radius, double n, double f, double rho0_bare, double LoVInRadians,
double x, double y,
double* latDeg, double* lonDeg)
{
x /= radius;
y /= radius;
y = rho0_bare - y;
double rho = hypot(x, y);
if (rho != 0.0) {
if (n < 0.0) {
rho = -rho;
x = -x;
y = -y;
}
double lp_phi = 2. * atan(pow(f / rho, 1.0/n)) - M_PI_2;
double lp_lam = atan2(x, y) / n;
*lonDeg = lp_lam*RAD2DEG + LoVInRadians*RAD2DEG;
*latDeg = lp_phi*RAD2DEG;
}
else {
*lonDeg = 0.0;
*latDeg = (n > 0.0 ? M_PI_2 : -M_PI_2) * RAD2DEG;
}
}
static int init_sphere(const grib_handle* h,
grib_iterator_lambert_conformal* self,
size_t nv, long nx, long ny,
double LoVInDegrees,
@ -171,9 +196,7 @@ static int init_sphere(grib_handle* h,
double LoVInRadians, double Latin1InRadians, double Latin2InRadians,
double LaDInRadians)
{
int i, j;
double f, n, rho, rho0, angle, x0, y0, x, y, tmp, tmp2;
double latDeg, lonDeg, lonDiff;
double n, angle, x0, y0, x, y;
if (fabs(Latin1InRadians - Latin2InRadians) < 1E-09) {
n = sin(Latin1InRadians);
@ -182,12 +205,14 @@ static int init_sphere(grib_handle* h,
log(tan(M_PI_4 + Latin2InRadians / 2.0) / tan(M_PI_4 + Latin1InRadians / 2.0));
}
f = (cos(Latin1InRadians) * pow(tan(M_PI_4 + Latin1InRadians / 2.0), n)) / n;
rho = radius * f * pow(tan(M_PI_4 + latFirstInRadians / 2.0), -n);
rho0 = radius * f * pow(tan(M_PI_4 + LaDInRadians / 2.0), -n);
if (n < 0) /* adjustment for southern hemisphere */
rho0 = -rho0;
lonDiff = lonFirstInRadians - LoVInRadians;
double f = (cos(Latin1InRadians) * pow(tan(M_PI_4 + Latin1InRadians / 2.0), n)) / n;
double rho = radius * f * pow(tan(M_PI_4 + latFirstInRadians / 2.0), -n);
double rho0_bare = f * pow(tan(M_PI_4 + LaDInRadians / 2.0), -n);
double rho0 = radius * rho0_bare;
//if (n < 0) /* adjustment for southern hemisphere */
// rho0 = -rho0;
double lonDiff = lonFirstInRadians - LoVInRadians;
/* Adjust longitude to range -180 to 180 */
if (lonDiff > M_PI)
@ -213,20 +238,33 @@ static int init_sphere(grib_handle* h,
return GRIB_OUT_OF_MEMORY;
}
double latDeg = 0, lonDeg = 0;
for (long j = 0; j < ny; j++) {
y = y0 + j * Dy;
for (long i = 0; i < nx; i++) {
const long index = i + j * nx;
x = x0 + i * Dx;
xy2latlon(radius, n, f, rho0_bare, LoVInRadians, x, y, &latDeg, &lonDeg);
self->lons[index] = lonDeg;
self->lats[index] = latDeg;
}
}
#if 0
/* Populate our arrays */
for (j = 0; j < ny; j++) {
y = y0 + j * Dy;
if (n < 0) { /* adjustment for southern hemisphere */
y = -y;
}
//if (n < 0) { /* adjustment for southern hemisphere */
// y = -y;
//}
tmp = rho0 - y;
tmp2 = tmp * tmp;
for (i = 0; i < nx; i++) {
int index = i + j * nx;
x = x0 + i * Dx;
if (n < 0) { /* adjustment for southern hemisphere */
x = -x;
}
//printf("j=%d i=%d xy= %.6f %.6f\t",j,i,x,y);
//if (n < 0) { /* adjustment for southern hemisphere */
// x = -x;
//}
angle = atan2(x, tmp); /* See ECC-524 */
rho = sqrt(x * x + tmp2);
if (n <= 0) rho = -rho;
@ -237,12 +275,12 @@ static int init_sphere(grib_handle* h,
self->lats[index] = latDeg;
}
}
#endif
return GRIB_SUCCESS;
}
/* Oblate spheroid */
static int init_oblate(grib_handle* h,
static int init_oblate(const grib_handle* h,
grib_iterator_lambert_conformal* self,
size_t nv, long nx, long ny,
double LoVInDegrees,

1
tests/CMakeLists.txt
View File

@ -254,6 +254,7 @@ if( HAVE_BUILD_TOOLS )
grib_ecc-1000
grib_ecc-1001
grib_ecc-1030
grib_ecc-1364
grib_ecc-1397
grib_ecc-1425
grib_ecc-1467

52
tests/grib_ecc-1364.sh Executable file
View File

@ -0,0 +1,52 @@
#!/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.ctest.sh
# ECC-1364: GRIB: Geoiterator for Lambert Conformal in the southern hemisphere
label="grib_ecc-1364_test"
tempGrib=temp.$label.grib
tempFilt=temp.$label.filt
tempLog=temp.$label.log
tempRef=temp.$label.ref
sample_grib1=$ECCODES_SAMPLES_PATH/GRIB1.tmpl
# Create a GRIB with a similar grid to the one in the JIRA issue
cat >$tempFilt<<EOF
set dataRepresentationType = 3; # lambert conformal conic
set resolutionAndComponentFlags = 136;
set Nx = 999;
set Ny = 1249;
set latitudeOfFirstGridPoint = -54387;
set longitudeOfFirstGridPoint = 265669;
set LoV = 295000;
set DxInMetres = 4000;
set DyInMetres = 4000;
set projectionCentreFlag = 128;
set Latin1 = -35000;
set Latin2 = -35000;
set latitudeOfSouthernPole = -90000;
set scanningMode = 64; # +i +j
write;
EOF
${tools_dir}/grib_filter -o $tempGrib $tempFilt $sample_grib1
${tools_dir}/grib_get_data $tempGrib > $tempLog
${tools_dir}/grib_ls -l -11.6277,-47.9583,1 $tempGrib > $tempLog
grep -q "Grid Point chosen #1 index=1247750 " $tempLog
grep -q "index=1247750 .* distance=0.01 " $tempLog
# Clean up
rm -f $tempGrib $tempFilt $tempLog $tempRef