Commit 983055f4 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added CDI_PROJ_STERE.

parent 8b45248d
......@@ -150,65 +150,66 @@ extern "C" {
/* Chunks */
#define CDI_CHUNK_AUTO 1 /* use default chunk size */
#define CDI_CHUNK_AUTO 1 // use default chunk size
#define CDI_CHUNK_GRID 2
#define CDI_CHUNK_LINES 3
/* GRID types */
#define GRID_GENERIC 1 /* Generic grid */
#define GRID_GAUSSIAN 2 /* Regular Gaussian lon/lat grid */
#define GRID_GAUSSIAN_REDUCED 3 /* Reduced Gaussian lon/lat grid */
#define GRID_LONLAT 4 /* Regular longitude/latitude grid */
#define GRID_SPECTRAL 5 /* Spherical harmonic coefficients */
#define GRID_FOURIER 6 /* Fourier coefficients */
#define GRID_GME 7 /* Icosahedral-hexagonal GME grid */
#define GRID_TRAJECTORY 8 /* Trajectory */
#define GRID_UNSTRUCTURED 9 /* General unstructured grid */
#define GRID_CURVILINEAR 10 /* Curvilinear grid */
#define GRID_LCC 11 /* Lambert Conformal Conic (GRIB) */
#define GRID_PROJECTION 12 /* Projected coordinates */
#define GRID_CHARXY 13 /* One horizontal character dimension */
#define CDI_PROJ_RLL 21 /* Rotated Latitude Longitude */
#define CDI_PROJ_LCC 22 /* Lambert Conformal Conic */
#define CDI_PROJ_LAEA 23 /* Lambert Azimuthal Equal Area */
#define CDI_PROJ_SINU 24 /* Sinusoidal */
#define GRID_GENERIC 1 // Generic grid
#define GRID_GAUSSIAN 2 // Regular Gaussian lon/lat grid
#define GRID_GAUSSIAN_REDUCED 3 // Reduced Gaussian lon/lat grid
#define GRID_LONLAT 4 // Regular longitude/latitude grid
#define GRID_SPECTRAL 5 // Spherical harmonic coefficients
#define GRID_FOURIER 6 // Fourier coefficients
#define GRID_GME 7 // Icosahedral-hexagonal GME grid
#define GRID_TRAJECTORY 8 // Trajectory
#define GRID_UNSTRUCTURED 9 // General unstructured grid
#define GRID_CURVILINEAR 10 // Curvilinear grid
#define GRID_LCC 11 // Lambert Conformal Conic (GRIB)
#define GRID_PROJECTION 12 // Projected coordinates
#define GRID_CHARXY 13 // One horizontal character dimension
#define CDI_PROJ_RLL 21 // Rotated Latitude Longitude
#define CDI_PROJ_LCC 22 // Lambert Conformal Conic
#define CDI_PROJ_LAEA 23 // Lambert Azimuthal Equal Area
#define CDI_PROJ_SINU 24 // Sinusoidal
#define CDI_PROJ_STERE 25 // Polar stereographic
/* ZAXIS types */
#define ZAXIS_SURFACE 0 /* Surface level */
#define ZAXIS_GENERIC 1 /* Generic level */
#define ZAXIS_HYBRID 2 /* Hybrid level */
#define ZAXIS_HYBRID_HALF 3 /* Hybrid half level */
#define ZAXIS_PRESSURE 4 /* Isobaric pressure level in Pascal */
#define ZAXIS_HEIGHT 5 /* Height above ground */
#define ZAXIS_DEPTH_BELOW_SEA 6 /* Depth below sea level in meters */
#define ZAXIS_DEPTH_BELOW_LAND 7 /* Depth below land surface in centimeters */
#define ZAXIS_ISENTROPIC 8 /* Isentropic */
#define ZAXIS_TRAJECTORY 9 /* Trajectory */
#define ZAXIS_ALTITUDE 10 /* Altitude above mean sea level in meters */
#define ZAXIS_SIGMA 11 /* Sigma level */
#define ZAXIS_MEANSEA 12 /* Mean sea level */
#define ZAXIS_TOA 13 /* Norminal top of atmosphere */
#define ZAXIS_SEA_BOTTOM 14 /* Sea bottom */
#define ZAXIS_ATMOSPHERE 15 /* Entire atmosphere */
#define ZAXIS_CLOUD_BASE 16 /* Cloud base level */
#define ZAXIS_CLOUD_TOP 17 /* Level of cloud tops */
#define ZAXIS_ISOTHERM_ZERO 18 /* Level of 0o C isotherm */
#define ZAXIS_SNOW 19 /* Snow level */
#define ZAXIS_LAKE_BOTTOM 20 /* Lake or River Bottom */
#define ZAXIS_SEDIMENT_BOTTOM 21 /* Bottom Of Sediment Layer */
#define ZAXIS_SEDIMENT_BOTTOM_TA 22 /* Bottom Of Thermally Active Sediment Layer */
#define ZAXIS_SEDIMENT_BOTTOM_TW 23 /* Bottom Of Sediment Layer Penetrated By Thermal Wave */
#define ZAXIS_MIX_LAYER 24 /* Mixing Layer */
#define ZAXIS_REFERENCE 25 /* zaxis reference number */
#define ZAXIS_CHAR 26 /* Area types */
#define ZAXIS_SURFACE 0 // Surface level
#define ZAXIS_GENERIC 1 // Generic level
#define ZAXIS_HYBRID 2 // Hybrid level
#define ZAXIS_HYBRID_HALF 3 // Hybrid half level
#define ZAXIS_PRESSURE 4 // Isobaric pressure level in Pascal
#define ZAXIS_HEIGHT 5 // Height above ground
#define ZAXIS_DEPTH_BELOW_SEA 6 // Depth below sea level in meters
#define ZAXIS_DEPTH_BELOW_LAND 7 // Depth below land surface in centimeters
#define ZAXIS_ISENTROPIC 8 // Isentropic
#define ZAXIS_TRAJECTORY 9 // Trajectory
#define ZAXIS_ALTITUDE 10 // Altitude above mean sea level in meters
#define ZAXIS_SIGMA 11 // Sigma level
#define ZAXIS_MEANSEA 12 // Mean sea level
#define ZAXIS_TOA 13 // Norminal top of atmosphere
#define ZAXIS_SEA_BOTTOM 14 // Sea bottom
#define ZAXIS_ATMOSPHERE 15 // Entire atmosphere
#define ZAXIS_CLOUD_BASE 16 // Cloud base level
#define ZAXIS_CLOUD_TOP 17 // Level of cloud tops
#define ZAXIS_ISOTHERM_ZERO 18 // Level of 0o C isotherm
#define ZAXIS_SNOW 19 // Snow level
#define ZAXIS_LAKE_BOTTOM 20 // Lake or River Bottom
#define ZAXIS_SEDIMENT_BOTTOM 21 // Bottom Of Sediment Layer
#define ZAXIS_SEDIMENT_BOTTOM_TA 22 // Bottom Of Thermally Active Sediment Layer
#define ZAXIS_SEDIMENT_BOTTOM_TW 23 // Bottom Of Sediment Layer Penetrated By Thermal Wave
#define ZAXIS_MIX_LAYER 24 // Mixing Layer
#define ZAXIS_REFERENCE 25 // zaxis reference number
#define ZAXIS_CHAR 26 // Area types
/* SUBTYPE types */
enum {
SUBTYPE_TILES = 0 /* Tiles variable */
SUBTYPE_TILES = 0 // Tiles variable
};
#define MAX_KV_PAIRS_MATCH 10
......
......@@ -32,26 +32,27 @@
#define GRIB2_LTYPE_ISENTROPIC 107
#define GRIB2_LTYPE_SNOW 114
#define GRIB2_LTYPE_REFERENCE 150
#define GRIB2_LTYPE_SEADEPTH 160 /* Depth Below Sea Level */
#define GRIB2_LTYPE_LAKE_BOTTOM 162 /* Lake or River Bottom */
#define GRIB2_LTYPE_SEDIMENT_BOTTOM 163 /* Bottom Of Sediment Layer */
#define GRIB2_LTYPE_SEDIMENT_BOTTOM_TA 164 /* Bottom Of Thermally Active Sediment Layer */
#define GRIB2_LTYPE_SEDIMENT_BOTTOM_TW 165 /* Bottom Of Sediment Layer Penetrated By Thermal Wave */
#define GRIB2_LTYPE_MIX_LAYER 166 /* Mixing Layer */
#define GRIB2_LTYPE_SEADEPTH 160 // Depth Below Sea Level
#define GRIB2_LTYPE_LAKE_BOTTOM 162 // Lake or River Bottom
#define GRIB2_LTYPE_SEDIMENT_BOTTOM 163 // Bottom Of Sediment Layer
#define GRIB2_LTYPE_SEDIMENT_BOTTOM_TA 164 // Bottom Of Thermally Active Sediment Layer
#define GRIB2_LTYPE_SEDIMENT_BOTTOM_TW 165 // Bottom Of Sediment Layer Penetrated By Thermal Wave
#define GRIB2_LTYPE_MIX_LAYER 166 // Mixing Layer
/* GRIB2 Data representation type (Grid Type) */
#define GRIB2_GTYPE_LATLON 0 /* latitude/longitude */
#define GRIB2_GTYPE_LATLON_ROT 1 /* rotated latitude/longitude */
#define GRIB2_GTYPE_LATLON_STR 2 /* stretched latitude/longitude */
#define GRIB2_GTYPE_LATLON_ROTSTR 3 /* rotated and stretched latitude/longitude */
#define GRIB2_GTYPE_GAUSSIAN 40 /* gaussian grid */
#define GRIB2_GTYPE_GAUSSIAN_ROT 41 /* rotated gaussian grid */
#define GRIB2_GTYPE_GAUSSIAN_STR 42 /* stretched gaussian grid */
#define GRIB2_GTYPE_GAUSSIAN_ROTSTR 43 /* rotated and stretched gaussian grid */
#define GRIB2_GTYPE_LCC 30 /* Lambert conformal */
#define GRIB2_GTYPE_SPECTRAL 50 /* spherical harmonics */
#define GRIB2_GTYPE_GME 100 /* hexagonal GME grid */
#define GRIB2_GTYPE_UNSTRUCTURED 101 /* General Unstructured Grid */
#define GRIB2_GTYPE_LATLON 0 // Latitude/longitude (or equidistant cylindrical, or Plate Carree)
#define GRIB2_GTYPE_LATLON_ROT 1 // Rotated Latitude/longitude
#define GRIB2_GTYPE_LATLON_STR 2 // Stretched Latitude/longitude
#define GRIB2_GTYPE_LATLON_ROTSTR 3 // Stretched and Rotated Latitude/longitude
#define GRIB2_GTYPE_STERE 20 // Polar stereographic projection
#define GRIB2_GTYPE_LCC 30 // Lambert conformal
#define GRIB2_GTYPE_GAUSSIAN 40 // Gaussian latitude/longitude
#define GRIB2_GTYPE_GAUSSIAN_ROT 41 // Rotated Gaussian latitude/longitude
#define GRIB2_GTYPE_GAUSSIAN_STR 42 // Stretched Gaussian latitude/longitude
#define GRIB2_GTYPE_GAUSSIAN_ROTSTR 43 // Stretched and rotated Gaussian latitude/longitude
#define GRIB2_GTYPE_SPECTRAL 50 // Spherical harmonic coefficients
#define GRIB2_GTYPE_GME 100 // Triangular grid based on an icosahedron (GME)
#define GRIB2_GTYPE_UNSTRUCTURED 101 // General Unstructured Grid
const char *gribapiLibraryVersionString(void);
void gribContainersNew(stream_t * streamptr);
......
This diff is collapsed.
......@@ -3930,6 +3930,135 @@ int gridVerifyGribParamLCC(double missval, double *lon_0, double *lat_0, double
return 0;
}
/*
@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)
@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 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).
@Item y_0 False northing (optional).
@Description
The function @func{gridDefParamSTERE} defines the parameter of a Polar stereographic grid.
@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)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);
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);
if ( IS_NOT_EQUAL(yval_0, missval) ) cdiDefAttFlt(gridID, CDI_GLOBAL, "latitudeOfFirstGridPointInDegrees", CDI_DATATYPE_FLT64, 1, &yval_0);
grid_t *gridptr = grid_to_pointer(gridID);
gridptr->projtype = CDI_PROJ_STERE;
gridVerifyProj(gridID);
}
/*
@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)
@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 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).
@Item y_0 False northing (optional).
@Description
The function @func{gridInqParamSTERE} returns the parameter of a Polar stereographic grid.
@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)
{
*a = 0; *rf = 0;
*lon_0 = missval; *lat_0 = missval, *lat_1 = missval, *lat_2 = missval;
*xval_0 = missval; *yval_0 = missval; *x_0 = missval, *y_0 = missval;
int status = -1;
if ( gridInqType(gridID) != GRID_PROJECTION ) return status;
status = -2;
const char *projection = "polar_stereographic";
char mapname[CDI_MAX_NAME]; mapname[0] = 0;
cdiGridInqKeyStr(gridID, CDI_KEY_MAPNAME, CDI_MAX_NAME, mapname);
if ( mapname[0] && strcmp(mapname, projection) == 0 )
{
int atttype, attlen;
char attname[CDI_MAX_NAME+1];
int natts;
cdiInqNatts(gridID, CDI_GLOBAL, &natts);
if ( natts ) status = 0;
for ( int iatt = 0; iatt < natts; ++iatt )
{
cdiInqAtt(gridID, CDI_GLOBAL, iatt, attname, &atttype, &attlen);
if ( attlen > 2 ) continue;
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];
}
}
}
}
return status;
}
void gridDefComplexPacking(int gridID, int lcomplex)
{
......@@ -4860,8 +4989,7 @@ struct addIfNewRes cdiVlistAddGridIfNew(int vlistID, grid_t *grid, int mode)
}
}
else
Error("Internal problem: undefined gridID in vlist "
"%d, position %u!", vlistID, index);
Error("Internal problem: undefined gridID in vlist %d, position %u!", vlistID, index);
}
if ( ! griddefined )
......
......@@ -489,6 +489,89 @@ void gribapiGetKeys(grib_handle *gh, int varID)
}
}
static
void gribapiDefProjRLL(grib_handle *gh, int gridID)
{
double xpole = 0, ypole = 0, angle = 0;
grib_get_double(gh, "latitudeOfSouthernPoleInDegrees", &ypole);
grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &xpole);
grib_get_double(gh, "angleOfRotation", &angle);
xpole -= 180;
if ( fabs(ypole) > 0 ) ypole = -ypole; // change from south to north pole
if ( fabs(angle) > 0 ) angle = -angle;
gridDefParamRLL(gridID, xpole, ypole, angle);
}
static
void gribapiDefProjLCC(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;
grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &xval_0);
grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &yval_0);
grib_get_double(gh, "LoVInDegrees", &lon_0);
grib_get_double(gh, "Latin1InDegrees", &lat_1);
grib_get_double(gh, "Latin2InDegrees", &lat_2);
grib_get_long(gh, "projectionCentreFlag", &projflag);
bool lsouth = gribbyte_get_bit((int)projflag, 1);
if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }
double lat_0 = lat_2;
double x_0 = grid_missval;
double y_0 = grid_missval;
if ( proj_lonlat_to_lcc_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);
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);
}
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;
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_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 = lat_2;
double x_0 = grid_missval;
double y_0 = grid_missval;
/*
if ( proj_lonlat_to_lcc_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);
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);
*/
}
static
void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
size_t recsize, off_t position, int datatype, int comptype, const char *varname,
......@@ -535,49 +618,10 @@ void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh,
struct addIfNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
int gridID = gridAdded.Id;
if ( !gridAdded.isNew ) Free(grid);
else if ( grid->projtype == CDI_PROJ_RLL )
{
double xpole = 0, ypole = 0, angle = 0;
grib_get_double(gh, "latitudeOfSouthernPoleInDegrees", &ypole);
grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &xpole);
grib_get_double(gh, "angleOfRotation", &angle);
xpole -= 180;
if ( fabs(ypole) > 0 ) ypole = -ypole; // change from south to north pole
if ( fabs(angle) > 0 ) angle = -angle;
gridDefParamRLL(gridID, xpole, ypole, angle);
}
else if ( grid->projtype == CDI_PROJ_LCC )
{
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;
grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &xval_0);
grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &yval_0);
grib_get_double(gh, "LoVInDegrees", &lon_0);
grib_get_double(gh, "Latin1InDegrees", &lat_1);
grib_get_double(gh, "Latin2InDegrees", &lat_2);
grib_get_long(gh, "projectionCentreFlag", &projflag);
bool lsouth = gribbyte_get_bit((int)projflag, 1);
if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }
double lat_0 = lat_2;
double x_0 = grid_missval;
double y_0 = grid_missval;
if ( proj_lonlat_to_lcc_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);
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);
}
if ( !gridAdded.isNew ) Free(grid);
else if ( grid->projtype == CDI_PROJ_RLL ) gribapiDefProjRLL(gh, gridID);
else if ( grid->projtype == CDI_PROJ_LCC ) gribapiDefProjLCC(gh, gridID);
else if ( grid->projtype == CDI_PROJ_STERE ) gribapiDefProjSTERE(gh, gridID);
int zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1);
......
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