Skip to content

GitLab

Commit 554fd4c9 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added gridDefParamSTERE() and gridInqParamSTERE().

parent 24fba901
Hide whitespace changes
Inline Side-by-side
Showing with 123 additions and 80 deletions
src/cdi.h
View file @ 554fd4c9
......@@ -959,10 +959,14 @@ void gridInqParamRLL(int gridID, double *xpole, double *ypole, double *angle)
void gridDefParamGME(int gridID, int nd, int ni, int ni2, int ni3);
void gridInqParamGME(int gridID, int *nd, int *ni, int *ni2, int *ni3);
/* Lambert Conformal Conic grid (GRIB version) */
/* Lambert Conformal Conic grid */
void gridDefParamLCC(int gridID, double missval, double lon_0, double lat_0, double lat_1, double lat_2, double a, double rf, double xval_0, double yval_0, double x_0, double y_0);
int gridInqParamLCC(int gridID, double missval, double *lon_0, double *lat_0, double *lat_1, double *lat_2, double *a, double *rf, double *xval_0, double *yval_0, double *x_0, double *y_0);
/* Polar stereographic grid */
void gridDefParamSTERE(int gridID, double missval, double lon_0, double lat_ts, double lat_0, double a, double xval_0, double yval_0, double x_0, double y_0);
int gridInqParamSTERE(int gridID, double missval, double *lon_0, double *lat_ts, double *lat_0, double *a, double *xval_0, double *yval_0, double *x_0, double *y_0);
void gridDefArea(int gridID, const double area[]);
void gridInqArea(int gridID, double area[]);
int gridHasArea(int gridID);
......
src/grid.c
View file @ 554fd4c9
......@@ -21,6 +21,8 @@
double grid_missval = -9999.;
int (*proj_lonlat_to_lcc_func)() = NULL;
int (*proj_lcc_to_lonlat_func)() = NULL;
int (*proj_lonlat_to_stere_func)() = NULL;
int (*proj_stere_to_lonlat_func)() = NULL;
/* the value in the second pair of brackets must match the length of
* the longest string (including terminating NUL) */
......@@ -3780,9 +3782,9 @@ const char **gridInqYCvalsPtr(int gridID)
@Prototype void gridDefParamLCC(int gridID, double missval, double lon_0, double lat_0, double lat_1, double lat_2, double a, double rf, double xval_0, double yval_0, double x_0, double y_0)
@Parameter
@Item gridID Grid ID, from a previous call to @fref{gridCreate}.
@Item missval Missing value
@Item missval Missing value.
@Item lon_0 The East longitude of the meridian which is parallel to the Y-axis.
@Item lat_0 Latitude of the projection origin
@Item lat_0 Latitude of the projection origin.
@Item lat_1 First latitude from the pole at which the secant cone cuts the sphere.
@Item lat_2 Second latitude at which the secant cone cuts the sphere.
@Item a Earth radius in metres (optional).
......@@ -3835,7 +3837,7 @@ void gridDefParamLCC(int gridID, double missval, double lon_0, double lat_0, dou
@Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
@Item missval Missing value
@Item lon_0 The East longitude of the meridian which is parallel to the Y-axis.
@Item lat_0 Latitude of the projection origin
@Item lat_0 Latitude of the projection origin.
@Item lat_1 First latitude from the pole at which the secant cone cuts the sphere.
@Item lat_2 Second latitude at which the secant cone cuts the sphere.
@Item a Earth radius in metres (optional).
......@@ -3934,16 +3936,14 @@ int gridVerifyGribParamLCC(double missval, double *lon_0, double *lat_0, double
@Function gridDefParamSTERE
@Title Define the parameter of a Polar stereographic grid
@Prototype void gridDefParamSTERE(int gridID, double missval, double lon_0, double lat_0, double lat_1, double lat_2, double a, double rf, double xval_0, double yval_0, double x_0, double y_0)
@Prototype void gridDefParamSTERE(int gridID, double missval, double lon_0, double lat_ts, double lat_0, double a, double xval_0, double yval_0, double x_0, double y_0)
@Parameter
@Item gridID Grid ID, from a previous call to @fref{gridCreate}.
@Item missval Missing value
@Item lon_0 The East longitude of the meridian which is parallel to the Y-axis.
@Item lat_0 Latitude of the projection origin
@Item lat_1 First latitude from the pole at which the secant cone cuts the sphere.
@Item lat_2 Second latitude at which the secant cone cuts the sphere.
@Item missval Missing value.
@Item lon_0 Longitude at natural origin.
@Item lat_ts Latitude of the projection origin.
@Item lat_0 First latitude from the pole at which the secant cone cuts the sphere.
@Item a Earth radius in metres (optional).
@Item rf Inverse flattening (1/f) (optional).
@Item xval_0 Longitude of the first grid point in degree (optional).
@Item yval_0 Latitude of the first grid point in degree (optional).
@Item x_0 False easting (optional).
......@@ -3954,24 +3954,18 @@ The function @func{gridDefParamSTERE} defines the parameter of a Polar stereogra
@EndFunction
*/
void gridDefParamSTERE(int gridID, double missval, double lon_0, double lat_0, double lat_1, double lat_2,
double a, double rf, double xval_0, double yval_0, double x_0, double y_0)
void gridDefParamSTERE(int gridID, double missval, double lon_0, double lat_ts, double lat_0,
double a, double xval_0, double yval_0, double x_0, double y_0)
{
(void)lat_0;
cdiGridDefKeyStr(gridID, CDI_KEY_MAPPING, CDI_MAX_NAME, "Polar_Stereographic");
const char *mapname = "polar_stereographic";
cdiGridDefKeyStr(gridID, CDI_KEY_MAPNAME, CDI_MAX_NAME, mapname);
cdiDefAttTxt(gridID, CDI_GLOBAL, "grid_mapping_name", (int)(strlen(mapname)), mapname);
int nlats = 0;
double lats[2];
lats[nlats++] = lat_1;
if ( IS_NOT_EQUAL(lat_1, lat_2) ) lats[nlats++] = lat_2;
cdiDefAttFlt(gridID, CDI_GLOBAL, "standard_parallel", CDI_DATATYPE_FLT64, nlats, lats);
cdiDefAttFlt(gridID, CDI_GLOBAL, "longitude_of_central_meridian", CDI_DATATYPE_FLT64, 1, &lon_0);
cdiDefAttFlt(gridID, CDI_GLOBAL, "latitude_of_projection_origin", CDI_DATATYPE_FLT64, 1, &lat_2);
cdiDefAttFlt(gridID, CDI_GLOBAL, "standard_parallel", CDI_DATATYPE_FLT64, 1, &lat_ts);
cdiDefAttFlt(gridID, CDI_GLOBAL, "straight_vertical_longitude_from_pole", CDI_DATATYPE_FLT64, 1, &lon_0);
cdiDefAttFlt(gridID, CDI_GLOBAL, "latitude_of_projection_origin", CDI_DATATYPE_FLT64, 1, &lat_0);
if ( a > 0 ) cdiDefAttFlt(gridID, CDI_GLOBAL, "earth_radius", CDI_DATATYPE_FLT64, 1, &a);
if ( rf > 0 ) cdiDefAttFlt(gridID, CDI_GLOBAL, "inverse_flattening", CDI_DATATYPE_FLT64, 1, &rf);
if ( IS_NOT_EQUAL(x_0, missval) ) cdiDefAttFlt(gridID, CDI_GLOBAL, "false_easting", CDI_DATATYPE_FLT64, 1, &x_0);
if ( IS_NOT_EQUAL(y_0, missval) ) cdiDefAttFlt(gridID, CDI_GLOBAL, "false_northing", CDI_DATATYPE_FLT64, 1, &y_0);
if ( IS_NOT_EQUAL(xval_0, missval) ) cdiDefAttFlt(gridID, CDI_GLOBAL, "longitudeOfFirstGridPointInDegrees", CDI_DATATYPE_FLT64, 1, &xval_0);
......@@ -3987,16 +3981,14 @@ void gridDefParamSTERE(int gridID, double missval, double lon_0, double lat_0, d
@Function gridInqParamSTERE
@Title Get the parameter of a Polar stereographic grid
@Prototype void gridInqParamSTERE(int gridID, double missval, double *lon_0, double *lat_0, double *lat_1, double *lat_2, double *a, double *rf, double *xval_0, double *yval_0, double *x_0, double *y_0)
@Prototype void gridInqParamSTERE(int gridID, double missval, double *lon_0, double *lat_ts, double *lat_0, double *a, double *xval_0, double *yval_0, double *x_0, double *y_0)
@Parameter
@Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
@Item missval Missing value
@Item lon_0 The East longitude of the meridian which is parallel to the Y-axis.
@Item lat_0 Latitude of the projection origin
@Item lat_1 First latitude from the pole at which the secant cone cuts the sphere.
@Item lat_2 Second latitude at which the secant cone cuts the sphere.
@Item lon_0 Longitude at natural origin.
@Item lat_ts Latitude of the projection origin.
@Item lat_0 First latitude from the pole at which the secant cone cuts the sphere.
@Item a Earth radius in metres (optional).
@Item rf Inverse flattening (1/f) (optional).
@Item xval_0 Longitude of the first grid point in degree (optional).
@Item yval_0 Latitude of the first grid point in degree (optional).
@Item x_0 False easting (optional).
......@@ -4007,11 +3999,11 @@ The function @func{gridInqParamSTERE} returns the parameter of a Polar stereogra
@EndFunction
*/
int gridInqParamSTERE(int gridID, double missval, double *lon_0, double *lat_0, double *lat_1, double *lat_2,
double *a, double *rf, double *xval_0, double *yval_0, double *x_0, double *y_0)
int gridInqParamSTERE(int gridID, double missval, double *lon_0, double *lat_ts, double *lat_0,
double *a, double *xval_0, double *yval_0, double *x_0, double *y_0)
{
*a = 0; *rf = 0;
*lon_0 = missval; *lat_0 = missval, *lat_1 = missval, *lat_2 = missval;
*a = 0;
*lon_0 = missval; *lat_ts = missval, *lat_0 = missval;
*xval_0 = missval; *yval_0 = missval; *x_0 = missval, *y_0 = missval;
int status = -1;
......@@ -4039,19 +4031,14 @@ int gridInqParamSTERE(int gridID, double missval, double *lon_0, double *lat_0,
double attflt[2];
if ( cdiInqAttConvertedToFloat(gridID, atttype, attname, attlen, attflt) )
{
if ( strcmp(attname, "earth_radius") == 0 ) *a = attflt[0];
else if ( strcmp(attname, "inverse_flattening") == 0 ) *rf = attflt[0];
else if ( strcmp(attname, "longitude_of_central_meridian") == 0 ) *lon_0 = attflt[0];
else if ( strcmp(attname, "latitude_of_projection_origin") == 0 ) *lat_0 = attflt[0];
else if ( strcmp(attname, "false_easting") == 0 ) *x_0 = attflt[0];
else if ( strcmp(attname, "false_northing") == 0 ) *y_0 = attflt[0];
else if ( strcmp(attname, "longitudeOfFirstGridPointInDegrees") == 0 ) *xval_0 = attflt[0];
else if ( strcmp(attname, "latitudeOfFirstGridPointInDegrees") == 0 ) *yval_0 = attflt[0];
else if ( strcmp(attname, "standard_parallel") == 0 )
{
*lat_1 = attflt[0];
*lat_2 = (attlen == 2) ? attflt[1] : attflt[0];
}
if ( strcmp(attname, "earth_radius") == 0 ) *a = attflt[0];
else if ( strcmp(attname, "standard_parallel") == 0 ) *lat_ts = attflt[0];
else if ( strcmp(attname, "straight_vertical_longitude_from_pole") == 0 ) *lon_0 = attflt[0];
else if ( strcmp(attname, "latitude_of_projection_origin") == 0 ) *lat_0 = attflt[0];
else if ( strcmp(attname, "false_easting") == 0 ) *x_0 = attflt[0];
else if ( strcmp(attname, "false_northing") == 0 ) *y_0 = attflt[0];
else if ( strcmp(attname, "longitudeOfFirstGridPointInDegrees") == 0 ) *xval_0 = attflt[0];
else if ( strcmp(attname, "latitudeOfFirstGridPointInDegrees") == 0 ) *yval_0 = attflt[0];
}
}
}
......
src/grid.h
View file @ 554fd4c9
......@@ -9,6 +9,8 @@
extern double grid_missval;
extern int (*proj_lonlat_to_lcc_func)();
extern int (*proj_lcc_to_lonlat_func)();
extern int (*proj_lonlat_to_stere_func)();
extern int (*proj_stere_to_lonlat_func)();
typedef unsigned char mask_t;
......
src/stream_grb.c
View file @ 554fd4c9
......@@ -437,15 +437,20 @@ int grbGetGridtype(int gridID, size_t gridsize, bool *gridIsRotated, bool *gridI
if ( gridtype == GRID_PROJECTION )
{
if ( gridInqProjType(gridID) == CDI_PROJ_RLL )
int projtype = gridInqProjType(gridID);
if ( projtype == CDI_PROJ_RLL )
{
gridtype = GRID_LONLAT;
*gridIsRotated = true;
}
else if ( gridInqProjType(gridID) == CDI_PROJ_LCC )
else if ( projtype == CDI_PROJ_LCC )
{
gridtype = CDI_PROJ_LCC;
}
else if ( projtype == CDI_PROJ_STERE )
{
gridtype = CDI_PROJ_STERE;
}
}
return gridtype;
......
src/stream_gribapi.c
View file @ 554fd4c9
......@@ -540,36 +540,32 @@ void gribapiDefProjLCC(grib_handle *gh, int gridID)
static
void gribapiDefProjSTERE(grib_handle *gh, int gridID)
{
/*
double a = 6367470., rf = 0;
long earthIsOblate;
grib_get_long(gh, "earthIsOblate", &earthIsOblate);
if ( earthIsOblate ) { a = 6378160.; rf = 297.0; }
*/
double lon_0, lat_1, lat_2, xval_0, yval_0;
long projflag = 0, southPoleOnProjectionPlane;
long shapeOfTheEarth;
grib_get_long(gh, "shapeOfTheEarth", &shapeOfTheEarth);
double a = (shapeOfTheEarth == 6) ? 6371229 : 6367470;
double lon_0, lat_ts, xval_0, yval_0;
grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &xval_0);
grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &yval_0);
grib_get_double(gh, "LaDInDegrees", &lat_1);
grib_get_double(gh, "orientationOfTheGridInDegrees", &lat_2);
grib_get_long(gh, "projectionCentreFlag", &projflag);
grib_get_double(gh, "LaDInDegrees", &lat_ts);
grib_get_double(gh, "orientationOfTheGridInDegrees", &lon_0);
long southPoleOnProjectionPlane;
grib_get_long(gh, "southPoleOnProjectionPlane", &southPoleOnProjectionPlane);
bool lsouth = gribbyte_get_bit((int)projflag, 1);
if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }
double lat_0 = southPoleOnProjectionPlane ? -90 : 90;
double lat_0 = lat_2;
double x_0 = grid_missval;
double y_0 = grid_missval;
/*
if ( proj_lonlat_to_lcc_func )
if ( proj_lonlat_to_stere_func )
{
x_0 = xval_0; y_0 = yval_0;
proj_lonlat_to_lcc_func(grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, (size_t)1, &x_0, &y_0);
// proj_lonlat_to_stere_func(grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, (size_t)1, &x_0, &y_0);
if ( IS_NOT_EQUAL(x_0, grid_missval) && IS_NOT_EQUAL(y_0, grid_missval) )
{ x_0 = -x_0; y_0 = -y_0; }
}
gridDefParamLCC(gridID, grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0);
*/
gridDefParamSTERE(gridID, grid_missval, lon_0, lat_ts, lat_0, a, xval_0, yval_0, x_0, y_0);
}
static
......@@ -2042,14 +2038,10 @@ void gribapiDefGridRegular(grib_handle *gh, int gridID, int gridtype, bool gridI
GRIB_CHECK(my_grib_set_double(gh, "latitudeOfFirstGridPointInDegrees", yfirst), 0);
GRIB_CHECK(my_grib_set_double(gh, "latitudeOfLastGridPointInDegrees", ylast), 0);
{
long iscan = xfirst > xlast;
GRIB_CHECK(my_grib_set_long(gh, "iScansNegatively", iscan), 0);
}
{
long jscan = yfirst < ylast;
GRIB_CHECK(my_grib_set_long(gh, "jScansPositively", jscan), 0);
}
long xscan = xfirst > xlast;
GRIB_CHECK(my_grib_set_long(gh, "iScansNegatively", xscan), 0);
xscan = yfirst < ylast;
GRIB_CHECK(my_grib_set_long(gh, "jScansPositively", xscan), 0);
if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
{
......@@ -2080,10 +2072,10 @@ void gribapiDefGridRegular(grib_handle *gh, int gridID, int gridtype, bool gridI
}
static
void gribapiDefGridLambert(grib_handle *gh, int editionNumber, int gridID)
void gribapiDefGridLCC(grib_handle *gh, int editionNumber, int gridID)
{
int xsize = (int) gridInqXsize(gridID);
int ysize = (int) gridInqYsize(gridID);
long xsize = (long) gridInqXsize(gridID);
long ysize = (long) gridInqYsize(gridID);
double lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0;
gridInqParamLCC(gridID, grid_missval, &lon_0, &lat_0, &lat_1, &lat_2, &a, &rf, &xval_0, &yval_0, &x_0, &y_0);
......@@ -2123,6 +2115,45 @@ void gribapiDefGridLambert(grib_handle *gh, int editionNumber, int gridID)
GRIB_CHECK(my_grib_set_long(gh, "scanningMode", (long)scanflag), 0);
}
static
void gribapiDefGridSTERE(grib_handle *gh, int editionNumber, int gridID)
{
long xsize = (long) gridInqXsize(gridID);
long ysize = (long) gridInqYsize(gridID);
double lon_0, lat_ts, lat_0, a, xval_0, yval_0, x_0, y_0;
gridInqParamSTERE(gridID, grid_missval, &lon_0, &lat_ts, &lat_0, &a, &xval_0, &yval_0, &x_0, &y_0);
// gridVerifyGribParamSTERE(grid_missval, &lon_0, &lat_ts, &lat_0, &a, &xval_0, &yval_0, &x_0, &y_0);
if ( xval_0 < 0 ) xval_0 += 360;
int projflag = 0;
double xinc = gridInqXinc(gridID);
double yinc = gridInqYinc(gridID);
static const char mesg[] = "polar_stereographic";
size_t len = sizeof(mesg) -1;
GRIB_CHECK(my_grib_set_string(gh, "gridType", mesg, &len), 0);
GRIB_CHECK(my_grib_set_long(gh, "Nx", xsize), 0);
GRIB_CHECK(my_grib_set_long(gh, "Ny", ysize), 0);
GRIB_CHECK(my_grib_set_long(gh, "DxInMetres", lround(xinc)), 0);
GRIB_CHECK(my_grib_set_long(gh, "DyInMetres", lround(yinc)), 0);
GRIB_CHECK(my_grib_set_double(gh, "longitudeOfFirstGridPointInDegrees", xval_0), 0);
GRIB_CHECK(my_grib_set_double(gh, "latitudeOfFirstGridPointInDegrees", yval_0), 0);
GRIB_CHECK(my_grib_set_double(gh, "LaDInDegrees", lat_ts), 0);
GRIB_CHECK(my_grib_set_double(gh, "orientationOfTheGridInDegrees", lon_0), 0);
long southPoleOnProjectionPlane = IS_EQUAL(lat_0, -90.);
GRIB_CHECK(my_grib_set_double(gh, "southPoleOnProjectionPlane", southPoleOnProjectionPlane), 0);
GRIB_CHECK(my_grib_set_long(gh, "projectionCentreFlag", projflag), 0);
long shapeOfTheEarth = IS_EQUAL(a, 6371229.) ? 6 : 0;
if ( shapeOfTheEarth ) GRIB_CHECK(my_grib_set_long(gh, "shapeOfTheEarth", shapeOfTheEarth), 0);
GRIB_CHECK(my_grib_set_long(gh, "resolutionAndComponentFlags", 8), 0);
GRIB_CHECK(my_grib_set_long(gh, "iScansNegatively", 0), 0);
GRIB_CHECK(my_grib_set_long(gh, "jScansPositively", 1), 0);
}
static
void gribapiDefGridGME(grib_handle *gh, int gridID, long gridsize)
{
......@@ -2279,7 +2310,12 @@ void gribapiDefGrid(int editionNumber, grib_handle *gh, int gridID, int comptype
}
case CDI_PROJ_LCC:
{
gribapiDefGridLambert(gh, editionNumber, gridID);
gribapiDefGridLCC(gh, editionNumber, gridID);
break;
}
case CDI_PROJ_STERE:
{
gribapiDefGridSTERE(gh, editionNumber, gridID);
break;
}
case GRID_SPECTRAL:
......@@ -2665,7 +2701,8 @@ void gribapiSetScanningMode(grib_handle *gh, int scanningMode)
}
static void gribapiSetUvRelativeToGrid(grib_handle *gh, int mode)
static
void gribapiSetUvRelativeToGrid(grib_handle *gh, int mode)
{
long uvRelativeToGridMode = mode;
long uvRelativeToGridModeOld;
......@@ -2981,8 +3018,16 @@ void gribapiSetExtMode(grib_handle *gh, int gridID, size_t datasize, const doubl
(void)datasize;
#endif
int gridtype = gridInqType(gridID);
if ( gridtype == GRID_PROJECTION )
{
int projtype = gridInqProjType(gridID);
if ( projtype == CDI_PROJ_RLL ) gridtype = GRID_LONLAT;
else if ( projtype == CDI_PROJ_LCC ) gridtype = CDI_PROJ_LCC;
else if ( projtype == CDI_PROJ_STERE ) gridtype = CDI_PROJ_STERE;
}
if ( gridtype == GRID_GENERIC || gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN ||
gridtype == GRID_GAUSSIAN_REDUCED || gridtype == GRID_PROJECTION )
gridtype == GRID_GAUSSIAN_REDUCED || gridtype == CDI_PROJ_LCC )
{
#ifdef HIRLAM_EXTENSIONS
int scanModeIN = gridInqScanningMode(gridID);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment