Commit 26e6ca5e authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

add support for GRID type Lambert Conformal Conic

parent 450e824c
2009-02-?? Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* using CDI library version 1.3.1
* add support for GRID type Lambert Conformal Conic
* Sinusoidal and Lambert Azimuthal Equal Area grids with units [km]
* Version 1.3.1 released
......
#! /bin/sh
# Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.61 for cdo 1.3.0.
# Generated by GNU Autoconf 2.61 for cdo 1.3.1.
#
# Report bugs to <Uwe.Schulzweida@zmaw.de>.
#
......@@ -574,8 +574,8 @@ SHELL=${CONFIG_SHELL-/bin/sh}
# Identity of this package.
PACKAGE_NAME='cdo'
PACKAGE_TARNAME='cdo'
PACKAGE_VERSION='1.3.0'
PACKAGE_STRING='cdo 1.3.0'
PACKAGE_VERSION='1.3.1'
PACKAGE_STRING='cdo 1.3.1'
PACKAGE_BUGREPORT='Uwe.Schulzweida@zmaw.de'
# Factoring default headers for most tests.
......@@ -1222,7 +1222,7 @@ if test "$ac_init_help" = "long"; then
# Omit some internal or obsolete options to make the list less imposing.
# This message is too long to be a string in the A/UX 3.1 sh.
cat <<_ACEOF
\`configure' configures cdo 1.3.0 to adapt to many kinds of systems.
\`configure' configures cdo 1.3.1 to adapt to many kinds of systems.
Usage: $0 [OPTION]... [VAR=VALUE]...
......@@ -1293,7 +1293,7 @@ fi
if test -n "$ac_init_help"; then
case $ac_init_help in
short | recursive ) echo "Configuration of cdo 1.3.0:";;
short | recursive ) echo "Configuration of cdo 1.3.1:";;
esac
cat <<\_ACEOF
......@@ -1404,7 +1404,7 @@ fi
test -n "$ac_init_help" && exit $ac_status
if $ac_init_version; then
cat <<\_ACEOF
cdo configure 1.3.0
cdo configure 1.3.1
generated by GNU Autoconf 2.61
Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
......@@ -1418,7 +1418,7 @@ cat >config.log <<_ACEOF
This file contains any messages produced by compilers while
running configure, to aid debugging if configure makes a mistake.
It was created by cdo $as_me 1.3.0, which was
It was created by cdo $as_me 1.3.1, which was
generated by GNU Autoconf 2.61. Invocation command line was
$ $0 $@
......@@ -2242,7 +2242,7 @@ fi
# Define the identity of the package.
PACKAGE='cdo'
VERSION='1.3.0'
VERSION='1.3.1'
cat >>confdefs.h <<_ACEOF
......@@ -7772,7 +7772,7 @@ exec 6>&1
# report actual input values of CONFIG_FILES etc. instead of their
# values after options handling.
ac_log="
This file was extended by cdo $as_me 1.3.0, which was
This file was extended by cdo $as_me 1.3.1, which was
generated by GNU Autoconf 2.61. Invocation command line was
CONFIG_FILES = $CONFIG_FILES
......@@ -7825,7 +7825,7 @@ Report bugs to <bug-autoconf@gnu.org>."
_ACEOF
cat >>$CONFIG_STATUS <<_ACEOF
ac_cs_version="\\
cdo config.status 1.3.0
cdo config.status 1.3.1
configured by $0, generated by GNU Autoconf 2.61,
with options \\"`echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"
......
# Process this file with autoconf to produce a configure script.
AC_INIT(cdo, 1.3.0, Uwe.Schulzweida@zmaw.de)
AC_INIT(cdo, 1.3.1, Uwe.Schulzweida@zmaw.de)
CONFIG_ABORT=yes
......
......@@ -129,6 +129,6 @@ Say @file{ifile} contains fields on a regular quadrilateral grid.
To remap all fields bilinear to a T42 Gaussian grid use:
@BeginVerbatim
cdo genbil,t42grid ifile remapweights.nc
cdo remap,t42grid,remapweights ifile ofile
cdo remap,t42grid,remapweights.nc ifile ofile
@EndVerbatim
@EndExample
......@@ -42,6 +42,6 @@ Say @file{ifile} contains fields on a regular quadrilateral grid.
To remap all fields bilinear to a T42 Gaussian grid use:
@BeginVerbatim
cdo genbil,t42grid ifile remapweights.nc
cdo remap,t42grid,remapweights ifile ofile
cdo remap,t42grid,remapweights.nc ifile ofile
@EndVerbatim
@EndExample
......@@ -157,7 +157,7 @@ void *Ensstat(void *argument)
streamID = ef[fileID].streamID;
nrecs = streamInqTimestep(streamID, tsID);
if ( nrecs != nrecs0 )
cdoAbort("Number of records change from %d to %d", nrecs0, nrecs);
cdoAbort("Number of records changed from %d to %d", nrecs0, nrecs);
}
taxisCopyTimestep(taxisID2, taxisID1);
......
......@@ -129,6 +129,7 @@ static void printGridInfo(int vlistID)
gridID+1, gridNamePtr(gridtype));
if ( gridtype == GRID_LONLAT ||
gridtype == GRID_LCC2 ||
gridtype == GRID_LAEA ||
gridtype == GRID_SINUSOIDAL ||
gridtype == GRID_GAUSSIAN ||
......@@ -165,7 +166,8 @@ static void printGridInfo(int vlistID)
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "%-9s : first = %.9g last = %.9g", yname, yfirst, ylast);
if ( !DBL_IS_EQUAL(yinc, 0) &&
(gridtype == GRID_LONLAT || gridtype == GRID_SINUSOIDAL || gridtype == GRID_LAEA) )
(gridtype == GRID_LONLAT || gridtype == GRID_SINUSOIDAL ||
gridtype == GRID_LCC2 || gridtype == GRID_LAEA) )
fprintf(stdout, " inc = %.9g", yinc);
fprintf(stdout, " %s", yunits);
fprintf(stdout, "\n");
......@@ -195,6 +197,15 @@ static void printGridInfo(int vlistID)
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "projpar : a = %g lon_0 = %g lat_0 = %g\n", a, lon_0, lat_0);
}
if ( gridtype == GRID_LCC2 )
{
double a, lon_0, lat_0, lat_1, lat_2;
gridInqLcc2(gridID, &a, &lon_0, &lat_0, &lat_1, &lat_2);
fprintf(stdout, "%*s", nbyte0, "");
fprintf(stdout, "projpar : a = %7.0f lon_0 = %g lat_0 = %g lat_1 = %g lat_2 = %g\n",
a, lon_0, lat_0, lat_1, lat_2);
}
}
else if ( gridtype == GRID_SPECTRAL )
{
......
......@@ -107,9 +107,10 @@ extern "C" {
#define GRID_TRAJECTORY 8
#define GRID_CELL 9
#define GRID_CURVILINEAR 10
#define GRID_LCC 11 /* Lambert Conformal Conic */
#define GRID_LAEA 12 /* Lambert Azimuthal Equal Area */
#define GRID_SINUSOIDAL 13 /* Sinusoidal */
#define GRID_LCC 11 /* Lambert Conformal Conic (GRIB) */
#define GRID_LCC2 12 /* Lambert Conformal Conic (PROJ) */
#define GRID_LAEA 13 /* Lambert Azimuthal Equal Area */
#define GRID_SINUSOIDAL 14 /* Sinusoidal */
/* ZAXIS types */
......@@ -544,10 +545,14 @@ void gridDefGMEni2(int gridID, int ni2);
int gridInqGMEni3(int gridID);
void gridDefGMEni3(int gridID, int ni3);
/* Lambert Conformal Conic grid */
/* Lambert Conformal Conic grid (GRIB version) */
void gridDefLCC(int gridID, double originLon, double originLat, double lonParY, double lat1, double lat2, double xinc, double yinc, int projflag, int scanflag);
void gridInqLCC(int gridID, double *originLon, double *originLat, double *lonParY, double *lat1, double *lat2, double *xinc, double *yinc, int *projflag, int *scanflag);
/* Lambert Conformal Conic 2 grid (PROJ version) */
void gridDefLcc2(int gridID, double earth_radius, double lon_0, double lat_0, double lat_1, double lat_2);
void gridInqLcc2(int gridID, double *earth_radius, double *lon_0, double *lat_0, double *lat_1, double *lat_2);
/* Lambert Azimuthal Equal Area grid */
void gridDefLaea(int gridID, double earth_radius, double lon_0, double lat_0);
void gridInqLaea(int gridID, double *earth_radius, double *lon_0, double *lat_0);
......
......@@ -308,6 +308,7 @@ int gridToCurvilinear(int gridID1)
case GRID_LONLAT:
case GRID_GAUSSIAN:
case GRID_LCC:
case GRID_LCC2:
case GRID_LAEA:
case GRID_SINUSOIDAL:
{
......@@ -472,6 +473,51 @@ int gridToCurvilinear(int gridID1)
cdoAbort("proj4 support not compiled in!");
#endif
}
else if ( gridtype == GRID_LCC2 )
{
#if defined (HAVE_LIBPROJ)
int i;
PJ *libProj;
char *params[20];
int nbpar=0;
projUV data, res;
double a, lon_0, lat_0, lat_1, lat_2;
gridInqLcc2(gridID1, &a , &lon_0, &lat_0, &lat_1, &lat_2);
nbpar = 0;
params[nbpar++] = gen_param("proj=lcc");
if ( a > 0 ) params[nbpar++] = gen_param("a=%g", a);
params[nbpar++] = gen_param("lon_0=%g", lon_0);
params[nbpar++] = gen_param("lat_0=%g", lat_0);
params[nbpar++] = gen_param("lat_1=%g", lat_1);
params[nbpar++] = gen_param("lat_2=%g", lat_2);
if ( cdoVerbose )
for ( i = 0; i < nbpar; ++i )
cdoPrint("Proj.param[%d] = %s", i+1, params[i]);
libProj = pj_init(nbpar, &params[0]);
if ( !libProj )
cdoAbort("proj error: %s", pj_strerrno(pj_errno));
for ( i = 0; i < nbpar; ++i ) free(params[i]);
/* libProj->over = 1; */ // allow longitude > 180
for ( j = 0; j < ny; j++ )
for ( i = 0; i < nx; i++ )
{
data.u = xscale*xvals[i];
data.v = yscale*yvals[j];
res = pj_inv(data, libProj);
xvals2D[j*nx+i] = res.u*rad2deg;
yvals2D[j*nx+i] = res.v*rad2deg;
}
#else
cdoAbort("proj4 support not compiled in!");
#endif
}
else
{
for ( j = 0; j < ny; j++ )
......@@ -562,7 +608,7 @@ int gridToCurvilinear(int gridID1)
else if ( ny > 1 )
{
ybounds = (double *) malloc(2*ny*sizeof(double));
if ( gridtype == GRID_SINUSOIDAL || gridtype == GRID_LAEA )
if ( gridtype == GRID_SINUSOIDAL || gridtype == GRID_LAEA || gridtype == GRID_LCC2 )
gridGenYboundsM(ny, yvals, ybounds);
else
gridGenYbounds(ny, yvals, ybounds);
......@@ -709,6 +755,75 @@ int gridToCurvilinear(int gridID1)
cdoAbort("proj4 support not compiled in!");
#endif
}
else if ( gridtype == GRID_LCC2 )
{
#if defined (HAVE_LIBPROJ)
int index;
int i;
PJ *libProj;
char *params[20];
int nbpar=0;
projUV data, res;
double a, lon_0, lat_0, lat_1, lat_2;
gridInqLcc2(gridID1, &a , &lon_0, &lat_0, &lat_1, &lat_2);
nbpar = 0;
params[nbpar++] = gen_param("proj=lcc");
if ( a > 0 ) params[nbpar++] = gen_param("a=%g", a);
params[nbpar++] = gen_param("lon_0=%g", lon_0);
params[nbpar++] = gen_param("lat_0=%g", lat_0);
params[nbpar++] = gen_param("lat_1=%g", lat_1);
params[nbpar++] = gen_param("lat_2=%g", lat_2);
if ( cdoVerbose )
for ( i = 0; i < nbpar; ++i )
cdoPrint("Proj.param[%d] = %s", i+1, params[i]);
libProj = pj_init(nbpar, params);
if ( !libProj )
cdoAbort("proj error: %s", pj_strerrno(pj_errno));
for ( i = 0; i < nbpar; ++i ) free(params[i]);
/* libProj->over = 1; */ // allow longitude > 180
/*
for ( j = 0; j < 2*ny; j++ ) printf("ybounds %d %g\n", j, ybounds[j]);
for ( i = 0; i < 2*nx; i++ ) printf("xbounds %d %g\n", i, xbounds[i]);
*/
for ( j = 0; j < ny; j++ )
for ( i = 0; i < nx; i++ )
{
index = j*4*nx + 4*i;
data.u = xscale*xbounds[2*i];
data.v = yscale*ybounds[2*j];
res = pj_inv(data, libProj);
xbounds2D[index+0] = res.u*rad2deg;
ybounds2D[index+0] = res.v*rad2deg;
data.u = xscale*xbounds[2*i];
data.v = yscale*ybounds[2*j+1];
res = pj_inv(data, libProj);
xbounds2D[index+1] = res.u*rad2deg;
ybounds2D[index+1] = res.v*rad2deg;
data.u = xscale*xbounds[2*i+1];
data.v = yscale*ybounds[2*j+1];
res = pj_inv(data, libProj);
xbounds2D[index+2] = res.u*rad2deg;
ybounds2D[index+2] = res.v*rad2deg;
data.u = xscale*xbounds[2*i+1];
data.v = yscale*ybounds[2*j];
res = pj_inv(data, libProj);
xbounds2D[index+3] = res.u*rad2deg;
ybounds2D[index+3] = res.v*rad2deg;
}
#else
cdoAbort("proj4 support not compiled in!");
#endif
}
else
{
gridGenXbounds2D(nx, ny, xbounds, xbounds2D);
......
......@@ -2,7 +2,7 @@
This file is part of CDO. CDO is a collection of Operators to
manipulate and analyse Climate model Data.
Copyright (C) 2003-2008 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
Copyright (C) 2003-2009 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
See COPYING file for copying and redistribution conditions.
This program is free software; you can redistribute it and/or modify
......@@ -85,6 +85,16 @@ void gridInit(GRID *grid)
grid->def_lat1 = FALSE;
grid->def_lat2 = FALSE;
grid->a = 0;
grid->lon_0 = 0;
grid->lat_0 = 0;
grid->lat_1 = 0;
grid->lat_2 = 0;
grid->def_lon_0 = FALSE;
grid->def_lat_0 = FALSE;
grid->def_lat_1 = FALSE;
grid->def_lat_2 = FALSE;
grid->def_xfirst = FALSE;
grid->def_yfirst = FALSE;
grid->def_xlast = FALSE;
......@@ -334,6 +344,55 @@ int gridDefine(GRID grid)
gridDefLCC(gridID, grid.originLon, grid.originLat, grid.lonParY,
grid.lat1, grid.lat2, grid.xinc, grid.yinc, grid.projflag, grid.scanflag);
break;
}
case GRID_LCC2:
{
if ( grid.xsize == 0 ) Error(func, "xsize undefined!");
if ( grid.ysize == 0 ) Error(func, "ysize undefined!");
if ( grid.size == 0 ) grid.size = grid.xsize*grid.ysize;
gridID = gridCreate(grid.type, grid.size);
gridDefPrec(gridID, grid.prec);
gridDefXsize(gridID, grid.xsize);
gridDefYsize(gridID, grid.ysize);
if ( grid.def_xfirst && grid.def_xinc && grid.xvals == NULL )
{
grid.xvals = (double *) malloc(grid.xsize*sizeof(double));
for ( i = 0; i < grid.xsize; ++i )
grid.xvals[i] = grid.xfirst + i*grid.xinc;
}
if ( grid.def_yfirst && grid.def_yinc && grid.yvals == NULL )
{
grid.yvals = (double *) malloc(grid.ysize*sizeof(double));
for ( i = 0; i < grid.ysize; ++i )
grid.yvals[i] = grid.yfirst + i*grid.yinc;
}
if ( grid.xvals )
{
gridDefXvals(gridID, grid.xvals);
free(grid.xvals);
}
if ( grid.yvals )
{
gridDefYvals(gridID, grid.yvals);
free(grid.yvals);
}
if ( grid.def_lon_0 == FALSE ) Error(func, "lon_0 undefined!");
if ( grid.def_lat_0 == FALSE ) Error(func, "lat_0 undefined!");
if ( grid.def_lat_1 == FALSE ) Error(func, "lat_1 undefined!");
if ( grid.def_lat_2 == FALSE ) grid.def_lat_2 = grid.def_lat_1;
gridDefLcc2(gridID, grid.a, grid.lon_0, grid.lat_0, grid.lat_1, grid.lat_2);
break;
}
case GRID_SPECTRAL:
......@@ -567,10 +626,12 @@ int gridFromFile(FILE *gfp, const char *dname)
grid.type = GRID_CELL;
else if ( strncmp(pline, "gme", 3) == 0 )
grid.type = GRID_GME;
else if ( strncmp(pline, "lambert", 7) == 0 )
grid.type = GRID_LCC;
else if ( strncmp(pline, "lcc2", 4) == 0 )
grid.type = GRID_LCC2;
else if ( strncmp(pline, "lcc", 3) == 0 )
grid.type = GRID_LCC;
else if ( strncmp(pline, "lambert", 7) == 0 )
grid.type = GRID_LCC;
else if ( strncmp(pline, "sinusoidal", 10) == 0 )
grid.type = GRID_SINUSOIDAL;
else if ( strncmp(pline, "laea", 4) == 0 )
......@@ -736,6 +797,30 @@ int gridFromFile(FILE *gfp, const char *dname)
else
Warning(func, "Invalid projection : %s", pline);
}
else if ( strncmp(pline, "a", 1) == 0 )
{
grid.a = atof(skipSeparator(pline + 1));
}
else if ( strncmp(pline, "lon_0", 5) == 0 )
{
grid.lon_0 = atof(skipSeparator(pline + 5));
grid.def_lon_0 = TRUE;
}
else if ( strncmp(pline, "lat_0", 5) == 0 )
{
grid.lat_0 = atof(skipSeparator(pline + 5));
grid.def_lat_0 = TRUE;
}
else if ( strncmp(pline, "lat_1", 5) == 0 )
{
grid.lat_1 = atof(skipSeparator(pline + 5));
grid.def_lat_1 = TRUE;
}
else if ( strncmp(pline, "lat_2", 5) == 0 )
{
grid.lat_2 = atof(skipSeparator(pline + 5));
grid.def_lat_2 = TRUE;
}
else if ( strncmp(pline, "xnpole", 6) == 0 )
{
grid.xpole = atof(skipSeparator(pline + 6));
......
......@@ -11,7 +11,7 @@ typedef struct {
double xlast, ylast;
double xinc, yinc;
double xpole, ypole, angle; /* rotated north pole */
double originLon; /* lambert */
double originLon; /* lambert */
double originLat;
double lonParY;
double lat1;
......@@ -23,6 +23,15 @@ typedef struct {
int def_lonParY;
int def_lat1;
int def_lat2;
double a;
double lon_0;
double lat_0;
double lat_1;
double lat_2;
int def_lon_0;
int def_lat_0;
int def_lat_1;
int def_lat_2;
int prec;
int isRotated; /* TRUE for rotated grids */
int type;
......
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