Commit 7dabeb17 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Importcmsaf.c: added support for CDI_PROJ_LAEA.

parent 720da445
...@@ -30,6 +30,7 @@ ...@@ -30,6 +30,7 @@
#include <cdi.h> #include <cdi.h>
#include "cdo.h" #include "cdo.h"
#include "cdo_int.h" #include "cdo_int.h"
#include "grid.h"
#include "pstream.h" #include "pstream.h"
...@@ -266,24 +267,25 @@ int defLaeaGrid(int nx, int ny, double xmin, double xmax, double ymin, double ym ...@@ -266,24 +267,25 @@ int defLaeaGrid(int nx, int ny, double xmin, double xmax, double ymin, double ym
double *yvals = (double*) Malloc(ny*sizeof(double)); double *yvals = (double*) Malloc(ny*sizeof(double));
for ( int i = 0; i < nx; ++i ) for ( int i = 0; i < nx; ++i )
{ xvals[i] = xmin + i*dx + dx/2;
xvals[i] = xmin + i*dx + dx/2;
/* printf("x[%d]=%g\n", i, xvals[i]); */
}
for ( int i = 0; i < ny; ++i ) for ( int i = 0; i < ny; ++i )
{ yvals[i] = ymax - i*dy - dy/2;
yvals[i] = ymax - i*dy - dy/2;
/* printf("y[%d]=%g\n", i, yvals[i]); */
}
#ifdef TEST_PROJECTION
int gridID = gridCreate(GRID_PROJECTION, nx*ny);
#else
int gridID = gridCreate(GRID_LAEA, nx*ny); int gridID = gridCreate(GRID_LAEA, nx*ny);
#endif
gridDefXsize(gridID, nx); gridDefXsize(gridID, nx);
gridDefYsize(gridID, ny); gridDefYsize(gridID, ny);
gridDefXvals(gridID, xvals); gridDefXvals(gridID, xvals);
gridDefYvals(gridID, yvals); gridDefYvals(gridID, yvals);
#ifdef TEST_PROJECTION
grid_def_param_laea(gridID, a, lon0, lat0);
#else
gridDefLaea(gridID, a, lon0, lat0); gridDefLaea(gridID, a, lon0, lat0);
#endif
Free(xvals); Free(xvals);
Free(yvals); Free(yvals);
......
...@@ -491,8 +491,21 @@ void sinusoidal_to_geo(int gridsize, double *xvals, double *yvals) ...@@ -491,8 +491,21 @@ void sinusoidal_to_geo(int gridsize, double *xvals, double *yvals)
#endif #endif
} }
void grid_def_param_laea(int gridID, double a, double lon_0, double lat_0)
{
const char *projection = "lambert_azimuthal_equal_area";
cdiGridDefKeyStr(gridID, CDI_KEY_MAPPING, (int)(strlen(projection)+1), projection);
const char *mapvarname = "Lambert_AEA";
cdiGridDefKeyStr(gridID, CDI_KEY_MAPNAME, (int)(strlen(mapvarname)+1), mapvarname);
vlistDefAttFlt(gridID, CDI_GLOBAL, "earth_radius", DATATYPE_FLT64, 1, &a);
vlistDefAttFlt(gridID, CDI_GLOBAL, "longitude_of_projection_origin", DATATYPE_FLT64, 1, &lon_0);
vlistDefAttFlt(gridID, CDI_GLOBAL, "latitude_of_projection_origin", DATATYPE_FLT64, 1, &lat_0);
}
static static
void grid_get_param_laea(int gridID, double *a, double *lon_0, double *lat_0) void grid_inq_param_laea(int gridID, double *a, double *lon_0, double *lat_0)
{ {
*a = 0; *lon_0 = 0; *lat_0 = 0; *a = 0; *lon_0 = 0; *lat_0 = 0;
...@@ -541,7 +554,7 @@ void laea_to_geo(int gridID, int gridsize, double *xvals, double *yvals) ...@@ -541,7 +554,7 @@ void laea_to_geo(int gridID, int gridsize, double *xvals, double *yvals)
projUV data, res; projUV data, res;
double a, lon_0, lat_0; double a, lon_0, lat_0;
grid_get_param_laea(gridID, &a , &lon_0, &lat_0); grid_inq_param_laea(gridID, &a , &lon_0, &lat_0);
int nbpar = 0; int nbpar = 0;
params[nbpar++] = gen_param("proj=laea"); params[nbpar++] = gen_param("proj=laea");
......
...@@ -21,6 +21,7 @@ ...@@ -21,6 +21,7 @@
#define DEG2RAD (M_PI/180.) /* conversion for deg to rad */ #define DEG2RAD (M_PI/180.) /* conversion for deg to rad */
#endif #endif
void grid_def_param_laea(int gridID, double a, double lon0, double lat0);
bool grid_is_distance_generic(int gridID); bool grid_is_distance_generic(int gridID);
......
Supports Markdown
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