From 3d45637fe72956912acea62fb72cde3919e504a8 Mon Sep 17 00:00:00 2001 From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de> Date: Sun, 18 Oct 2009 19:35:32 +0000 Subject: [PATCH] --- ChangeLog | 7 +- configure | 20 ++-- configure.ac | 2 +- src/Importobs.c | 258 +++++++++++++++++++++++------------------------- src/Sethalo.c | 53 +++++++--- src/pstream.c | 4 +- 6 files changed, 183 insertions(+), 161 deletions(-) diff --git a/ChangeLog b/ChangeLog index 78a5493c3..15496082e 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,6 +1,9 @@ -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> diff --git a/configure b/configure index d096fbc13..c2e6cd722 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /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'`\\" diff --git a/configure.ac b/configure.ac index 168c8b203..26e48dc46 100644 --- a/configure.ac +++ b/configure.ac @@ -1,6 +1,6 @@ # 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 diff --git a/src/Importobs.c b/src/Importobs.c index 30eaddcf8..c79886c0a 100644 --- a/src/Importobs.c +++ b/src/Importobs.c @@ -1,147 +1,74 @@ +#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); diff --git a/src/Sethalo.c b/src/Sethalo.c index 57826b33f..ed0cffc5c 100644 --- a/src/Sethalo.c +++ b/src/Sethalo.c @@ -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++ ) diff --git a/src/pstream.c b/src/pstream.c index 24ad449c7..7817140ba 100644 --- a/src/pstream.c +++ b/src/pstream.c @@ -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!"); -- GitLab