Commit 072e3ce3 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Use CDO function proj_lonlat_to_lcc() to convert lonlat to lcc coordinates.

parent 316b9569
......@@ -21,6 +21,7 @@
#include "vlist.h"
double grid_missval = -9999.;
int (*proj_lonlat_to_lcc_func)() = NULL;
/* the value in the second pair of brackets must match the length of
* the longest string (including terminating NUL) */
......@@ -3627,10 +3628,10 @@ void gridDefParamLCC(int gridID, double missval, double lon_0, double lat_0, dou
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(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);
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_LCC;
......
......@@ -7,6 +7,7 @@
#include "cdi_att.h"
extern double grid_missval;
extern int (*proj_lonlat_to_lcc_func)();
typedef unsigned char mask_t;
......
......@@ -406,6 +406,14 @@ void cgribexAddRecord(stream_t * streamptr, int param, int *isec1, int *isec2, d
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, 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);
}
}
......
......@@ -492,7 +492,7 @@ void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
long earthIsOblate;
grib_get_long(gh, "earthIsOblate", &earthIsOblate);
if ( earthIsOblate ) { a = 6378160.; rf = 297.0; }
double lon_0, lat_0, lat_1, lat_2, xval_0, yval_0, x_0 = grid_missval, y_0 = grid_missval;
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);
......@@ -503,7 +503,17 @@ void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
bool lsouth = gribbyte_get_bit(projflag, 1);
if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }
lat_0 = 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, 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);
}
......
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