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

ECC-529: clean up

This commit is contained in:
Shahram Najm 2017-08-22 11:39:34 +01:00
parent 041396bd94
commit 21ca0ec6a0
1 changed files with 30 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

78
src/grib_iterator_class_polar_stereographic.c
View File

@ -96,29 +96,17 @@ static int next(grib_iterator* i, double *lat, double *lon, double *val)
return 1;
}
/* Data struct for Forward Projection i.e. from lat,lon to x,y */
typedef struct forward_proj_data {
double centre_lon; /* centre longitude */
double centre_lat; /* centre latitude */
double fac; /* sign variable */
double ind; /* flag variable */
double mcs; /* small m */
double tcs; /* small t */
double false_northing; /* y offset in meters */
double false_easting; /* x offset in meters */
} forward_proj_data;
/* Data struct for Inverse Projection i.e. from x,y to lat,lon */
typedef struct inverse_proj_data {
double centre_lon; /* centre longitude */
double centre_lat; /* centre latitude */
double fac; /* sign variable */
/* Data struct for Forward and Inverse Projections */
typedef struct proj_data_t {
double centre_lon; /* central longitude */
double centre_lat; /* central latitude */
double sign; /* sign variable */
double ind; /* flag variable */
double mcs; /* small m */
double tcs; /* small t */
double false_northing; /* y offset in meters */
double false_easting; /* x offset in meters */
} inverse_proj_data;
} proj_data_t;
#define RAD2DEG 57.29577951308232087684 /* 180 over pi */
#define DEG2RAD 0.01745329251994329576 /* pi over 180 */
@ -128,22 +116,19 @@ typedef struct inverse_proj_data {
static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
{
int ret=0;
double *lats,*lons;
double lonFirstInDegrees, latFirstInDegrees, radius=0;
long nx,ny,centralLongitudeInDegrees,centralLatitudeInDegrees;
double *lats, *lons; /* arrays for latitudes and longitudes */
double lonFirstInDegrees, latFirstInDegrees, radius;
double x, y, Dx, Dy;
long alternativeRowScanning, iScansNegatively;
long nx, ny, centralLongitudeInDegrees, centralLatitudeInDegrees;
long alternativeRowScanning, iScansNegatively, i, j;
long jScansPositively, jPointsAreConsecutive, southPoleOnPlane;
double cosphi;
double centralLongitude, centralLatitude; /* in radians */
double con1; /* temporary angle */
double con2; /* adjusted latitude */
double ts; /* value of small t */
double height; /* height above ellipsoid */
double x0, y0, lonFirst, latFirst;
forward_proj_data fwd_proj_data = {0,};
inverse_proj_data inv_proj_data = {0,};
long i,j;
proj_data_t fwd_proj_data = {0,};
proj_data_t inv_proj_data = {0,};
grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)iter;
@ -203,27 +188,26 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
fwd_proj_data.false_easting = 0;
fwd_proj_data.centre_lon = centralLongitude;
fwd_proj_data.centre_lat = centralLatitude;
if (centralLatitude < 0) fwd_proj_data.fac = -1.0;
else fwd_proj_data.fac = +1.0;
if (centralLatitude < 0) fwd_proj_data.sign = -1.0;
else fwd_proj_data.sign = +1.0;
fwd_proj_data.ind = 0;
if (fabs(fabs(centralLatitude) - PI_OVER_2) > EPSILON) {
/* central latitude different from 90 i.e. not north/south polar */
fwd_proj_data.ind = 1;
con1 = fwd_proj_data.fac * centralLatitude;
cosphi = cos(con1);
fwd_proj_data.mcs = cosphi;
con1 = fwd_proj_data.sign * centralLatitude;
fwd_proj_data.mcs = cos(con1);
fwd_proj_data.tcs = tan(0.5 * (PI_OVER_2 - con1));
}
/* Forward projection from initial lat,lon to initial x,y */
con1 = fwd_proj_data.fac * (lonFirst - fwd_proj_data.centre_lon);
con2 = fwd_proj_data.fac * latFirst;
ts = tan(0.5 * (PI_OVER_2 - con2));
con1 = fwd_proj_data.sign * (lonFirst - fwd_proj_data.centre_lon);
ts = tan(0.5 * (PI_OVER_2 - fwd_proj_data.sign * latFirst));
if (fwd_proj_data.ind)
height = radius * fwd_proj_data.mcs * ts / fwd_proj_data.tcs;
else
height = 2.0 * radius * ts;
x0 = fwd_proj_data.fac * height * sin(con1) + fwd_proj_data.false_easting;
y0 = -fwd_proj_data.fac * height * cos(con1) + fwd_proj_data.false_northing;
x0 = fwd_proj_data.sign * height * sin(con1) + fwd_proj_data.false_easting;
y0 = -fwd_proj_data.sign * height * cos(con1) + fwd_proj_data.false_northing;
x0 = -x0;
y0 = -y0;
@ -233,14 +217,13 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
inv_proj_data.false_northing= y0;
inv_proj_data.centre_lon = centralLongitude;
inv_proj_data.centre_lat = centralLatitude;
if (centralLatitude < 0) inv_proj_data.fac = -1.0;
else inv_proj_data.fac = +1.0;
if (centralLatitude < 0) inv_proj_data.sign = -1.0;
else inv_proj_data.sign = +1.0;
inv_proj_data.ind = 0;
if (fabs(fabs(centralLatitude) - PI_OVER_2) > EPSILON) {
inv_proj_data.ind = 1;
con1 = inv_proj_data.fac * inv_proj_data.centre_lat;
cosphi = cos(con1);
inv_proj_data.mcs = cosphi;
con1 = inv_proj_data.sign * inv_proj_data.centre_lat;
inv_proj_data.mcs = cos(con1);
inv_proj_data.tcs = tan(0.5 * (PI_OVER_2 - con1));
}
self->lats = (double*)grib_context_malloc(h->context,iter->nv*sizeof(double));
@ -264,25 +247,24 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
for (i=0;i<nx;i++) {
/* Inverse projection from x,y to lat,lon */
/* int index =i+j*nx; */
double _x = (x - inv_proj_data.false_easting) * inv_proj_data.fac;
double _y = (y - inv_proj_data.false_northing) * inv_proj_data.fac;
double _x = (x - inv_proj_data.false_easting) * inv_proj_data.sign;
double _y = (y - inv_proj_data.false_northing) * inv_proj_data.sign;
double rh = sqrt(_x * _x + _y * _y);
if (inv_proj_data.ind)
ts = rh * inv_proj_data.tcs/(radius * inv_proj_data.mcs);
else
ts = rh / (radius * 2.0);
*lats = inv_proj_data.fac * (PI_OVER_2 - 2 * atan(ts));
*lats = inv_proj_data.sign * (PI_OVER_2 - 2 * atan(ts));
if (rh == 0) {
*lons = inv_proj_data.fac * inv_proj_data.centre_lon;
*lons = inv_proj_data.sign * inv_proj_data.centre_lon;
} else {
double temp = atan2(_x, -_y);
*lons = inv_proj_data.fac * temp + inv_proj_data.centre_lon;
*lons = inv_proj_data.sign * temp + inv_proj_data.centre_lon;
}
*lats = *lats * RAD2DEG;
*lons = *lons * RAD2DEG;
while (*lons<0) *lons += 360;
while (*lons>360) *lons -= 360;
/* printf("DBK: llat[%d] = %g \t llon[%d] = %g\n", index,*lats, index,*lons); */
lons++;
lats++;