mirror of https://github.com/ecmwf/eccodes.git
ECC-383: Implement iterator for 'space view'. Alternative impl based on MSGnavigation
This commit is contained in:
Shahram Najm
parent
ac4e34d9d3
commit
2bd72b8dad
|
|
@ -99,6 +99,260 @@ 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 */
|
||||
|
||||
#if 0
|
||||
static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
|
||||
{
|
||||
int ret = GRIB_SUCCESS;
|
||||
double *lats,*lons;
|
||||
double latOfSubSatellitePointInDegrees,lonOfSubSatellitePointInDegrees;
|
||||
double orientationInDegrees, nrInRadiusOfEarth;
|
||||
double radius=0,xpInGridLengths=0,ypInGridLengths=0;
|
||||
long nx, ny, earthIsOblate;
|
||||
long alternativeRowScanning,iScansNegatively;
|
||||
long Xo, Yo;
|
||||
long jScansPositively,jPointsAreConsecutive;
|
||||
long i;
|
||||
|
||||
double major, minor, r_eq, r_pol, sat_height;
|
||||
double lap, lop, orient_angle, angular_size;
|
||||
double xp, yp, dx, dy, rx, ry, x, y;
|
||||
double cos_x, cos_y, sin_x, sin_y;
|
||||
double factor_1, factor_2, tmp1, Sd, Sn, Sxy, S1, S2, S3;
|
||||
int x0, y0, ix, iy;
|
||||
double *s_x, *c_x;
|
||||
double const1, const2, const3;
|
||||
|
||||
grib_iterator_space_view* self = (grib_iterator_space_view*)iter;
|
||||
|
||||
const char* sradius = grib_arguments_get_name(h,args,self->carg++);
|
||||
const char* sEarthIsOblate = grib_arguments_get_name(h,args,self->carg++);
|
||||
const char* sMajorAxisInMetres = grib_arguments_get_name(h,args,self->carg++);
|
||||
const char* sMinorAxisInMetres = 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* sLatOfSubSatellitePointInDegrees = grib_arguments_get_name(h,args,self->carg++);
|
||||
const char* sLonOfSubSatellitePointInDegrees = 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* sXpInGridLengths = grib_arguments_get_name(h,args,self->carg++);
|
||||
const char* sYpInGridLengths = grib_arguments_get_name(h,args,self->carg++);
|
||||
const char* sOrientationInDegrees = grib_arguments_get_name(h,args,self->carg++);
|
||||
const char* sNrInRadiusOfEarth = grib_arguments_get_name(h,args,self->carg++);
|
||||
const char* sXo = grib_arguments_get_name(h,args,self->carg++);
|
||||
const char* sYo = 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++);
|
||||
|
||||
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((ret = grib_get_long_internal(h, sEarthIsOblate, &earthIsOblate)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
if (earthIsOblate) {
|
||||
if((ret = grib_get_double_internal(h, sMajorAxisInMetres, &major)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(h, sMinorAxisInMetres, &minor)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
} else {
|
||||
if((ret = grib_get_double_internal(h, sradius, &radius)) != 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, sLatOfSubSatellitePointInDegrees,&latOfSubSatellitePointInDegrees)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(h, sLonOfSubSatellitePointInDegrees,&lonOfSubSatellitePointInDegrees)) != 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_double_internal(h, sXpInGridLengths,&xpInGridLengths)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(h, sYpInGridLengths,&ypInGridLengths)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(h, sOrientationInDegrees,&orientationInDegrees)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
/* Orthographic not supported. This happens when Nr (camera altitude) is missing */
|
||||
if (grib_is_missing(h, sNrInRadiusOfEarth, &ret)) {
|
||||
grib_context_log(h->context,GRIB_LOG_ERROR, "Orthographic view (Nr missing) not supported");
|
||||
return GRIB_NOT_IMPLEMENTED;
|
||||
}
|
||||
if((ret = grib_get_double_internal(h, sNrInRadiusOfEarth,&nrInRadiusOfEarth)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
if((ret = grib_get_long_internal(h, sXo,&Xo)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_long_internal(h, sYo,&Yo)) != 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;
|
||||
|
||||
if (earthIsOblate) {
|
||||
r_eq = major; /* In km */
|
||||
r_pol = minor;
|
||||
} else {
|
||||
r_eq = r_pol = radius * 0.001;/*conv to km*/
|
||||
}
|
||||
|
||||
double as = 2 * asin(1.0 / nrInRadiusOfEarth); /* angular_size */
|
||||
double cfac = dx / as;
|
||||
double lfac = dy / as;
|
||||
// use km, so scale by the earth radius
|
||||
double scale_factor = (Nr - 1) * majorAxis / 1000; // this sets the units of the projection x,y coords in km
|
||||
double scale_x = scale_factor; // LOOK fake neg need scan value
|
||||
double scale_y = -scale_factor; // LOOK fake neg need scan value
|
||||
|
||||
double startx = scale_factor * (1 - Xp) / cfac;
|
||||
double incrx = scale_factor / cfac;
|
||||
boolean yscanPositive = GribUtils.scanModeYisPositive(scanMode);
|
||||
double starty, incry;
|
||||
if (yscanPositive) {
|
||||
starty = scale_factor * (Yp - getNy()) / lfac;
|
||||
incry = (scale_factor / lfac);
|
||||
} else {
|
||||
incry = -(scale_factor / lfac);
|
||||
starty = scale_factor * (Yp - getNy()) / lfac - incry * (getNy() - 1);
|
||||
}
|
||||
|
||||
//MSGnavigation proj = new MSGnavigation(LaP, LoP, majorAxis, minorAxis, Nr * majorAxis, scale_x, scale_y);
|
||||
sat_height = nrInRadiusOfEarth * r_eq; /* km */
|
||||
const1 = major / minor;
|
||||
const1 *= const1;
|
||||
const2 = 1.0 - (minor * minor) / (major * major);
|
||||
const3 = sat_height * sat_height - major * major;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
///////OLD
|
||||
angular_size = 2.0 * asin(1.0 / nrInRadiusOfEarth);
|
||||
sat_height = nrInRadiusOfEarth * r_eq;
|
||||
|
||||
lap = latOfSubSatellitePointInDegrees;
|
||||
lop = lonOfSubSatellitePointInDegrees;
|
||||
/* apply default scaling factor */
|
||||
lap *= 1e-6;
|
||||
lop *= 1e-6;
|
||||
lap *= DEG2RAD;
|
||||
lop *= DEG2RAD;
|
||||
|
||||
orient_angle = orientationInDegrees;
|
||||
/* apply default scaling factor */
|
||||
if (orient_angle != 0.0) return GRIB_NOT_IMPLEMENTED;
|
||||
|
||||
xp = xpInGridLengths;
|
||||
yp = ypInGridLengths;
|
||||
x0 = Xo;
|
||||
y0 = Yo;
|
||||
|
||||
rx = angular_size / dx;
|
||||
ry = (r_pol/r_eq) * angular_size / dy;
|
||||
|
||||
/*dx = iScansNegatively == 0 ? dx : -dx;
|
||||
* dy = jScansPositively == 1 ? dy : -dy;*/
|
||||
self->lats = (double*)grib_context_malloc(h->context,iter->nv*sizeof(double));
|
||||
if (!self->lats) {
|
||||
grib_context_log(h->context,GRIB_LOG_ERROR, "unable to allocate %ld bytes",iter->nv*sizeof(double));
|
||||
return GRIB_OUT_OF_MEMORY;
|
||||
}
|
||||
self->lons = (double*)grib_context_malloc(h->context,iter->nv*sizeof(double));
|
||||
if (!self->lats) {
|
||||
grib_context_log(h->context,GRIB_LOG_ERROR, "unable to allocate %ld bytes",iter->nv*sizeof(double));
|
||||
return GRIB_OUT_OF_MEMORY;
|
||||
}
|
||||
lats=self->lats;
|
||||
lons=self->lons;
|
||||
|
||||
if (!iScansNegatively) {
|
||||
xp = xp - x0;
|
||||
} else {
|
||||
xp = (nx-1) - (xp - x0);
|
||||
}
|
||||
if (jScansPositively) {
|
||||
yp = yp - y0;
|
||||
}
|
||||
else {
|
||||
yp = (ny-1) - (yp - y0);
|
||||
}
|
||||
i = 0;
|
||||
factor_2 = (r_eq/r_pol)*(r_eq/r_pol);
|
||||
factor_1 = sat_height * sat_height - r_eq * r_eq;
|
||||
|
||||
s_x = (double *) grib_context_malloc(h->context, nx*sizeof(double));
|
||||
if (!s_x) {
|
||||
grib_context_log(h->context,GRIB_LOG_ERROR, "unable to allocate %ld bytes",nx*sizeof(double));
|
||||
return GRIB_OUT_OF_MEMORY;
|
||||
}
|
||||
c_x = (double *) grib_context_malloc(h->context, nx*sizeof(double));
|
||||
if (!c_x) {
|
||||
grib_context_log(h->context,GRIB_LOG_ERROR, "unable to allocate %ld bytes",nx*sizeof(double));
|
||||
return GRIB_OUT_OF_MEMORY;
|
||||
}
|
||||
|
||||
for (ix = 0; ix < nx; ix++) {
|
||||
x = (ix - xp) * rx;
|
||||
s_x[ix] = sin(x);
|
||||
c_x[ix] = sqrt(1.0 - s_x[ix]*s_x[ix]);
|
||||
}
|
||||
|
||||
for (iy = 0; iy < ny; iy++) {
|
||||
y = (iy - yp) * ry;
|
||||
sin_y = sin(y);
|
||||
cos_y = sqrt(1.0 - sin_y*sin_y);
|
||||
|
||||
tmp1 = (1 + (factor_2-1.0)*sin_y*sin_y);
|
||||
|
||||
for (ix = 0; ix < nx; ix++, i++) {
|
||||
/* x = (ix - xp) * rx; */
|
||||
sin_x = s_x[ix];
|
||||
cos_x = c_x[ix];
|
||||
|
||||
Sd = sat_height * cos_x * cos_y;
|
||||
Sd = Sd * Sd - tmp1*factor_1;
|
||||
if (Sd <= 0.0) { // outside of view
|
||||
lats[i] = lons[i] = 0; /* TODO: error? */
|
||||
}
|
||||
else {
|
||||
Sd = sqrt(Sd);
|
||||
Sn = (sat_height*cos_x*cos_y - Sd) / tmp1;
|
||||
S1 = sat_height - Sn * cos_x * cos_y;
|
||||
S2 = Sn * sin_x * cos_y;
|
||||
S3 = Sn * sin_y;
|
||||
Sxy = sqrt(S1*S1 + S2*S2);
|
||||
lons[i] = atan(S2/S1)*(RAD2DEG) + lop;
|
||||
lats[i] = atan(factor_2*S3/Sxy)*(RAD2DEG);
|
||||
/*printf("lat=%g lon=%g\n", lats[i], lons[i]);*/
|
||||
}
|
||||
/*while (lons[i]<0) lons[i] += 360;
|
||||
* while (lons[i]>360) lons[i] -= 360;*/
|
||||
}
|
||||
}
|
||||
grib_context_free(h->context, s_x);
|
||||
grib_context_free(h->context, c_x);
|
||||
iter->e = -1;
|
||||
|
||||
return ret;
|
||||
}
|
||||
#endif
|
||||
#if 1
|
||||
static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
|
||||
{
|
||||
int ret=0;
|
||||
|
|
@ -142,7 +396,7 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
|
|||
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 char* sAlternativeRowScanning = grib_arguments_get_name(h,args,self->carg++);
|
||||
|
||||
if((ret = grib_get_long_internal(h, snx,&nx)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
|
@ -179,8 +433,15 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
|
|||
return ret;
|
||||
if((ret = grib_get_double_internal(h, sOrientationInDegrees,&orientationInDegrees)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
/* Orthographic not supported. This happens when Nr (camera altitude) is missing */
|
||||
if (grib_is_missing(h, sNrInRadiusOfEarth, &ret)) {
|
||||
grib_context_log(h->context,GRIB_LOG_ERROR, "Orthographic view (Nr missing) not supported");
|
||||
return GRIB_NOT_IMPLEMENTED;
|
||||
}
|
||||
if((ret = grib_get_double_internal(h, sNrInRadiusOfEarth,&nrInRadiusOfEarth)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
if((ret = grib_get_long_internal(h, sXo,&Xo)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_long_internal(h, sYo,&Yo)) != GRIB_SUCCESS)
|
||||
|
|
@ -191,15 +452,14 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
|
|||
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)
|
||||
if((ret = grib_get_long_internal(h, sAlternativeRowScanning,&alternativeRowScanning)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
if (earthIsOblate) {
|
||||
r_eq = major ; /* In m */
|
||||
r_eq = major; /* In km */
|
||||
r_pol = minor;
|
||||
} else {
|
||||
return GRIB_NOT_IMPLEMENTED;
|
||||
r_eq = r_pol = radius * 0.001;
|
||||
r_eq = r_pol = radius * 0.001;/*conv to km*/
|
||||
}
|
||||
angular_size = 2.0 * asin(1.0 / nrInRadiusOfEarth);
|
||||
sat_height = nrInRadiusOfEarth * r_eq;
|
||||
|
|
@ -309,6 +569,7 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
|
|||
|
||||
return ret;
|
||||
}
|
||||
#endif
|
||||
|
||||
static int destroy(grib_iterator* i)
|
||||
{
|
||||
|
|
|
|||
Loading…
Reference in New Issue