Commit 3d45637f authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

No commit message

No commit message
parent d8b6293a
2009-10-08 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
2009-10-18 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* IEG format: bug fix for lonlat grids
* using CDI library version 1.4.0.1
* IEG format: bug fix for identification of lonlat grids
* GRIB format: bug fix for decoding of missing values (scalar version only)
* Version 1.4.0.1 released
2009-10-07 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
......
#! /bin/sh
# Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.61 for cdo 1.4.1.
# Generated by GNU Autoconf 2.61 for cdo 1.4.0.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.4.1'
PACKAGE_STRING='cdo 1.4.1'
PACKAGE_VERSION='1.4.0.1'
PACKAGE_STRING='cdo 1.4.0.1'
PACKAGE_BUGREPORT='Uwe.Schulzweida@zmaw.de'
# Factoring default headers for most tests.
......@@ -1221,7 +1221,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.4.1 to adapt to many kinds of systems.
\`configure' configures cdo 1.4.0.1 to adapt to many kinds of systems.
Usage: $0 [OPTION]... [VAR=VALUE]...
......@@ -1292,7 +1292,7 @@ fi
if test -n "$ac_init_help"; then
case $ac_init_help in
short | recursive ) echo "Configuration of cdo 1.4.1:";;
short | recursive ) echo "Configuration of cdo 1.4.0.1:";;
esac
cat <<\_ACEOF
......@@ -1412,7 +1412,7 @@ fi
test -n "$ac_init_help" && exit $ac_status
if $ac_init_version; then
cat <<\_ACEOF
cdo configure 1.4.1
cdo configure 1.4.0.1
generated by GNU Autoconf 2.61
Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
......@@ -1426,7 +1426,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.4.1, which was
It was created by cdo $as_me 1.4.0.1, which was
generated by GNU Autoconf 2.61. Invocation command line was
$ $0 $@
......@@ -2252,7 +2252,7 @@ fi
# Define the identity of the package.
PACKAGE='cdo'
VERSION='1.4.1'
VERSION='1.4.0.1'
cat >>confdefs.h <<_ACEOF
......@@ -8006,7 +8006,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.4.1, which was
This file was extended by cdo $as_me 1.4.0.1, which was
generated by GNU Autoconf 2.61. Invocation command line was
CONFIG_FILES = $CONFIG_FILES
......@@ -8059,7 +8059,7 @@ Report bugs to <bug-autoconf@gnu.org>."
_ACEOF
cat >>$CONFIG_STATUS <<_ACEOF
ac_cs_version="\\
cdo config.status 1.4.1
cdo config.status 1.4.0.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.4.1, Uwe.Schulzweida@zmaw.de)
AC_INIT(cdo, 1.4.0.1, Uwe.Schulzweida@zmaw.de)
CONFIG_ABORT=yes
......
#include <ctype.h>
#include "cdi.h"
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#define NLON 1440
#define NLAT 720
#define MAX_VARS 6
static
void init_amsr_day(int vlistID, int gridID, int zaxisID, int nvars)
{
/*
Version-5 RSS AMSR-E or AMSR-J daily files
filename with path in form satname_yyyymmdd_v5.gz
where satname = name of satellite (amsre or amsr)
yyyy = year
mm = month
dd = day of month
1:time time of measurement in fractional hours GMT
2:sst sea surface temperature in deg Celcius
3:wind 10m surface wind in meters/second
4:vapor columnar water vapor in millimeters
5:cloud cloud liquid water in millimeters
6:rain rain rate in millimeters/hour
*/
char *name[] = {"hours", "sst", "wind", "vapor", "cloud", "rain"};
char *units[] = {"h", "deg Celcius", "m/s", "mm", "mm", "mm/h"};
double xscale[] = {0.1, 0.15, 0.2, 0.3, 0.01, 0.1};
double xminval[] = {0., -3., 0., 0., 0., 0.};
int i, varID;
for ( i = 0; i < nvars; ++i )
{
varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARIABLE);
vlistDefVarName(vlistID, varID, name[i]);
vlistDefVarUnits(vlistID, varID, units[i]);
vlistDefVarDatatype(vlistID, varID, DATATYPE_INT16);
vlistDefVarMissval(vlistID, varID, 254);
vlistDefVarScalefactor(vlistID, varID, xscale[i]);
vlistDefVarAddoffset(vlistID, varID, xminval[i]);
}
}
static
void init_amsr_averaged(int vlistID, int gridID, int zaxisID, int nvars)
void init_vars(int vlistID, int gridID, int zaxisID, int nvars)
{
/*
Version-5 AMSR-E or AMSR-J time-averaged files including:
3-day (average of 3 days ending on file date)
weekly (average of 7 days ending on Saturday of file date)
monthly (average of all days in month)
filename
format of file names are:
3-day satname_yyyymmddv5_d3d.gz
weekly satname_yyyymmddv5.gz
monthly satname_yyyymmv5.gz
where satname =name of satellite (amsre or amsr)
yyyy =year
mm =month
dd =day of month
1:sst sea surface temperature in deg Celcius
2:wind 10m surface wind in meters/second
3:vapor columnar water vapor in millimeters
4:cloud cloud liquid water in millimeters
5:rain rain rate in millimeters/hour
*/
char *name[] = {"sst", "wind", "vapor", "cloud", "rain"};
char *units[] = {"deg Celcius", "m/s", "mm", "mm", "mm/h"};
double xscale[] = {0.15, 0.2, 0.3, 0.01, 0.1};
double xminval[] = {-3., 0., 0., 0., 0.};
int code[] = {11, 17, 33, 34, 1, 2/*, 3*/};
char *name[] = {"temp", "unknown", "u", "v", "height", "pressure" /*, "station"*/};
char *units[] = {"Celsius", "", "m/s", "m/s", "m", "hPa" /*, ""*/};
int i, varID;
for ( i = 0; i < nvars; ++i )
{
/* varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_CONSTANT); */
varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARIABLE);
vlistDefVarCode(vlistID, varID, code[i]);
vlistDefVarName(vlistID, varID, name[i]);
vlistDefVarUnits(vlistID, varID, units[i]);
vlistDefVarDatatype(vlistID, varID, DATATYPE_INT16);
vlistDefVarMissval(vlistID, varID, 254);
vlistDefVarScalefactor(vlistID, varID, xscale[i]);
vlistDefVarAddoffset(vlistID, varID, xminval[i]);
vlistDefVarDatatype(vlistID, varID, DATATYPE_FLT32);
}
}
static
void read_amsr(FILE *fp, int vlistID, int nvars, double *data[], int *nmiss)
void init_data(int vlistID, int nvars, double *data[])
{
static char func[] = "read_amsr_averaged";
static char func[] = "initdata";
int varID, i, gridsize;
unsigned char *amsr_data = NULL;
double xminval, xscale, missval;
size_t size;
double missval;
for ( varID = 0; varID < nvars; ++varID )
{
gridsize = gridInqSize(vlistInqVarGrid(vlistID, varID));
amsr_data = (unsigned char *) realloc(amsr_data, gridsize);
size = fread(amsr_data, 1, gridsize, fp);
if ( (int)size != gridsize ) cdoAbort("Read error!");
missval = vlistInqVarMissval(vlistID, varID);
xminval = vlistInqVarAddoffset(vlistID, varID);
xscale = vlistInqVarScalefactor(vlistID, varID);
missval = vlistInqVarMissval(vlistID, varID);
nmiss[varID] = 0;
for ( i = 0; i < gridsize; ++i )
{
if ( amsr_data[i] <= 250 )
{
data[varID][i] = amsr_data[i]*xscale + xminval;
}
else
{
data[varID][i] = missval;
nmiss[varID]++;
}
data[varID][i] = missval;
}
}
free(amsr_data);
}
static
void write_data(int streamID, int nvars, double *data[], int *nmiss)
void write_data(int streamID, int vlistID, int nvars, double *data[])
{
int i;
int varID;
int nmiss;
int gridsize;
double missval;
for ( varID = 0; varID < nvars; ++varID )
{
gridsize = gridInqSize(vlistInqVarGrid(vlistID, varID));
missval = vlistInqVarMissval(vlistID, varID);
streamDefRecord(streamID, varID, 0);
streamWriteRecord(streamID, data[varID], nmiss[varID]);
nmiss = 0;
for ( i = 0; i < gridsize; ++i )
if ( DBL_IS_EQUAL(data[varID][i], missval) ) nmiss++;
streamWriteRecord(streamID, data[varID], nmiss);
}
}
......@@ -167,25 +94,55 @@ int getDate(const char *name)
void *Importobs(void *argument)
{
static char func[] = "Importobs";
int operatorID;
char line[MAX_LINE_LEN];
int streamID;
int tsID;
int gridID, zaxisID, taxisID, vlistID;
int gridsize;
int i;
int nvars;
int varID, gridID, zaxisID, taxisID, vlistID;
int i, j;
int nvars = MAX_VARS;
int vdate = 0, vtime = 0;
double xvals[NLON], yvals[NLAT];
int vdate0 = 0, vtime0 = 0;
int gridsize, xsize, ysize;
double *xvals = NULL, *yvals = NULL;
double *data[MAX_VARS];
int nmiss[MAX_VARS];
int nmiss;
FILE *fp;
size_t fsize;
char dummy[32], station[32], datetime[32];
float lat, lon, height1, pressure, height2, value;
double missval;
double latmin = 90, latmax = -90, lonmin = 360, lonmax = -360;
int code;
int index;
double dx, dy;
char *pstation;
cdoInitialize(argument);
cdoOperatorAdd("import_obs", 0, 0, "grid description file or name");
operatorID = cdoOperatorID();
operatorInputArg(cdoOperatorEnter(operatorID));
gridID = cdoDefineGrid(operatorArgv()[0]);
if ( gridInqType(gridID) != GRID_LONLAT )
cdoAbort("%s grid unsupported!", gridNamePtr(gridInqType(gridID)));
gridsize = gridInqSize(gridID);
xsize = gridInqXsize(gridID);
ysize = gridInqYsize(gridID);
// printf("gridsize=%d, xsize=%d, ysize=%d\n", gridsize, xsize, ysize);
xvals = (double *) malloc(gridsize*sizeof(double));
yvals = (double *) malloc(gridsize*sizeof(double));
gridInqXvals(gridID, xvals);
gridInqYvals(gridID, yvals);
fp = fopen(cdoStreamName(0), "r");
if ( fp == NULL ) { perror(cdoStreamName(0)); exit(EXIT_FAILURE); }
......@@ -195,19 +152,6 @@ void *Importobs(void *argument)
streamID = streamOpenWrite(cdoStreamName(1), cdoFiletype());
if ( streamID < 0 ) cdiError(streamID, "Open failed on %s", cdoStreamName(1));
/*
Longitude is 0.25*xdim-0.125 degrees east
Latitude is 0.25*ydim-90.125
*/
gridsize = NLON*NLAT;
gridID = gridCreate(GRID_LONLAT, gridsize);
gridDefXsize(gridID, NLON);
gridDefYsize(gridID, NLAT);
for ( i = 0; i < NLON; ++i ) xvals[i] = 0.25*(i+1) - 0.125;
for ( i = 0; i < NLAT; ++i ) yvals[i] = 0.25*(i+1) - 90.125;
gridDefXvals(gridID, xvals);
gridDefYvals(gridID, yvals);
zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
taxisID = taxisCreate(TAXIS_ABSOLUTE);
......@@ -216,38 +160,83 @@ void *Importobs(void *argument)
vlistDefTaxis(vlistID, taxisID);
{
nvars = 5;
for ( i = 0; i < nvars; ++i ) data[i] = (double *) malloc(gridsize*sizeof(double));
init_amsr_averaged(vlistID, gridID, zaxisID, nvars);
init_vars(vlistID, gridID, zaxisID, nvars);
/* vlistDefNtsteps(vlistID, 0);*/
streamDefVlist(streamID, vlistID);
taxisDefVdate(taxisID, vdate);
taxisDefVtime(taxisID, vtime);
tsID = 0;
streamDefTimestep(streamID, tsID);
processDefTimesteps(streamID);
// read_amsr(fp, vlistID, nvars, data, nmiss);
vdate = 0;
vtime = 0;
vdate0 = 0;
vtime0 = 0;
//ntime = 0;
tsID = 0;
while ( readline(fp, line, MAX_LINE_LEN) )
{
sscanf(line, "%s %s %s %g %g %g %d %g %g %g",
dummy, station, datetime, &lat, &lon, &height1, &code, &pressure, &height2, &value);
sscanf(datetime, "%d_%d", &vdate, &vtime);
if ( vdate != vdate0 || vtime != vtime0 )
{
if ( tsID > 0 ) write_data(streamID, vlistID, nvars, data);
vdate0 = vdate;
vtime0 = vtime;
/*
printf("%s %d %d %g %g %g %d %g %g %g\n",
station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
*/
taxisDefVdate(taxisID, vdate);
taxisDefVtime(taxisID, vtime);
streamDefTimestep(streamID, tsID);
init_data(vlistID, nvars, data);
tsID++;
}
if ( lon < lonmin ) lonmin = lon;
if ( lon > lonmax ) lonmax = lon;
if ( lat < latmin ) latmin = lat;
if ( lat > latmax ) latmax = lat;
dy = yvals[1] - yvals[0];
for ( j = 0; j < ysize; ++j )
if ( lat >= (yvals[j]-dy/2) && lat < (yvals[j]+dy/2) ) break;
dx = xvals[1] - xvals[0];
for ( i = 0; i < xsize; ++i )
if ( lon >= (xvals[i]-dx/2) && lon < (xvals[i]+dx/2) ) break;
index = -1;
if ( code == 11 ) index = 0;
if ( code == 17 ) index = 1;
if ( code == 33 ) index = 2;
if ( code == 34 ) index = 3;
//printf("%d %d %d %g %g %g %g\n", i, j, index, dx, dy, lon, lat);
if ( i < xsize && j < ysize && index >= 0 )
{
pstation = station;
while (isalpha(*pstation)) *pstation++;
// printf("station %s %d\n", pstation, atoi(pstation));
data[index][j*xsize+i] = value;
data[ 4][j*xsize+i] = height1;
data[ 5][j*xsize+i] = pressure;
// data[ 6][j*xsize+i] = atoi(pstation);
}
/*
printf("%s %d %d %g %g %g %d %g %g %g\n",
station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
*/
}
write_data(streamID, nvars, data, nmiss);
write_data(streamID, vlistID, nvars, data);
for ( i = 0; i < nvars; ++i ) free(data[i]);
}
printf("lonmin=%g, lonmax=%g, latmin=%g, latmax=%g\n", lonmin, lonmax, latmin, latmax);
processDefVarNum(vlistNvars(vlistID), streamID);
......@@ -260,6 +249,9 @@ void *Importobs(void *argument)
zaxisDestroy(zaxisID);
taxisDestroy(taxisID);
free(xvals);
free(yvals);
cdoFinish();
return (0);
......
......@@ -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
......@@ -185,6 +185,7 @@ static int gengrid(int gridID1, int lhalo, int rhalo)
int gridtype, gridID2;
int nlon1, nlat1;
int nlon2, nlat2;
int nmin, nmax;
int i;
int prec;
int ilat, ilon;
......@@ -203,6 +204,14 @@ static int gengrid(int gridID1, int lhalo, int rhalo)
nlon2 = nlon1 + lhalo + rhalo;
nlat2 = nlat1;
nmin = 0;
nmax = nlon1;
if ( lhalo < 0 ) nmin = -lhalo;
if ( rhalo < 0 ) nmax += rhalo;
/*
printf("nlon1=%d, nlon2=%d, lhalo=%d, rhalo=%d nmin=%d, nmax=%d\n",
nlon1, nlon2, lhalo, rhalo, nmin, nmax);
*/
gridtype = gridInqType(gridID1);
prec = gridInqPrec(gridID1);
......@@ -259,7 +268,7 @@ static int gengrid(int gridID1, int lhalo, int rhalo)
*pyvals2++ = yvals1[ilat*nlon1 + ilon];
}
for ( ilon = 0; ilon < nlon1; ilon++ )
for ( ilon = nmin; ilon < nmax; ilon++ )
{
*pxvals2++ = xvals1[ilat*nlon1 + ilon];
*pyvals2++ = yvals1[ilat*nlon1 + ilon];
......@@ -275,7 +284,7 @@ static int gengrid(int gridID1, int lhalo, int rhalo)
else
{
for ( i = nlon1-lhalo; i < nlon1; i++ ) *pxvals2++ = xvals1[i] - 360;
for ( i = 0; i < nlon1; i++ ) *pxvals2++ = xvals1[i];
for ( i = nmin; i < nmax; i++ ) *pxvals2++ = xvals1[i];
for ( i = 0; i < rhalo; i++ ) *pxvals2++ = xvals1[i] + 360;
for ( i = 0; i < nlat1; i++ ) yvals2[i] = yvals1[i];
......@@ -327,7 +336,7 @@ static int gengrid(int gridID1, int lhalo, int rhalo)
*pybounds2++ = ybounds1[4*ilat*nlon1 + ilon];
}
for ( ilon = 0; ilon < 4*nlon1; ilon++ )
for ( ilon = nmin; ilon < 4*nmax; ilon++ )
{
*pxbounds2++ = xbounds1[4*ilat*nlon1 + ilon];
*pybounds2++ = ybounds1[4*ilat*nlon1 + ilon];
......@@ -344,7 +353,7 @@ static int gengrid(int gridID1, int lhalo, int rhalo)
{
gridDefNvertex(gridID2, 2);
for ( i = 2*(nlon1-lhalo); i < 2*nlon1; i++ ) *pxbounds2++ = xbounds1[i] - 360;
for ( i = 0; i < 2*nlon1; i++ ) *pxbounds2++ = xbounds1[i];
for ( i = 2*nmin; i < 2*nmax; i++ ) *pxbounds2++ = xbounds1[i];
for ( i = 0; i < 2*rhalo; i++ ) *pxbounds2++ = xbounds1[i] + 360;
for ( i = 0; i < 2*nlat2; i++ ) ybounds2[i] = ybounds1[i];
......@@ -375,16 +384,28 @@ static int genindexgrid(int gridID1, int *lhalo, int *rhalo)
nlon1 = gridInqXsize(gridID1);
if ( *lhalo < 0 || *lhalo > nlon1 )
if ( *lhalo > nlon1 )
{
*lhalo = nlon1;
cdoWarning("left halo out of range. Set to %d.", *lhalo);
}
if ( *lhalo < 0 && -(*lhalo) > nlon1/2 )
{
cdoWarning("left halo out of range. Set to 0.");
*lhalo = 0;
*lhalo = -nlon1/2;
cdoWarning("left halo out of range. Set to %d.", *lhalo);
}
if ( *rhalo < 0 || *rhalo > nlon1 )
if ( *rhalo > nlon1 )
{
cdoWarning("right halo out of range. Set to 0.");
*rhalo = 0;
*rhalo = nlon1;
cdoWarning("right halo out of range. Set to %d.", rhalo);
}
if ( *rhalo < 0 && -(*rhalo) > nlon1/2 )
{
*rhalo = -nlon1/2;
cdoWarning("right halo out of range. Set to %d.", rhalo);
}
gridID2 = gengrid(gridID1, *lhalo, *rhalo);
......@@ -397,16 +418,22 @@ static void halo(double *array1, int gridID1, double *array2, int lhalo, int rha
{
int nlon1, nlat;
int ilat, ilon;
int nmin, nmax;
nlon1 = gridInqXsize(gridID1);
nlat = gridInqYsize(gridID1);
nlat = gridInqYsize(gridID1);
nmin = 0;
nmax = nlon1;
if ( lhalo < 0 ) nmin = -lhalo;
if ( rhalo < 0 ) nmax += rhalo;
for ( ilat = 0; ilat < nlat; ilat++ )
{
for ( ilon = nlon1-lhalo; ilon < nlon1; ilon++ )
*array2++ = array1[ilat*nlon1 + ilon];
for ( ilon = 0; ilon < nlon1; ilon++ )
for ( ilon = nmin; ilon < nmax; ilon++ )
*array2++ = array1[ilat*nlon1 + ilon];
for ( ilon = 0; ilon < rhalo; ilon++ )
......
......@@ -706,8 +706,8 @@ int pstreamOpenWrite(const char *argument, int filetype)
streamDefZtype(fileID, cdoZtype);
streamDefZlevel(fileID, cdoZlevel);
if ( cdoZtype == COMPRESS_SZIP && filetype != FILETYPE_GRB )
cdoWarning("SZIP compression not available for non GRIB data!");
if ( cdoZtype == COMPRESS_SZIP && (filetype != FILETYPE_GRB && filetype != FILETYPE_NC4) )
cdoWarning("SZIP compression not available for non GRIB/netCDF4 data!");
if ( cdoZtype == COMPRESS_ZIP && filetype != FILETYPE_NC4 )
cdoWarning("Deflate compression not available for non netCDF4 data!");
......
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