/* v5d.c */
/* Vis5D version 5.2 */
/*
Vis5D system for visualizing five dimensional gridded data sets
Copyright (C) 1990 - 2000 Bill Hibbard, Johan Kellum, Brian Paul,
Dave Santek, and Andre Battaiola.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 1, or (at your option)
any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
/* this should be updated when the file version changes */
#define FILE_VERSION "4.3"
/*
* New grid file format for VIS-5D:
*
* The header is a list of tagged items. Each item has 3 parts:
* 1. A tag which is a 4-byte integer identifying the type of item.
* 2. A 4-byte integer indicating how many bytes of data follow.
* 3. The binary data.
*
* If we need to add new information to a file header we just create a
* new tag and add the code to read/write the information.
*
* If we're reading a header and find an unknown tag, we can use the
* length field to skip ahead to the next tag. Therefore, the file
* format is forward (and backward) compatible.
*
* Grid data is stored as either:
* 1-byte unsigned integers (255=missing)
* 2-byte unsigned integers (65535=missing)
* 4-byte IEEE floats ( >1.0e30 = missing)
*
* All numeric values are stored in big endian order. All floating point
* values are in IEEE format.
*/
/*
* Updates:
*
* April 13, 1995, brianp
* finished Cray support for 2-byte and 4-byte compress modes
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "binio.h"
#include "v5d.h"
#include "vis5d.h"
#ifndef SEEK_SET
# define SEEK_SET 0
#endif
#ifndef SEEK_CUR
# define SEEK_CUR 1
#endif
#ifndef SEEK_END
# define SEEK_END 2
#endif
/*
* Currently defined tags:
* Note: the notation a[i] doesn't mean a is an array of i elements,
* rather it just refers to the ith element of a[].
*
* Tags marked as PHASED OUT should be readable but are no longer written.
* Old tag numbers can't be reused!
*
*/
/* TAG NAME VALUE DATA (comments) */
/*----------------------------------------------------------------------*/
#define TAG_ID 0x5635440a /* hex encoding of "V5D\n" */
/* general stuff 1000+ */
#define TAG_VERSION 1000 /* char*10 FileVersion */
#define TAG_NUMTIMES 1001 /* int*4 NumTimes */
#define TAG_NUMVARS 1002 /* int*4 NumVars */
#define TAG_VARNAME 1003 /* int*4 var; char*10 VarName[var] */
#define TAG_NR 1004 /* int*4 Nr */
#define TAG_NC 1005 /* int*4 Nc */
#define TAG_NL 1006 /* int*4 Nl (Nl for all vars) */
#define TAG_NL_VAR 1007 /* int*4 var; int*4 Nl[var] */
#define TAG_LOWLEV_VAR 1008 /* int*4 var; int*4 LowLev[var] */
#define TAG_TIME 1010 /* int*4 t; int*4 TimeStamp[t] */
#define TAG_DATE 1011 /* int*4 t; int*4 DateStamp[t] */
#define TAG_MINVAL 1012 /* int*4 var; real*4 MinVal[var] */
#define TAG_MAXVAL 1013 /* int*4 var; real*4 MaxVal[var] */
#define TAG_COMPRESS 1014 /* int*4 CompressMode; (#bytes/grid)*/
#define TAG_UNITS 1015 /* int *4 var; char*20 Units[var] */
/* vertical coordinate system 2000+ */
#define TAG_VERTICAL_SYSTEM 2000 /* int*4 VerticalSystem */
#define TAG_VERT_ARGS 2100 /* int*4 n; real*4 VertArgs[0..n-1]*/
#define TAG_BOTTOMBOUND 2001 /* real*4 BottomBound (PHASED OUT) */
#define TAG_LEVINC 2002 /* real*4 LevInc (PHASED OUT) */
#define TAG_HEIGHT 2003 /* int*4 l; real*4 Height[l] (PHASED OUT) */
/* projection 3000+ */
#define TAG_PROJECTION 3000 /* int*4 projection: */
/* 0 = generic linear */
/* 1 = cylindrical equidistant */
/* 2 = Lambert conformal/Polar Stereo */
/* 3 = rotated equidistant */
#define TAG_PROJ_ARGS 3100 /* int *4 n; real*4 ProjArgs[0..n-1] */
#define TAG_NORTHBOUND 3001 /* real*4 NorthBound (PHASED OUT) */
#define TAG_WESTBOUND 3002 /* real*4 WestBound (PHASED OUT) */
#define TAG_ROWINC 3003 /* real*4 RowInc (PHASED OUT) */
#define TAG_COLINC 3004 /* real*4 ColInc (PHASED OUT) */
#define TAG_LAT1 3005 /* real*4 Lat1 (PHASED OUT) */
#define TAG_LAT2 3006 /* real*4 Lat2 (PHASED OUT) */
#define TAG_POLE_ROW 3007 /* real*4 PoleRow (PHASED OUT) */
#define TAG_POLE_COL 3008 /* real*4 PoleCol (PHASED OUT) */
#define TAG_CENTLON 3009 /* real*4 CentralLon (PHASED OUT) */
#define TAG_CENTLAT 3010 /* real*4 CentralLat (PHASED OUT) */
#define TAG_CENTROW 3011 /* real*4 CentralRow (PHASED OUT) */
#define TAG_CENTCOL 3012 /* real*4 CentralCol (PHASED OUT) */
#define TAG_ROTATION 3013 /* real*4 Rotation (PHASED OUT) */
#define TAG_ROWINCKM 3014 /* real*4 RowIncKm (PHASED OUT) */
#define TAG_COLINCKM 3015 /* real*4 ColIncKm (PHASED OUT) */
#define TAG_END 9999
/**********************************************************************/
/***** Miscellaneous Functions *****/
/**********************************************************************/
float pressure_to_height(float pressure)
{
return (float) DEFAULT_LOG_EXP * log((double) pressure / DEFAULT_LOG_SCALE);
}
float height_to_pressure(float height)
{
return (float) DEFAULT_LOG_SCALE * exp((double) height / DEFAULT_LOG_EXP);
}
/*
* Return current file position.
* Input: f - file descriptor
*/
static off_t ltell( int f )
{
return lseek( f, 0, SEEK_CUR );
}
/*
* Copy up to maxlen characters from src to dst stopping upon whitespace
* in src. Terminate dst with null character.
* Return: length of dst.
*/
static int copy_string2( char *dst, const char *src, int maxlen )
{
int i;
for (i=0;i<maxlen;i++) dst[i] = src[i];
for (i=maxlen-1; i>=0; i--) {
if (dst[i]==' ' || i==maxlen-1) dst[i] = 0;
else break;
}
return strlen(dst);
}
/*
* Copy up to maxlen characters from src to dst stopping upon whitespace
* in src. Terminate dst with null character.
* Return: length of dst.
*/
static int copy_string( char *dst, const char *src, int maxlen )
{
int i;
for (i=0;i<maxlen;i++) {
if (src[i]==' ' || i==maxlen-1) {
dst[i] = 0;
break;
}
else {
dst[i] = src[i];
}
}
return i;
}
/*
* Convert a date from YYDDD format to days since Jan 1, 1900.
*/
int v5dYYDDDtoDays( int yyddd )
{
int iy, id, idays;
iy = yyddd / 1000;
id = yyddd - 1000*iy;
if (iy < 50) iy += 100; /* WLH 31 July 96 << 31 Dec 99 */
idays = 365*iy + (iy-1)/4 + id;
return idays;
}
/*
* Convert a time from HHMMSS format to seconds since midnight.
*/
int v5dHHMMSStoSeconds( int hhmmss )
{
int h, m, s;
h = hhmmss / 10000;
m = (hhmmss / 100) % 100;
s = hhmmss % 100;
return s + m*60 + h*60*60;
}
/*
* Convert a day since Jan 1, 1900 to YYDDD format.
*/
int v5dDaysToYYDDD( int days )
{
int iy, id, iyyddd;
iy = (4*days)/1461;
id = days-(365*iy+(iy-1)/4);
if (iy > 99) iy = iy - 100; /* WLH 31 July 96 << 31 Dec 99 */
/* iy = iy + 1900; is the right way to fix this, but requires
changing all places where dates are printed - procrastinate */
iyyddd = iy*1000+id;
return iyyddd;
}
/*
* Convert a time in seconds since midnight to HHMMSS format.
*/
int v5dSecondsToHHMMSS( int seconds )
{
int hh, mm, ss;
hh = seconds / (60*60);
mm = (seconds / 60) % 60;
ss = seconds % 60;
return hh*10000 + mm * 100 + ss;
}
void v5dPrintStruct( const v5dstruct *v )
{
static char day[7][10] = { "Sunday", "Monday", "Tuesday", "Wednesday",
"Thursday", "Friday", "Saturday" };
int time, var, i;
int maxnl;
maxnl = 0;
for (var=0;var<v->NumVars;var++) {
if (v->Nl[var]+v->LowLev[var]>maxnl) {
maxnl = v->Nl[var]+v->LowLev[var];
}
}
if (v->FileFormat==0) {
if (v->FileVersion[0] == 0) {
printf("File format: v5d version: (4.0 or 4.1)\n");
}
else {
printf("File format: v5d version: %s\n", v->FileVersion);
}
}
else {
printf("File format: comp5d (VIS-5D 3.3 or older)\n");
}
if (v->CompressMode==1) {
printf("Compression: 1 byte per gridpoint.\n");
}
else {
printf("Compression: %d bytes per gridpoint.\n", v->CompressMode);
}
printf("header size=%d\n", v->FirstGridPos);
printf("sizeof(v5dstruct)=%d\n", sizeof(v5dstruct) );
printf("\n");
printf("NumVars = %d\n", v->NumVars );
printf("Var Name Units Rows Cols Levels LowLev MinVal MaxVal\n");
for (var=0;var<v->NumVars;var++) {
printf("%3d %-10s %-10s %3d %3d %3d %3d",
var+1, v->VarName[var], v->Units[var],
v->Nr, v->Nc, v->Nl[var], v->LowLev[var] );
if (v->MinVal[var] > v->MaxVal[var]) {
printf(" MISSING MISSING\n");
}
else {
printf(" %-12g %-12g\n", v->MinVal[var], v->MaxVal[var] );
}
}
printf("\n");
printf("NumTimes = %d\n", v->NumTimes );
printf("Step Date(YYDDD) Time(HH:MM:SS) Day\n");
for (time=0;time<v->NumTimes;time++) {
int i = v->TimeStamp[time];
printf("%3d %05d %5d:%02d:%02d %s\n",
time+1,
v->DateStamp[time],
i/10000, (i/100)%100, i%100,
day[ v5dYYDDDtoDays(v->DateStamp[time]) % 7 ]);
}
printf("\n");
switch (v->VerticalSystem) {
case 0:
printf("Generic linear vertical coordinate system:\n");
printf("\tBottom Bound: %f\n", v->VertArgs[0] );
printf("\tIncrement between levels: %f\n", v->VertArgs[1] );
break;
case 1:
printf("Equally spaced levels in km:\n");
printf("\tBottom Bound: %f\n", v->VertArgs[0] );
printf("\tIncrement: %f\n", v->VertArgs[1] );
break;
case 2:
printf("Unequally spaced levels in km:\n");
printf("Level\tHeight(km)\n");
for (i=0;i<maxnl;i++) {
printf("%3d %10.3f\n", i+1, v->VertArgs[i] );
}
break;
case 3:
printf("Unequally spaced levels in mb:\n");
printf("Level\tPressure(mb)\n");
for (i=0;i<maxnl;i++) {
printf("%3d %10.3f\n", i+1, height_to_pressure(v->VertArgs[i]) );
}
break;
default:
printf("Bad VerticalSystem value: %d\n", v->VerticalSystem );
}
printf("\n");
switch (v->Projection) {
case 0:
printf("Generic linear projection:\n");
printf("\tNorth Boundary: %f\n", v->ProjArgs[0] );
printf("\tWest Boundary: %f\n", v->ProjArgs[1] );
printf("\tRow Increment: %f\n", v->ProjArgs[2] );
printf("\tColumn Increment: %f\n", v->ProjArgs[3] );
break;
case 1:
printf("Cylindrical Equidistant projection:\n");
printf("\tNorth Boundary: %f degrees\n", v->ProjArgs[0] );
printf("\tWest Boundary: %f degrees\n", v->ProjArgs[1] );
printf("\tRow Increment: %f degrees\n", v->ProjArgs[2] );
printf("\tColumn Increment: %f degrees\n", v->ProjArgs[3] );
/*
printf("\tSouth Boundary: %f degrees\n",
v->NorthBound - v->RowInc * (v->Nr-1) );
printf("\tEast Boundary: %f degrees\n",
v->WestBound - v->ColInc * (v->Nc-1) );
*/
break;
case 2:
printf("Lambert Conformal projection:\n");
printf("\tStandard Latitude 1: %f\n", v->ProjArgs[0] );
printf("\tStandard Latitude 2: %f\n", v->ProjArgs[1] );
printf("\tNorth/South Pole Row: %f\n", v->ProjArgs[2] );
printf("\tNorth/South Pole Column: %f\n", v->ProjArgs[3] );
printf("\tCentral Longitude: %f\n", v->ProjArgs[4] );
printf("\tColumn Increment: %f km\n", v->ProjArgs[5] );
break;
case 3:
printf("Stereographic:\n");
printf("\tCenter Latitude: %f\n", v->ProjArgs[0] );
printf("\tCenter Longitude: %f\n", v->ProjArgs[1] );
printf("\tCenter Row: %f\n", v->ProjArgs[2] );
printf("\tCenter Column: %f\n", v->ProjArgs[3] );
printf("\tColumn Spacing: %f\n", v->ProjArgs[4] );
break;
case 4:
/* WLH 4-21-95 */
printf("Rotated equidistant projection:\n");
printf("\tLatitude of grid(0,0): %f\n", v->ProjArgs[0] );
printf("\tLongitude of grid(0,0): %f\n", v->ProjArgs[1] );
printf("\tRow Increment: %f degress\n", v->ProjArgs[2] );
printf("\tColumn Increment: %f degrees\n", v->ProjArgs[3] );
printf("\tCenter Latitude: %f\n", v->ProjArgs[4] );
printf("\tCenter Longitude: %f\n", v->ProjArgs[5] );
printf("\tRotation: %f degrees\n", v->ProjArgs[6] );
break;
case 5:
printf("Mercator:\n");
printf("\tCenter Latitude: %f\n", v->ProjArgs[0] );
printf("\tCenter Longitude: %f\n", v->ProjArgs[1] );
printf("\tRow Increment: %f Kilometers\n", v->ProjArgs[2] );
printf("\tCol Increment: %f Kilometers\n", v->ProjArgs[3] );
break;
default:
printf("Bad projection number: %d\n", v->Projection );
}
}
/*
* Compute the location of a compressed grid within a file.
* Input: v - pointer to v5dstruct describing the file header.
* time, var - which timestep and variable.
* Return: file offset in bytes
*/
static int grid_position( const v5dstruct *v, int time, int var )
{
int pos, i;
assert( time >= 0 );
assert( var >= 0 );
assert( time < v->NumTimes );
assert( var < v->NumVars );
pos = v->FirstGridPos + time * v->SumGridSizes;
for (i=0;i<var;i++) {
pos += v->GridSize[i];
}
return pos;
}
/*
* Compute the ga and gb (de)compression values for a grid.
* Input: nr, nc, nl - size of grid
* data - the grid data
* ga, gb - arrays to store results.
* minval, maxval - pointer to floats to return min, max values
* compressmode - 1, 2 or 4 bytes per grid point
* Output: ga, gb - the (de)compression values
* minval, maxval - the min and max grid values
* Side effect: the MinVal[var] and MaxVal[var] fields in g may be
* updated with new values.
*/
static void compute_ga_gb( int nr, int nc, int nl,
const float data[], int compressmode,
float ga[], float gb[],
float *minval, float *maxval )
{
#ifdef SIMPLE_COMPRESSION
/*
* Compute ga, gb values for whole grid.
*/
int i, lev, allmissing, num;
float min, max, a, b;
min = 1.0e30;
max = -1.0e30;
num = nr * nc * nl;
allmissing = 1;
for (i=0;i<num;i++) {
if (!IS_MISSING(data[i])) {
if (data[i]<min) min = data[i];
if (data[i]>max) max = data[i];
allmissing = 0;
}
}
if (allmissing) {
a = 1.0;
b = 0.0;
}
else {
a = (max-min) / 254.0;
b = min;
}
/* return results */
for (i=0;i<nl;i++) {
ga[i] = a;
gb[i] = b;
}
*minval = min;
*maxval = max;
#else
/*
* Compress grid on level-by-level basis.
*/
# define SMALLVALUE -1.0e30
# define BIGVALUE 1.0e30
# define ABS(x) ( ((x) < 0.0) ? -(x) : (x) )
float gridmin, gridmax;
float levmin[MAXLEVELS], levmax[MAXLEVELS];
float d[MAXLEVELS], dmax;
float ival, mval;
int j, k, lev, nrnc;
nrnc = nr * nc;
/* find min and max for each layer and the whole grid */
gridmin = BIGVALUE;
gridmax = SMALLVALUE;
j = 0;
for (lev=0;lev<nl;lev++) {
float ave, var;
float min, max;
min = BIGVALUE;
max = SMALLVALUE;
ave = 0.0;
var = 0.0;
for (k=0;k<nrnc;k++) {
if (!IS_MISSING(data[j]) && data[j]<min)
min = data[j];
if (!IS_MISSING(data[j]) && data[j]>max)
max = data[j];
j++;
}
if (min<gridmin)
gridmin = min;
if (max>gridmax)
gridmax = max;
levmin[lev] = min;
levmax[lev] = max;
}
/* WLH 2-2-95 */
#ifdef KLUDGE
/* if the grid minimum is within delt of 0.0, fudge all values */
/* within delt of 0.0 to delt, and recalculate mins and maxes */
{
float delt;
int nrncnl = nrnc * nl;
delt = (gridmax - gridmin)/100000.0;
if ( ABS(gridmin) < delt && gridmin!=0.0 && compressmode != 4 ) {
float min, max;
for (j=0; j<nrncnl; j++) {
if (!IS_MISSING(data[j]) && data[j]<delt)
data[j] = delt;
}
/* re-calculate min and max for each layer and the whole grid */
gridmin = delt;
for (lev=0;lev<nl;lev++) {
if (ABS(levmin[lev]) < delt)
levmin[lev] = delt;
if (ABS(levmax[lev]) < delt)
levmax[lev] = delt;
}
}
}
#endif
/* find d[lev] and dmax = MAX( d[0], d[1], ... d[nl-1] ) */
dmax = 0.0;
for (lev=0;lev<nl;lev++) {
if (levmin[lev]>=BIGVALUE && levmax[lev]<=SMALLVALUE) {
/* all values in the layer are MISSING */
d[lev] = 0.0;
}
else {
d[lev] = levmax[lev]-levmin[lev];
}
if (d[lev]>dmax)
dmax = d[lev];
}
/*** Compute ga (scale) and gb (bias) for each grid level */
if (dmax==0.0) {
/*** Special cases ***/
if (gridmin==gridmax) {
/*** whole grid is of same value ***/
for (lev=0; lev<nl; lev++) {
ga[lev] = gridmin;
gb[lev] = 0.0;
}
}
else {
/*** every layer is of a single value ***/
for (lev=0; lev<nl; lev++) {
ga[lev] = levmin[lev];
gb[lev] = 0.0;
}
}
}
else {
/*** Normal cases ***/
if (compressmode == 1) {
#define ORIGINAL
#ifdef ORIGINAL
ival = dmax / 254.0;
mval = gridmin;
for (lev=0; lev<nl; lev++) {
ga[lev] = ival;
/* WLH 2 Aug 2000 */
if (levmin[lev] > levmax[lev]) {
gb[lev] = 0.0;
}
else {
gb[lev] = mval + ival * (int) ( (levmin[lev]-mval) / ival );
}
}
#else
for (lev=0; lev<nl; lev++) {
if (d[lev]==0.0) {
ival = 1.0;
}
else {
ival = d[lev] / 254.0;
}
ga[lev] = ival;
gb[lev] = levmin[lev];
}
#endif
}
else if (compressmode == 2) {
ival = dmax / 65534.0;
mval = gridmin;
for (lev=0; lev<nl; lev++) {
ga[lev] = ival;
/* WLH 2 Aug 2000 */
if (levmin[lev] > levmax[lev]) {
gb[lev] = 0.0;
}
else {
gb[lev] = mval + ival * (int) ( (levmin[lev]-mval) / ival );
}
}
}
else {
assert( compressmode==4 );
for (lev=0; lev<nl; lev++) {
ga[lev] = 1.0;
gb[lev] = 0.0;
}
}
}
/* update min, max values */
*minval = gridmin;
*maxval = gridmax;
#endif
}
/*
* Compress a 3-D grid from floats to 1-byte unsigned integers.
* Input: nr, nc, nl - size of grid
* compressmode - 1, 2 or 4 bytes per grid point
* data - array of [nr*nc*nl] floats
* compdata - pointer to array of [nr*nc*nl*compressmode] bytes
* to put results into.
* ga, gb - pointer to arrays to put ga and gb decompression values
* minval, maxval - pointers to float to return min & max values
* Output: compdata - the compressed grid data
* ga, gb - the decompression values
* minval, maxval - the min and max grid values
*/
void v5dCompressGrid( int nr, int nc, int nl, int compressmode,
const float data[],
void *compdata, float ga[], float gb[],
float *minval, float *maxval )
{
int nrnc = nr * nc;
int nrncnl = nr * nc * nl;
V5Dubyte *compdata1 = (V5Dubyte *) compdata;
V5Dushort *compdata2 = (V5Dushort *) compdata;
/* compute ga, gb values */
compute_ga_gb( nr, nc, nl, data, compressmode, ga, gb, minval, maxval );
/* compress the data */
if (compressmode==1) {
int i, lev, p;
p = 0;
for (lev=0;lev<nl;lev++) {
float one_over_a, b;
/* WLH 5 Nov 98
b = gb[lev] - 0.0001;
*/
/* WLH 5 Nov 98 */
b = gb[lev];
/* subtract an epsilon so the int((d-b)/a) */
/* expr below doesn't get mis-truncated. */
if (ga[lev]==0.0) {
one_over_a = 1.0;
}
else {
one_over_a = 1.0 / ga[lev];
}
for (i=0;i<nrnc;i++,p++) {
if (IS_MISSING(data[p])) {
compdata1[p] = 255;
}
else {
/* MJK 1.19.99
compdata1[p] = (V5Dubyte) (int) ((data[p]-b) * one_over_a);
*/
compdata1[p] = (V5Dubyte) rint((data[p]-b) * one_over_a);
if (compdata1[p] >= 255){
compdata1[p] = (V5Dubyte) (int) (255.0 - .0001);
}
}
}
}
}
else if (compressmode == 2) {
int i, lev, p;
p = 0;
for (lev=0;lev<nl;lev++) {
float one_over_a, b;
/* WLH 5 Nov 98
b = gb[lev] - 0.0001;
*/
/* WLH 5 Nov 98 */
b = gb[lev];
if (ga[lev]==0.0) {
one_over_a = 1.0;
}
else {
one_over_a = 1.0 / ga[lev];
}
#ifdef _CRAY
/* this is tricky because sizeof(V5Dushort)==8, not 2 */
for (i=0;i<nrnc;i++,p++) {
V5Dushort compvalue;
if (IS_MISSING(data[p])) {
compvalue = 65535;
}
else {
/* MJK 3.2.99
compvalue = (V5Dushort) (int) ((data[p]-b) * one_over_a);
*/
compvalue = (V5Dushort) rint((data[p]-b) * one_over_a);
}
compdata1[p*2+0] = compvalue >> 8; /* upper byte */
compdata1[p*2+1] = compvalue & 0xffu; /* lower byte */
}
#else
for (i=0;i<nrnc;i++,p++) {
if (IS_MISSING(data[p])) {
compdata2[p] = 65535;
}
else {
compdata2[p] = (V5Dushort) rint((data[p]-b) * one_over_a);
/*
compdata2[p] = (V5Dushort) (int) ((data[p]-b) * one_over_a);
*/
/* MJK 3.24.99 I put this here so if the value is close
to the missing value and get's rounded up it won't come out
as missing data */
if (compdata2[p] == 65535){
compdata2[p] = 65534;
}
}
}
/* TODO: byte-swapping on little endian??? */
#endif
}
}
else {
/* compressmode==4 */
#ifdef _CRAY
cray_to_ieee_array( compdata, data, nrncnl );
#else
/* other machines: just copy 4-byte IEEE floats */
assert( sizeof(float)==4 );
memcpy( compdata, data, nrncnl*4 );
/* TODO: byte-swapping on little endian??? */
#endif
}
}
/*
* Decompress a 3-D grid from 1-byte integers to 4-byte floats.
* Input: nr, nc, nl - size of grid
* compdata - array of [nr*nr*nl*compressmode] bytes
* ga, gb - arrays of decompression factors
* compressmode - 1, 2 or 4 bytes per grid point
* data - address to put decompressed values
* Output: data - uncompressed floating point data values
*/
void v5dDecompressGrid( int nr, int nc, int nl, int compressmode,
void *compdata, float ga[], float gb[],
float data[] )
{
int nrnc = nr * nc;
int nrncnl = nr * nc * nl;
V5Dubyte *compdata1 = (V5Dubyte *) compdata;
V5Dushort *compdata2 = (V5Dushort *) compdata;
if (compressmode == 1) {
int p, i, lev;
p = 0;
for (lev=0;lev<nl;lev++) {
float a = ga[lev];
float b = gb[lev];
/* WLH 2-2-95 */
float d, aa;
int id;
if (a > 0.0000000001) {
d = b / a;
id = floor(d);
d = d - id;
aa = a * 0.000001;
}
else {
id = 1;
}
if (-254 <= id && id <= 0 && d < aa) {
for (i=0;i<nrnc;i++,p++) {
if (compdata1[p]==255) {
data[p] = MISSING;
}
else {
data[p] = (float) (int) compdata1[p] * a + b;
if (fabs(data[p]) < aa) data[p] = aa;
}
}
}
else {
for (i=0;i<nrnc;i++,p++) {
if (compdata1[p]==255) {
data[p] = MISSING;
}
else {
data[p] = (float) (int) compdata1[p] * a + b;
}
}
}
/* end of WLH 2-2-95 */
}
}
else if (compressmode == 2) {
int p, i, lev;
p = 0;
for (lev=0;lev<nl;lev++) {
float a = ga[lev];
float b = gb[lev];
#ifdef _CRAY
/* this is tricky because sizeof(V5Dushort)==8, not 2 */
for (i=0;i<nrnc;i++,p++) {
int compvalue;
compvalue = (compdata1[p*2] << 8) | compdata1[p*2+1];
if (compvalue==65535) {
data[p] = MISSING;
}
else {
data[p] = (float) compvalue * a + b;
}
}
#else
/* sizeof(V5Dushort)==2! */
for (i=0;i<nrnc;i++,p++) {
if (compdata2[p]==65535) {
data[p] = MISSING;
}
else {
data[p] = (float) (int) compdata2[p] * a + b;
}
}
#endif
}
}
else {
/* compressmode==4 */
#ifdef _CRAY
ieee_to_cray_array( data, compdata, nrncnl );
#else
/* other machines: just copy 4-byte IEEE floats */
assert( sizeof(float)==4 );
memcpy( data, compdata, nrncnl*4 );
#endif
}
}
/*
* Return the size (in bytes) of the 3-D grid specified by time and var.
* Input: v - pointer to v5dstruct describing the file
* time, var - which timestep and variable
* Return: number of data points.
*/
int v5dSizeofGrid( const v5dstruct *v, int time, int var )
{
return v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
}
/*
* Initialize a v5dstructure to reasonable initial values.
* Input: v - pointer to v5dstruct.
*/
void v5dInitStruct( v5dstruct *v )
{
int i;
/* set everything to zero */
memset( v, 0, sizeof(v5dstruct) );
/* special cases */
v->Projection = -1;
v->VerticalSystem = -1;
for (i=0;i<MAXVARS;i++) {
v->MinVal[i] = MISSING;
v->MaxVal[i] = -MISSING;
v->LowLev[i] = 0;
}
/* set file version */
strcpy(v->FileVersion, FILE_VERSION);
v->CompressMode = 1;
v->FileDesc = -1;
}
/*
* Return a pointer to a new, initialized v5dstruct.
*/
v5dstruct *v5dNewStruct( void )
{
v5dstruct *v;
v = (v5dstruct *) malloc( sizeof(v5dstruct) );
if (v) {
v5dInitStruct(v);
}
return v;
}
/*
* Free an initialized v5dstruct. (Todd Plessel)
*/
void v5dFreeStruct( v5dstruct* v )
{
/*assert( v5dVerifyStruct( v ) );*/
free( v );
v = 0;
}
/*
* Do some checking that the information in a v5dstruct is valid.
* Input: v - pointer to v5dstruct
* Return: 1 = g is ok, 0 = g is invalid
*/
int v5dVerifyStruct( const v5dstruct *v )
{
int var, i, invalid, maxnl;
invalid = 0;
if (!v)
return 0;
/* Number of variables */
if (v->NumVars<0) {
printf("Invalid number of variables: %d\n", v->NumVars );
invalid = 1;
}
else if (v->NumVars>MAXVARS) {
printf("Too many variables: %d (Maximum is %d)\n",
v->NumVars, MAXVARS);
invalid = 1;
}
/* Variable Names */
for (i=0;i<v->NumVars;i++) {
if (v->VarName[i][0]==0) {
printf("Missing variable name: VarName[%d]=\"\"\n", i );
invalid = 1;
}
}
/* Number of timesteps */
if (v->NumTimes<0) {
printf("Invalid number of timesteps: %d\n", v->NumTimes );
invalid = 1;
}
else if (v->NumTimes>MAXTIMES) {
printf("Too many timesteps: %d (Maximum is %d)\n",
v->NumTimes, MAXTIMES );
invalid = 1;
}
/* Make sure timestamps are increasing */
for (i=1;i<v->NumTimes;i++) {
int date0 = v5dYYDDDtoDays( v->DateStamp[i-1] );
int date1 = v5dYYDDDtoDays( v->DateStamp[i] );
int time0 = v5dHHMMSStoSeconds( v->TimeStamp[i-1] );
int time1 = v5dHHMMSStoSeconds( v->TimeStamp[i] );
if (time1<=time0 && date1<=date0) {
printf("Timestamp for step %d must be later than step %d\n", i, i-1);
invalid = 1;
}
}
/* Rows */
if (v->Nr<2) {
printf("Too few rows: %d (2 is minimum)\n", v->Nr );
invalid = 1;
}
else if (v->Nr>MAXROWS) {
printf("Too many rows: %d (%d is maximum)\n", v->Nr, MAXROWS );
invalid = 1;
}
/* Columns */
if (v->Nc<2) {
printf("Too few columns: %d (2 is minimum)\n", v->Nc );
invalid = 1;
}
else if (v->Nc>MAXCOLUMNS) {
printf("Too many columns: %d (%d is maximum)\n", v->Nc, MAXCOLUMNS );
invalid = 1;
}
/* Levels */
maxnl = 0;
for (var=0;var<v->NumVars;var++) {
if (v->LowLev[var] < 0) {
printf("Low level cannot be negative for var %s: %d\n",
v->VarName[var], v->LowLev[var] );
invalid = 1;
}
if (v->Nl[var]<1) {
printf("Too few levels for var %s: %d (1 is minimum)\n",
v->VarName[var], v->Nl[var] );
invalid = 1;
}
if (v->Nl[var]+v->LowLev[var]>MAXLEVELS) {
printf("Too many levels for var %s: %d (%d is maximum)\n",
v->VarName[var], v->Nl[var]+v->LowLev[var], MAXLEVELS );
invalid = 1;
}
if (v->Nl[var]+v->LowLev[var]>maxnl) {
maxnl = v->Nl[var]+v->LowLev[var];
}
}
if (v->CompressMode != 1 && v->CompressMode != 2 && v->CompressMode != 4) {
printf("Bad CompressMode: %d (must be 1, 2 or 4)\n", v->CompressMode );
invalid = 1;
}
switch (v->VerticalSystem) {
case 0:
case 1:
if (v->VertArgs[1]==0.0) {
printf("Vertical level increment is zero, must be non-zero\n");
invalid = 1;
}
break;
case 2:
/* Check that Height values increase upward */
for (i=1;i<maxnl;i++) {
if (v->VertArgs[i] <= v->VertArgs[i-1]) {
printf("Height[%d]=%f <= Height[%d]=%f, level heights must increase\n",
i, v->VertArgs[i], i-1, v->VertArgs[i-1] );
invalid = 1;
break;
}
}
break;
case 3:
/* Check that Pressure values decrease upward */
for (i=1;i<maxnl;i++) {
if (v->VertArgs[i] <= v->VertArgs[i-1]) {
printf("Pressure[%d]=%f >= Pressure[%d]=%f, level pressures must decrease\n",
i, height_to_pressure(v->VertArgs[i]),
i-1, height_to_pressure(v->VertArgs[i-1]) );
invalid = 1;
break;
}
}
break;
default:
printf("VerticalSystem = %d, must be in 0..3\n", v->VerticalSystem );
invalid = 1;
}
switch (v->Projection) {
case 0: /* Generic */
if (v->ProjArgs[2]==0.0) {
printf("Row Increment (ProjArgs[2]) can't be zero\n");
invalid = 1;
}
if (v->ProjArgs[3]==0.0) {
printf("Column increment (ProjArgs[3]) can't be zero\n");
invalid = 1;
}
break;
case 1: /* Cylindrical equidistant */
if (v->ProjArgs[2]<0.0) {
printf("Row Increment (ProjArgs[2]) = %g (must be >=0.0)\n",
v->ProjArgs[2] );
invalid = 1;
}
if (v->ProjArgs[3]<=0.0) {
printf("Column Increment (ProjArgs[3]) = %g (must be >=0.0)\n",
v->ProjArgs[3] );
invalid = 1;
}
break;
case 2: /* Lambert Conformal */
if (v->ProjArgs[0]<-90.0 || v->ProjArgs[0]>90.0) {
printf("Lat1 (ProjArgs[0]) out of range: %g\n", v->ProjArgs[0] );
invalid = 1;
}
if (v->ProjArgs[1]<-90.0 || v->ProjArgs[1]>90.0) {
printf("Lat2 (ProjArgs[1] out of range: %g\n", v->ProjArgs[1] );
invalid = 1;
}
if (v->ProjArgs[5]<=0.0) {
printf("ColInc (ProjArgs[5]) = %g (must be >=0.0)\n",
v->ProjArgs[5] );
invalid = 1;
}
break;
case 3: /* Stereographic */
if (v->ProjArgs[0]<-90.0 || v->ProjArgs[0]>90.0) {
printf("Central Latitude (ProjArgs[0]) out of range: ");
printf("%g (must be in +/-90)\n", v->ProjArgs[0] );
invalid = 1;
}
if (v->ProjArgs[1]<-180.0 || v->ProjArgs[1]>180.0) {
printf("Central Longitude (ProjArgs[1]) out of range: ");
printf("%g (must be in +/-180)\n", v->ProjArgs[1] );
invalid = 1;
}
if (v->ProjArgs[4]<0) {
printf("Column spacing (ProjArgs[4]) = %g (must be positive)\n",
v->ProjArgs[4]);
invalid = 1;
}
break;
case 4: /* Rotated */
/* WLH 4-21-95 */
if (v->ProjArgs[2]<=0.0) {
printf("Row Increment (ProjArgs[2]) = %g (must be >=0.0)\n",
v->ProjArgs[2] );
invalid = 1;
}
if (v->ProjArgs[3]<=0.0) {
printf("Column Increment = (ProjArgs[3]) %g (must be >=0.0)\n",
v->ProjArgs[3] );
invalid = 1;
}
if (v->ProjArgs[4]<-90.0 || v->ProjArgs[4]>90.0) {
printf("Central Latitude (ProjArgs[4]) out of range: ");
printf("%g (must be in +/-90)\n", v->ProjArgs[4] );
invalid = 1;
}
if (v->ProjArgs[5]<-180.0 || v->ProjArgs[5]>180.0) {
printf("Central Longitude (ProjArgs[5]) out of range: ");
printf("%g (must be in +/-180)\n", v->ProjArgs[5] );
invalid = 1;
}
if (v->ProjArgs[6]<-180.0 || v->ProjArgs[6]>180.0) {
printf("Central Longitude (ProjArgs[6]) out of range: ");
printf("%g (must be in +/-180)\n", v->ProjArgs[6] );
invalid = 1;
}
break;
case 5: /*Mercator*/
if (v->ProjArgs[2] == 0.0){
printf("Row Increment(Km) can not be 0.0\n");
invalid = 1;
}
if (v->ProjArgs[3] == 0.0){
printf("Column Increment(Km) can not be 0.0\n");
invalid = 1;
}
break;
default:
printf("Projection = %d, must be in 0..4\n", v->Projection );
invalid = 1;
}
return !invalid;
}
/*
* Get the McIDAS file number and grid number associated with the grid
* identified by time and var.
* Input: v - v5d grid struct
* time, var - timestep and variable of grid
* Output: mcfile, mcgrid - McIDAS grid file number and grid number
*/
int v5dGetMcIDASgrid( v5dstruct *v, int time, int var,
int *mcfile, int *mcgrid )
{
if (time<0 || time>=v->NumTimes) {
printf("Bad time argument to v5dGetMcIDASgrid: %d\n", time );
return 0;
}
if (var<0 || var>=v->NumVars) {
printf("Bad var argument to v5dGetMcIDASgrid: %d\n", var );
return 0;
}
*mcfile = (int) v->McFile[time][var];
*mcgrid = (int) v->McGrid[time][var];
return 1;
}
/*
* Set the McIDAS file number and grid number associated with the grid
* identified by time and var.
* Input: v - v5d grid struct
* time, var - timestep and variable of grid
* mcfile, mcgrid - McIDAS grid file number and grid number
* Return: 1 = ok, 0 = error (bad time or var)
*/
int v5dSetMcIDASgrid( v5dstruct *v, int time, int var,
int mcfile, int mcgrid )
{
if (time<0 || time>=v->NumTimes) {
printf("Bad time argument to v5dSetMcIDASgrid: %d\n", time );
return 0;
}
if (var<0 || var>=v->NumVars) {
printf("Bad var argument to v5dSetMcIDASgrid: %d\n", var );
return 0;
}
v->McFile[time][var] = (short) mcfile;
v->McGrid[time][var] = (short) mcgrid;
return 1;
}
/**********************************************************************/
/***** Input Functions *****/
/**********************************************************************/
/*
* Read the header from a COMP* file and return results in the v5dstruct.
* Input: f - the file descriptor
* v - pointer to a v5dstruct.
* Return: 1 = ok, 0 = error.
*/
static int read_comp_header( int f, v5dstruct *v )
{
unsigned int id;
/* reset file position to start of file */
lseek( f, 0, SEEK_SET );
/* read file ID */
read_int4( f, (int *) &id );
if (id==0x80808080 || id==0x80808081) {
/* Older COMP5D format */
int gridtimes, gridparms;
int i, j, it, iv, nl;
int gridsize;
float hgttop, hgtinc;
/*char *compgrid;*/
if (id==0x80808080) {
/* 20 vars, 300 times */
gridtimes = 300;
gridparms = 20;
}
else {
/* 30 vars, 400 times */
gridtimes = 400;
gridparms = 30;
}
v->FirstGridPos = 12*4 + 8*gridtimes + 4*gridparms;
read_int4( f, &v->NumTimes );
read_int4( f, &v->NumVars );
read_int4( f, &v->Nr );
read_int4( f, &v->Nc );
read_int4( f, &nl );
for (i=0;i<v->NumVars;i++) {
v->Nl[i] = nl;
v->LowLev[i] = 0;
}
read_float4( f, &v->ProjArgs[0] );
read_float4( f, &v->ProjArgs[1] );
read_float4( f, &hgttop );
read_float4( f, &v->ProjArgs[2] );
read_float4( f, &v->ProjArgs[3] );
read_float4( f, &hgtinc );
/*
for (i=0;i<nl;i++) {
v->Height[nl-i-1] = hgttop - i * hgtinc;
}
*/
v->VerticalSystem = 1;
v->VertArgs[0] = hgttop - hgtinc * (nl-1);
v->VertArgs[1] = hgtinc;
/* read dates and times */
for (i=0;i<gridtimes;i++) {
read_int4( f, &j );
v->DateStamp[i] = v5dDaysToYYDDD( j );
}
for (i=0;i<gridtimes;i++) {
read_int4( f, &j );
v->TimeStamp[i] = v5dSecondsToHHMMSS( j );
}
/* read variable names */
for (i=0;i<gridparms;i++) {
char name[4];
read_bytes( f, name, 4 );
/* remove trailing spaces, if any */
for (j=3;j>0;j--) {
if (name[j]==' ' || name[j]==0)
name[j] = 0;
else
break;
}
strncpy( v->VarName[i], name, 4 );
v->VarName[i][4] = 0;
}
gridsize = ( (v->Nr * v->Nc * nl + 3) / 4) * 4;
for (i=0;i<v->NumVars;i++) {
v->GridSize[i] = 8 + gridsize;
}
v->SumGridSizes = (8+gridsize) * v->NumVars;
/* read the grids and their ga,gb values to find min and max values */
for (i=0;i<v->NumVars;i++) {
v->MinVal[i] = 999999.9;
v->MaxVal[i] = -999999.9;
}
/*compgrid = (char *) malloc( gridsize );*/
for (it=0; it<v->NumTimes; it++) {
for (iv=0; iv<v->NumVars; iv++) {
float ga, gb;
float min, max;
read_float4( f, &ga );
read_float4( f, &gb );
/* skip ahead by 'gridsize' bytes */
if (lseek( f, gridsize, SEEK_CUR )==-1) {
printf("Error: Unexpected end of file, ");
printf("file may be corrupted.\n");
return 0;
}
min = -(125.0+gb)/ga;
max = (125.0-gb)/ga;
if (min<v->MinVal[iv]) v->MinVal[iv] = min;
if (max>v->MaxVal[iv]) v->MaxVal[iv] = max;
}
}
/*free( compgrid );*/
/* done */
}
else if (id==0x80808082 || id==0x80808083) {
/* Newer COMP5D format */
int gridtimes, gridsize;
int it, iv, nl, i, j;
float delta;
read_int4( f, &gridtimes );
read_int4( f, &v->NumVars );
read_int4( f, &v->NumTimes );
read_int4( f, &v->Nr );
read_int4( f, &v->Nc );
read_int4( f, &nl );
for (i=0;i<v->NumVars;i++) {
v->Nl[i] = nl;
}
read_float4( f, &v->ProjArgs[2] );
read_float4( f, &v->ProjArgs[3] );
/* Read height and determine if equal spacing */
v->VerticalSystem = 1;
for (i=0;i<nl;i++) {
read_float4( f, &v->VertArgs[i] );
if (i==1) {
delta = v->VertArgs[1] - v->VertArgs[0];
}
else if (i>1) {
if (delta != (v->VertArgs[i] - v->VertArgs[i-1])) {
v->VerticalSystem = 2;
}
}
}
if (v->VerticalSystem==1) {
v->VertArgs[1] = delta;
}
/* read variable names */
for (iv=0; iv<v->NumVars; iv++) {
char name[8];
read_bytes( f, name, 8 );
/* remove trailing spaces, if any */
for (j=7;j>0;j--) {
if (name[j]==' ' || name[j]==0)
name[j] = 0;
else
break;
}
strncpy( v->VarName[iv], name, 8 );
v->VarName[iv][8] = 0;
}
for (iv=0;iv<v->NumVars;iv++) {
read_float4( f, &v->MinVal[iv] );
}
for (iv=0;iv<v->NumVars;iv++) {
read_float4( f, &v->MaxVal[iv] );
}
for (it=0;it<gridtimes;it++) {
read_int4( f, &j );
v->TimeStamp[it] = v5dSecondsToHHMMSS( j );
}
for (it=0;it<gridtimes;it++) {
read_int4( f, &j );
v->DateStamp[it] = v5dDaysToYYDDD( j );
}
for (it=0;it<gridtimes;it++) {
float nlat;
read_float4( f, &nlat );
if (it==0) v->ProjArgs[0] = nlat;
}
for (it=0;it<gridtimes;it++) {
float wlon;
read_float4( f, &wlon );
if (it==0) v->ProjArgs[1] = wlon;
}
/* calculate grid storage sizes */
if (id==0x80808082) {
gridsize = nl*2*4 + ( (v->Nr * v->Nc * nl + 3) / 4) * 4;
}
else {
/* McIDAS grid and file numbers present */
gridsize = 8 + nl*2*4 + ( (v->Nr * v->Nc * nl + 3) / 4) * 4;
}
for (i=0;i<v->NumVars;i++) {
v->GridSize[i] = gridsize;
}
v->SumGridSizes = gridsize * v->NumVars;
/* read McIDAS numbers??? */
/* size (in bytes) of all header info */
v->FirstGridPos = 9*4 + v->Nl[0]*4 + v->NumVars*16 + gridtimes*16;
}
v->CompressMode = 1; /* one byte per grid point */
v->Projection = 1; /* Cylindrical equidistant */
v->FileVersion[0] = 0;
return 1;
}
/*
* Read a compressed grid from a COMP* file.
* Return: 1 = ok, 0 = error.
*/
static int read_comp_grid( v5dstruct *v, int time, int var,
float *ga, float *gb, void *compdata )
{
unsigned int pos;
V5Dubyte bias;
int i, n, nl;
int f;
V5Dubyte *compdata1 = (V5Dubyte *) compdata;
f = v->FileDesc;
/* move to position in file */
pos = grid_position( v, time, var );
lseek( f, pos, SEEK_SET );
if (v->FileFormat==0x80808083) {
/* read McIDAS grid and file numbers */
int mcfile, mcgrid;
read_int4( f, &mcfile );
read_int4( f, &mcgrid );
v->McFile[time][var] = (short) mcfile;
v->McGrid[time][var] = (short) mcgrid;
}
nl = v->Nl[var];
if (v->FileFormat==0x80808080 || v->FileFormat==0x80808081) {
/* single ga,gb pair for whole grid */
float a, b;
read_float4( f, &a );
read_float4( f, &b );
/* convert a, b to new v5d ga, gb values */
for (i=0;i<nl;i++) {
if (a==0.0) {
ga[i] = gb[i] = 0.0;
}
else {
gb[i] = (b+128.0) / -a;
ga[i] = 1.0 / a;
}
}
bias = 128;
}
else {
/* read ga, gb arrays */
read_float4_array( f, ga, v->Nl[var] );
read_float4_array( f, gb, v->Nl[var] );
/* convert ga, gb values to v5d system */
for (i=0;i<nl;i++) {
if (ga[i]==0.0) {
ga[i] = gb[i] = 0.0;
}
else {
/*gb[i] = (gb[i]+125.0) / -ga[i];*/
gb[i] = (gb[i]+128.0) / -ga[i];
ga[i] = 1.0 / ga[i];
}
}
bias = 128; /* 125 ??? */
}
/* read compressed grid data */
n = v->Nr * v->Nc * v->Nl[var];
if (read_bytes( f, compdata1, n )!=n)
return 0;
/* convert data values to v5d system */
n = v->Nr * v->Nc * v->Nl[var];
for (i=0;i<n;i++) {
compdata1[i] += bias;
}
return 1;
}
/*
* Read a v5d file header.
* Input: f - file opened for reading.
* v - pointer to v5dstruct to store header info into.
* Return: 1 = ok, 0 = error.
*/
static int read_v5d_header( v5dstruct *v )
{
#define SKIP(N) lseek( f, N, SEEK_CUR )
int end_of_header = 0;
unsigned int id;
int idlen, var, numargs;
int f;
f = v->FileDesc;
/* first try to read the header id */
read_int4( f, (int*) &id );
read_int4( f, &idlen );
if (id==TAG_ID && idlen==0) {
/* this is a v5d file */
v->FileFormat = 0;
}
else if (id>=0x80808080 && id<=0x80808083) {
/* this is an old COMP* file */
v->FileFormat = id;
return read_comp_header( f, v );
}
else {
/* unknown file type */
printf("Error: not a v5d file\n");
return 0;
}
v->CompressMode = 1; /* default */
while (!end_of_header) {
int tag, length;
int i, var, time, nl, lev;
if (read_int4(f,&tag)<1 || read_int4(f,&length)<1) {
printf("Error while reading header, premature EOF\n");
return 0;
}
switch (tag) {
case TAG_VERSION:
assert( length==10 );
read_bytes( f, v->FileVersion, 10 );
/* Check if reading a file made by a future version of Vis5D */
if (strcmp(v->FileVersion, FILE_VERSION)>0) {
/* WLH 6 Oct 98 */
printf("Warning: Trying to read a version %s file,", v->FileVersion);
printf(" you should upgrade Vis5D.\n");
}
break;
case TAG_NUMTIMES:
assert( length==4 );
read_int4( f, &v->NumTimes );
break;
case TAG_NUMVARS:
assert( length==4 );
read_int4( f, &v->NumVars );
break;
case TAG_VARNAME:
assert( length==14 ); /* 1 int + 10 char */
read_int4( f, &var );
read_bytes( f, v->VarName[var], 10 );
break;
case TAG_NR:
/* Number of rows for all variables */
assert( length==4 );
read_int4( f, &v->Nr );
break;
case TAG_NC:
/* Number of columns for all variables */
assert( length==4 );
read_int4( f, &v->Nc );
break;
case TAG_NL:
/* Number of levels for all variables */
assert( length==4 );
read_int4( f, &nl );
for (i=0;i<v->NumVars;i++) {
v->Nl[i] = nl;
}
break;
case TAG_NL_VAR:
/* Number of levels for one variable */
assert( length==8 );
read_int4( f, &var );
read_int4( f, &v->Nl[var] );
break;
case TAG_LOWLEV_VAR:
/* Lowest level for one variable */
assert( length==8 );
read_int4( f, &var );
read_int4( f, &v->LowLev[var] );
break;
case TAG_TIME:
/* Time stamp for 1 timestep */
assert( length==8 );
read_int4( f, &time );
read_int4( f, &v->TimeStamp[time] );
break;
case TAG_DATE:
/* Date stamp for 1 timestep */
assert( length==8 );
read_int4( f, &time );
read_int4( f, &v->DateStamp[time] );
break;
case TAG_MINVAL:
/* Minimum value for a variable */
assert( length==8 );
read_int4( f, &var );
read_float4( f, &v->MinVal[var] );
break;
case TAG_MAXVAL:
/* Maximum value for a variable */
assert( length==8 );
read_int4( f, &var );
read_float4( f, &v->MaxVal[var] );
break;
case TAG_COMPRESS:
/* Compress mode */
assert( length==4 );
read_int4( f, &v->CompressMode );
break;
case TAG_UNITS:
/* physical units */
assert( length==24 );
read_int4( f, &var );
read_bytes( f, v->Units[var], 20 );
break;
/*
* Vertical coordinate system
*/
case TAG_VERTICAL_SYSTEM:
assert( length==4 );
read_int4( f, &v->VerticalSystem );
if (v->VerticalSystem<0 || v->VerticalSystem>3) {
printf("Error: bad vertical coordinate system: %d\n",
v->VerticalSystem );
}
break;
case TAG_VERT_ARGS:
read_int4( f, &numargs );
assert( numargs <= MAXVERTARGS );
read_float4_array( f, v->VertArgs, numargs );
assert( length==numargs*4+4 );
break;
case TAG_HEIGHT:
/* height of a grid level */
assert( length==8 );
read_int4( f, &lev );
read_float4( f, &v->VertArgs[lev] );
break;
case TAG_BOTTOMBOUND:
assert( length==4 );
read_float4( f, &v->VertArgs[0] );
break;
case TAG_LEVINC:
assert( length==4 );
read_float4( f, &v->VertArgs[1] );
break;
/*
* Map projection information
*/
case TAG_PROJECTION:
assert( length==4 );
read_int4( f, &v->Projection );
if (v->Projection<0 || v->Projection>5) { /* WLH 4-21-95 */
printf("Error while reading header, bad projection (%d)\n",
v->Projection );
return 0;
}
break;
case TAG_PROJ_ARGS:
read_int4( f, &numargs );
assert( numargs <= MAXPROJARGS );
read_float4_array( f, v->ProjArgs, numargs );
assert( length==4*numargs+4 );
break;
case TAG_NORTHBOUND:
assert( length==4 );
if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
read_float4( f, &v->ProjArgs[0] );
}
else {
SKIP( 4 );
}
break;
case TAG_WESTBOUND:
assert( length==4 );
if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
read_float4( f, &v->ProjArgs[1] );
}
else {
SKIP( 4 );
}
break;
case TAG_ROWINC:
assert( length==4 );
if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
read_float4( f, &v->ProjArgs[2] );
}
else {
SKIP( 4 );
}
break;
case TAG_COLINC:
assert( length==4 );
if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
read_float4( f, &v->ProjArgs[3] );
}
else if (v->Projection==2) {
read_float4( f, &v->ProjArgs[5] );
}
else if (v->Projection==3) {
read_float4( f, &v->ProjArgs[4] );
}
else {
SKIP( 4 );
}
break;
case TAG_ROWINCKM:
assert( length==4 );
if (v->Projection==5){
read_float4( f, &v->ProjArgs[2] );
}
else {
SKIP( 4 );
}
break;
case TAG_COLINCKM:
assert( length==4 );
if (v->Projection==5){
read_float4( f, &v->ProjArgs[3] );
}
else {
SKIP( 4 );
}
break;
case TAG_LAT1:
assert( length==4 );
if (v->Projection==2) {
read_float4( f, &v->ProjArgs[0] );
}
else {
SKIP( 4 );
}
break;
case TAG_LAT2:
assert( length==4 );
if (v->Projection==2) {
read_float4( f, &v->ProjArgs[1] );
}
else {
SKIP( 4 );
}
break;
case TAG_POLE_ROW:
assert( length==4 );
if (v->Projection==2) {
read_float4( f, &v->ProjArgs[2] );
}
else {
SKIP( 4 );
}
break;
case TAG_POLE_COL:
assert( length==4 );
if (v->Projection==2) {
read_float4( f, &v->ProjArgs[3] );
}
else {
SKIP( 4 );
}
break;
case TAG_CENTLON:
assert( length==4 );
if (v->Projection==2) {
read_float4( f, &v->ProjArgs[4] );
}
else if (v->Projection==3 || v->Projection==5) {
read_float4( f, &v->ProjArgs[1] );
}
else if (v->Projection==4) { /* WLH 4-21-95 */
read_float4( f, &v->ProjArgs[5] );
}
else {
SKIP( 4 );
}
break;
case TAG_CENTLAT:
assert( length==4 );
if (v->Projection==3 || v->Projection==5) {
read_float4( f, &v->ProjArgs[0] );
}
else if (v->Projection==4) { /* WLH 4-21-95 */
read_float4( f, &v->ProjArgs[4] );
}
else {
SKIP( 4 );
}
break;
case TAG_CENTROW:
assert( length==4 );
if (v->Projection==3) {
read_float4( f, &v->ProjArgs[2] );
}
else {
SKIP( 4 );
}
break;
case TAG_CENTCOL:
assert( length==4 );
if (v->Projection==3) {
read_float4( f, &v->ProjArgs[3] );
}
else {
SKIP( 4 );
}
break;
case TAG_ROTATION:
assert( length==4 );
if (v->Projection==4) { /* WLH 4-21-95 */
read_float4( f, &v->ProjArgs[6] );
}
else {
SKIP( 4 );
}
break;
case TAG_END:
/* end of header */
end_of_header = 1;
lseek( f, length, SEEK_CUR );
break;
default:
/* unknown tag, skip to next tag */
printf("Unknown tag: %d length=%d\n", tag, length );
lseek( f, length, SEEK_CUR );
break;
}
}
v5dVerifyStruct( v );
/* Now we're ready to read the grid data */
/* Save current file pointer */
v->FirstGridPos = ltell(f);
/* compute grid sizes */
v->SumGridSizes = 0;
for (var=0;var<v->NumVars;var++) {
v->GridSize[var] = 8 * v->Nl[var] + v5dSizeofGrid( v, 0, var );
v->SumGridSizes += v->GridSize[var];
}
return 1;
#undef SKIP
}
/*
* Open a v5d file for reading.
* Input: filename - name of v5d file to open
* v - pointer to a v5dstruct in which to put header info or NULL
* if a struct should be dynamically allocated.
* Return: NULL if error, else v or a pointer to a new v5dstruct if v was NULL
*/
v5dstruct *v5dOpenFile( const char *filename, v5dstruct *v )
{
int fd;
fd = open( filename, O_RDONLY );
if (fd==-1) {
/* error */
return 0;
}
if (v) {
v5dInitStruct( v );
}
else {
v = v5dNewStruct();
if (!v) {
return NULL;
}
}
v->FileDesc = fd;
v->Mode = 'r';
if (read_v5d_header( v )) {
return v;
}
else {
return NULL;
}
}
/*
* Read a compressed grid from a v5d file.
* Input: v - pointer to v5dstruct describing the file
* time, var - which timestep and variable
* ga, gb - arrays to store grid (de)compression values
* compdata - address of where to store compressed grid data.
* Return: 1 = ok, 0 = error.
*/
int v5dReadCompressedGrid( v5dstruct *v, int time, int var,
float *ga, float *gb, void *compdata )
{
int pos, n, k;
if (time<0 || time>=v->NumTimes) {
printf("Error in v5dReadCompressedGrid: bad timestep argument (%d)\n",
time);
return 0;
}
if (var<0 || var>=v->NumVars) {
printf("Error in v5dReadCompressedGrid: bad var argument (%d)\n",
var);
return 0;
}
if (v->FileFormat) {
/* old COMP* file */
return read_comp_grid( v, time, var, ga, gb, compdata );
}
/* move to position in file */
pos = grid_position( v, time, var );
lseek( v->FileDesc, pos, SEEK_SET );
/* read ga, gb arrays */
read_float4_array( v->FileDesc, ga, v->Nl[var] );
read_float4_array( v->FileDesc, gb, v->Nl[var] );
/* read compressed grid data */
n = v->Nr * v->Nc * v->Nl[var];
if (v->CompressMode==1) {
k = read_block( v->FileDesc, compdata, n, 1 )==n;
}
else if (v->CompressMode==2) {
k = read_block( v->FileDesc, compdata, n, 2 )==n;
}
else if (v->CompressMode==4) {
k = read_block( v->FileDesc, compdata, n, 4 )==n;
}
if (!k) {
/* error */
printf("Error in v5dReadCompressedGrid: read failed, bad file?\n");
}
return k;
/*
n = v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
if (read( v->FileDesc, compdata, n )==n)
return 1;
else
return 0;
*/
}
/*
* Read a grid from a v5d file, decompress it and return it.
* Input: v - pointer to v5dstruct describing file header
* time, var - which timestep and variable.
* data - address of buffer to put grid data
* Output: data - the grid data
* Return: 1 = ok, 0 = error.
*/
int v5dReadGrid( v5dstruct *v, int time, int var, float data[] )
{
float ga[MAXLEVELS], gb[MAXLEVELS];
void *compdata;
int bytes;
if (time<0 || time>=v->NumTimes) {
printf("Error in v5dReadGrid: bad timestep argument (%d)\n", time);
return 0;
}
if (var<0 || var>=v->NumVars) {
printf("Error in v5dReadGrid: bad variable argument (%d)\n", var);
return 0;
}
/* allocate compdata buffer */
if (v->CompressMode==1) {
bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned char);
}
else if (v->CompressMode==2) {
bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned short);
}
else if (v->CompressMode==4) {
bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(float);
}
compdata = (void *) malloc( bytes );
if (!compdata) {
printf("Error in v5dReadGrid: out of memory (needed %d bytes)\n", bytes);
return 0;
}
/* read the compressed data */
if (!v5dReadCompressedGrid( v, time, var, ga, gb, compdata )) {
return 0;
}
/* decompress the data */
v5dDecompressGrid( v->Nr, v->Nc, v->Nl[var], v->CompressMode,
compdata, ga, gb, data );
/* free compdata */
free( compdata );
return 1;
}
/**********************************************************************/
/***** Output Functions *****/
/**********************************************************************/
static int write_tag( v5dstruct *v, int tag, int length, int newfile )
{
if (!newfile) {
/* have to check that there's room in header to write this tagged item */
if (v->CurPos+8+length > v->FirstGridPos) {
printf("Error: out of header space!\n");
/* Out of header space! */
return 0;
}
}
if (write_int4( v->FileDesc, tag )==0) return 0;
if (write_int4( v->FileDesc, length )==0) return 0;
v->CurPos += 8 + length;
return 1;
}
/*
* Write the information in the given v5dstruct as a v5d file header.
* Note that the current file position is restored when this function
* returns normally.
* Input: f - file already open for writing
* v - pointer to v5dstruct
* Return: 1 = ok, 0 = error.
*/
static int write_v5d_header( v5dstruct *v )
{
int var, time, filler, maxnl;
int f;
int newfile;
if (v->FileFormat!=0) {
printf("Error: v5d library can't write comp5d format files.\n");
return 0;
}
f = v->FileDesc;
if (!v5dVerifyStruct( v ))
return 0;
/* Determine if we're writing to a new file */
if (v->FirstGridPos==0) {
newfile = 1;
}
else {
newfile = 0;
}
/* compute grid sizes */
v->SumGridSizes = 0;
for (var=0;var<v->NumVars;var++) {
v->GridSize[var] = 8 * v->Nl[var] + v5dSizeofGrid( v, 0, var );
v->SumGridSizes += v->GridSize[var];
}
/* set file pointer to start of file */
lseek( f, 0, SEEK_SET );
v->CurPos = 0;
/*
* Write the tagged header info
*/
#define WRITE_TAG( V, T, L ) if (!write_tag(V,T,L,newfile)) return 0;
/* ID */
WRITE_TAG( v, TAG_ID, 0 );
/* File Version */
WRITE_TAG( v, TAG_VERSION, 10 );
write_bytes( f, FILE_VERSION, 10 );
/* Number of timesteps */
WRITE_TAG( v, TAG_NUMTIMES, 4 );
write_int4( f, v->NumTimes );
/* Number of variables */
WRITE_TAG( v, TAG_NUMVARS, 4 );
write_int4( f, v->NumVars );
/* Names of variables */
for (var=0;var<v->NumVars;var++) {
WRITE_TAG( v, TAG_VARNAME, 14 );
write_int4( f, var );
write_bytes( f, v->VarName[var], 10 );
}
/* Physical Units */
for (var=0;var<v->NumVars;var++) {
WRITE_TAG( v, TAG_UNITS, 24 );
write_int4( f, var );
write_bytes( f, v->Units[var], 20 );
}
/* Date and time of each timestep */
for (time=0;time<v->NumTimes;time++) {
WRITE_TAG( v, TAG_TIME, 8 );
write_int4( f, time );
write_int4( f, v->TimeStamp[time] );
WRITE_TAG( v, TAG_DATE, 8 );
write_int4( f, time );
write_int4( f, v->DateStamp[time] );
}
/* Number of rows */
WRITE_TAG( v, TAG_NR, 4 );
write_int4( f, v->Nr );
/* Number of columns */
WRITE_TAG( v, TAG_NC, 4 );
write_int4( f, v->Nc );
/* Number of levels, compute maxnl */
maxnl = 0;
for (var=0;var<v->NumVars;var++) {
WRITE_TAG( v, TAG_NL_VAR, 8 );
write_int4( f, var );
write_int4( f, v->Nl[var] );
WRITE_TAG( v, TAG_LOWLEV_VAR, 8 );
write_int4( f, var );
write_int4( f, v->LowLev[var] );
if (v->Nl[var]+v->LowLev[var]>maxnl) {
maxnl = v->Nl[var]+v->LowLev[var];
}
}
/* Min/Max values */
for (var=0;var<v->NumVars;var++) {
WRITE_TAG( v, TAG_MINVAL, 8 );
write_int4( f, var );
write_float4( f, v->MinVal[var] );
WRITE_TAG( v, TAG_MAXVAL, 8 );
write_int4( f, var );
write_float4( f, v->MaxVal[var] );
}
/* Compress mode */
WRITE_TAG( v, TAG_COMPRESS, 4 );
write_int4( f, v->CompressMode );
/* Vertical Coordinate System */
WRITE_TAG( v, TAG_VERTICAL_SYSTEM, 4 );
write_int4( f, v->VerticalSystem );
WRITE_TAG( v, TAG_VERT_ARGS, 4+4*MAXVERTARGS );
write_int4( f, MAXVERTARGS );
write_float4_array( f, v->VertArgs, MAXVERTARGS );
/* Map Projection */
WRITE_TAG( v, TAG_PROJECTION, 4 );
write_int4( f, v->Projection );
WRITE_TAG( v, TAG_PROJ_ARGS, 4+4*MAXPROJARGS );
write_int4( f, MAXPROJARGS );
write_float4_array( f, v->ProjArgs, MAXPROJARGS );
/* write END tag */
if (newfile) {
/* We're writing to a brand new file. Reserve 10000 bytes */
/* for future header growth. */
WRITE_TAG( v, TAG_END, 10000 );
lseek( f, 10000, SEEK_CUR );
/* Let file pointer indicate where first grid is stored */
v->FirstGridPos = ltell( f );
}
else {
/* we're rewriting a header */
filler = v->FirstGridPos - ltell(f);
WRITE_TAG( v, TAG_END, filler-8 );
}
#undef WRITE_TAG
return 1;
}
/*
* Open a v5d file for writing. If the named file already exists,
* it will be deleted.
* Input: filename - name of v5d file to create.
* v - pointer to v5dstruct with the header info to write.
* Return: 1 = ok, 0 = error.
*/
int v5dCreateFile( const char *filename, v5dstruct *v )
{
mode_t mask;
int fd;
mask = 0666;
fd = open( filename, O_WRONLY | O_CREAT | O_TRUNC, mask );
if (fd==-1) {
printf("Error in v5dCreateFile: open failed\n");
v->FileDesc = -1;
v->Mode = 0;
return 0;
}
else {
/* ok */
v->FileDesc = fd;
v->Mode = 'w';
/* write header and return status */
return write_v5d_header(v);
}
}
/*
* Open a v5d file for updating/appending and read the header info.
* Input: filename - name of v5d file to open for updating.
* v - pointer to v5dstruct in which the file header info will be
* put. If v is NULL a v5dstruct will be allocated and returned.
* Return: NULL if error, else v or a pointer to a new v5dstruct if v as NULL
*/
v5dstruct *v5dUpdateFile( const char *filename, v5dstruct *v )
{
int fd;
fd = open( filename, O_RDWR );
if (fd==-1) {
return NULL;
}
if (!v) {
v = v5dNewStruct();
if (!v) {
return NULL;
}
}
v->FileDesc = fd;
v->Mode = 'w';
if (read_v5d_header( v )) {
return v;
}
else {
return NULL;
}
}
/*
* Write a compressed grid to a v5d file.
* Input: v - pointer to v5dstruct describing the file
* time, var - which timestep and variable
* ga, gb - the GA and GB (de)compression value arrays
* compdata - address of array of compressed data values
* Return: 1 = ok, 0 = error.
*/
int v5dWriteCompressedGrid( const v5dstruct *v, int time, int var,
const float *ga, const float *gb,
const void *compdata )
{
int pos, n, k;
/* simple error checks */
if (v->Mode!='w') {
printf("Error in v5dWriteCompressedGrid: file opened for reading,");
printf(" not writing.\n");
return 0;
}
if (time<0 || time>=v->NumTimes) {
printf("Error in v5dWriteCompressedGrid: bad timestep argument (%d)\n",
time);
return 0;
}
if (var<0 || var>=v->NumVars) {
printf("Error in v5dWriteCompressedGrid: bad variable argument (%d)\n",
var);
return 0;
}
/* move to position in file */
pos = grid_position( v, time, var );
if (lseek( v->FileDesc, pos, SEEK_SET )<0) {
/* lseek failed, return error */
printf("Error in v5dWrite[Compressed]Grid: seek failed, disk full?\n");
return 0;
}
/* write ga, gb arrays */
k = 0;
if (write_float4_array( v->FileDesc, ga, v->Nl[var] ) == v->Nl[var] &&
write_float4_array( v->FileDesc, gb, v->Nl[var] ) == v->Nl[var]) {
/* write compressed grid data (k=1=OK, k=0=Error) */
n = v->Nr * v->Nc * v->Nl[var];
if (v->CompressMode==1) {
k = write_block( v->FileDesc, compdata, n, 1 )==n;
}
else if (v->CompressMode==2) {
k = write_block( v->FileDesc, compdata, n, 2 )==n;
}
else if (v->CompressMode==4) {
k = write_block( v->FileDesc, compdata, n, 4 )==n;
}
}
if (k==0) {
/* Error while writing */
printf("Error in v5dWrite[Compressed]Grid: write failed, disk full?\n");
}
return k;
/*
n = v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
if (write_bytes( v->FileDesc, compdata, n )!=n) {
printf("Error in v5dWrite[Compressed]Grid: write failed, disk full?\n");
return 0;
}
else {
return 1;
}
*/
}
/*
* Compress a grid and write it to a v5d file.
* Input: v - pointer to v5dstruct describing the file
* time, var - which timestep and variable (starting at 0)
* data - address of uncompressed grid data
* Return: 1 = ok, 0 = error.
*/
int v5dWriteGrid( v5dstruct *v, int time, int var, const float data[] )
{
float ga[MAXLEVELS], gb[MAXLEVELS];
void *compdata;
int n, bytes;
float min, max;
if (v->Mode!='w') {
printf("Error in v5dWriteGrid: file opened for reading,");
printf(" not writing.\n");
return 0;
}
if (time<0 || time>=v->NumTimes) {
printf("Error in v5dWriteGrid: bad timestep argument (%d)\n", time);
return 0;
}
if (var<0 || var>=v->NumVars) {
printf("Error in v5dWriteGrid: bad variable argument (%d)\n", var);
return 0;
}
/* allocate compdata buffer */
if (v->CompressMode==1) {
bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned char);
}
else if (v->CompressMode==2) {
bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned short);
}
else if (v->CompressMode==4) {
bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(float);
}
compdata = (void *) malloc( bytes );
if (!compdata) {
printf("Error in v5dWriteGrid: out of memory (needed %d bytes)\n",
bytes );
return 0;
}
/* compress the grid data */
v5dCompressGrid( v->Nr, v->Nc, v->Nl[var], v->CompressMode, data,
compdata, ga, gb, &min, &max );
/* update min and max value */
if (min<v->MinVal[var])
v->MinVal[var] = min;
if (max>v->MaxVal[var])
v->MaxVal[var] = max;
/* write the compressed grid */
n = v5dWriteCompressedGrid( v, time, var, ga, gb, compdata );
/* free compdata */
free( compdata );
return n;
}
/*
* Close a v5d file which was opened with open_v5d_file() or
* create_v5d_file().
* Input: f - file descriptor
* Return: 1 = ok, 0 = error
*/
int v5dCloseFile( v5dstruct *v )
{
int status = 1;
if (v->Mode=='w') {
/* rewrite header because writing grids updates the minval and */
/* maxval fields */
lseek( v->FileDesc, 0, SEEK_SET );
status = write_v5d_header( v );
lseek( v->FileDesc, 0, SEEK_END );
close( v->FileDesc );
}
else if (v->Mode=='r') {
/* just close the file */
close(v->FileDesc);
}
else {
printf("Error in v5dCloseFile: bad v5dstruct argument\n");
return 0;
}
v->FileDesc = -1;
v->Mode = 0;
return status;
}
/**********************************************************************/
/***** Simple v5d file writing functions. *****/
/**********************************************************************/
static v5dstruct *Simple = NULL;
/*
* Create a new v5d file specifying both a map projection and vertical
* coordinate system. See README file for argument details.
* Return: 1 = ok, 0 = error.
*/
int v5dCreate( const char *name, int numtimes, int numvars,
int nr, int nc, const int nl[],
const char varname[MAXVARS][10],
const int timestamp[], const int datestamp[],
int compressmode,
int projection,
const float proj_args[],
int vertical,
const float vert_args[] )
{
int var, time, maxnl, i;
/* initialize the v5dstruct */
Simple = v5dNewStruct();
printf("Intialized v5dstruct\n");
Simple->NumTimes = numtimes;
Simple->NumVars = numvars;
Simple->Nr = nr;
Simple->Nc = nc;
maxnl = nl[0];
for (var=0;var<numvars;var++) {
if (nl[var]>maxnl) {
maxnl = nl[var];
}
Simple->Nl[var] = nl[var];
Simple->LowLev[var] = 0;
strncpy( Simple->VarName[var], varname[var], 10 );
Simple->VarName[var][9] = 0;
}
/* time and date for each timestep */
for (time=0;time<numtimes;time++) {
Simple->TimeStamp[time] = timestamp[time];
Simple->DateStamp[time] = datestamp[time];
}
Simple->CompressMode = compressmode;
/* Map projection and vertical coordinate system */
Simple->Projection = projection;
memcpy( Simple->ProjArgs, proj_args, MAXPROJARGS*sizeof(float) );
Simple->VerticalSystem = vertical;
if (vertical == 3) {
/* convert pressures to heights */
for (i=0; i<MAXVERTARGS; i++) {
if (vert_args[i] > 0.000001) {
Simple->VertArgs[i] = pressure_to_height(vert_args[i]);
}
else Simple->VertArgs[i] = 0.0;
}
}
else {
memcpy( Simple->VertArgs, vert_args, MAXVERTARGS*sizeof(float) );
}
printf("before create file \n");
/* create the file */
if (v5dCreateFile( name, Simple )==0) {
printf("Error in v5dCreateSimpleFile: unable to create %s\n", name );
return 0;
}
else {
printf("before return \n");
return 1;
}
}
/*
* Create a new v5d file using minimal information.
* Return: 1 = ok, 0 = error. See README file for argument details.
*/
int v5dCreateSimple( const char *name, int numtimes, int numvars,
int nr, int nc, int nl,
const char varname[MAXVARS][10],
const int timestamp[], const int datestamp[],
float northlat, float latinc,
float westlon, float loninc,
float bottomhgt, float hgtinc )
{
int nlvar[MAXVARS];
int compressmode, projection, vertical;
float proj_args[100], vert_args[MAXLEVELS];
int i;
for (i=0;i<numvars;i++) {
nlvar[i] = nl;
}
compressmode = 1;
projection = 1;
proj_args[0] = northlat;
proj_args[1] = westlon;
proj_args[2] = latinc;
proj_args[3] = loninc;
vertical = 1;
vert_args[0] = bottomhgt;
vert_args[1] = hgtinc;
return v5dCreate( name, numtimes, numvars, nr, nc, nlvar,
varname, timestamp, datestamp, compressmode,
projection, proj_args, vertical, vert_args );
}
/*
* Set lowest levels for each variable (other than default of 0).
* Input: lowlev - array [NumVars] of ints
* Return: 1 = ok, 0 = error
*/
int v5dSetLowLev( int lowlev[] )
{
int var;
if (Simple) {
for (var=0;var<Simple->NumVars;var++) {
Simple->LowLev[var] = lowlev[var];
}
return 1;
}
else {
printf("Error: must call v5dCreate before v5dSetLowLev\n");
return 0;
}
}
/*
* Set the units for a variable.
* Input: var - a variable in [1,NumVars]
* units - a string
* Return: 1 = ok, 0 = error
*/
int v5dSetUnits( int var, const char *units )
{
if (Simple) {
if (var>=1 && var<=Simple->NumVars) {
strncpy( Simple->Units[var-1], units, 19 );
Simple->Units[var-1][19] = 0;
return 1;
}
else {
printf("Error: bad variable number in v5dSetUnits\n");
return 0;
}
}
else {
printf("Error: must call v5dCreate before v5dSetUnits\n");
return 0;
}
}
/*
* Write a grid to a v5d file.
* Input: time - timestep in [1,NumTimes]
* var - timestep in [1,NumVars]
* data - array [nr*nc*nl] of floats
* Return: 1 = ok, 0 = error
*/
int v5dWrite( int time, int var, const float data[] )
{
if (Simple) {
if (time<1 || time>Simple->NumTimes) {
printf("Error in v5dWrite: bad timestep number: %d\n", time );
return 0;
}
if (var<1 || var>Simple->NumVars) {
printf("Error in v5dWrite: bad variable number: %d\n", var );
}
return v5dWriteGrid( Simple, time-1, var-1, data );
}
else {
printf("Error: must call v5dCreate before v5dWrite\n");
return 0;
}
}
/*
* Close a v5d file after the last grid has been written to it.
* Return: 1 = ok, 0 = error
*/
int v5dClose( void )
{
if (Simple) {
int ok = v5dCloseFile( Simple );
v5dFreeStruct( Simple );
return ok;
}
else {
printf("Error: v5dClose: no file to close\n");
return 0;
}
}
/**********************************************************************/
/***** FORTRAN-callable simple output *****/
/**********************************************************************/
/*
* Create a v5d file. See README file for argument descriptions.
* Return: 1 = ok, 0 = error.
*/
#ifdef UNDERSCORE
int v5dcreate_
#else
# ifdef _CRAY
int V5DCREATE
# else
int v5dcreate
# endif
#endif
( const char *name, const int *numtimes, const int *numvars,
const int *nr, const int *nc, const int nl[],
const char varname[][10],
const int timestamp[], const int datestamp[],
const int *compressmode,
const int *projection,
const float proj_args[],
const int *vertical,
const float vert_args[] )
{
char filename[100];
char names[MAXVARS][10];
int i, maxnl, args;
/* copy name to filename and remove trailing spaces if any */
copy_string( filename, name, 100 );
/*
* Check for uninitialized arguments
*/
if (*numtimes<1) {
printf("Error: numtimes invalid\n");
return 0;
}
if (*numvars<1) {
printf("Error: numvars invalid\n");
return 0;
}
if (*nr<2) {
printf("Error: nr invalid\n");
return 0;
}
if (*nc<2) {
printf("Error: nc invalid\n");
return 0;
}
maxnl = 0;
for (i=0;i<*numvars;i++) {
if (nl[i]<1) {
printf("Error: nl(%d) invalid\n", i+1);
return 0;
}
if (nl[i]>maxnl) {
maxnl = nl[i];
}
}
for (i=0;i<*numvars;i++) {
if (copy_string2( names[i], varname[i], 10)==0) {
printf("Error: unitialized varname(%d)\n", i+1);
return 0;
}
}
for (i=0;i<*numtimes;i++) {
if (timestamp[i]<0) {
printf("Error: times(%d) invalid\n", i+1);
return 0;
}
if (datestamp[i]<0) {
printf("Error: dates(%d) invalid\n", i+1);
return 0;
}
}
if (*compressmode != 1 && *compressmode != 2 && *compressmode != 4) {
printf("Error: compressmode invalid\n");
return 0;
}
switch (*projection) {
case 0:
args = 4;
break;
case 1:
args = 0;
if (IS_MISSING(proj_args[0])) {
printf("Error: northlat (proj_args(1)) invalid\n");
return 0;
}
if (IS_MISSING(proj_args[1])) {
printf("Error: westlon (proj_args(2)) invalid\n");
return 0;
}
if (IS_MISSING(proj_args[2])) {
printf("Error: latinc (proj_args(3)) invalid\n");
return 0;
}
if (IS_MISSING(proj_args[3])) {
printf("Error: loninc (proj_args(4)) invalid\n");
return 0;
}
break;
case 2:
args = 6;
break;
case 3:
args = 5;
break;
case 4:
args = 7;
break;
/* MJK 12.12.99 */
case 5:
args = 4;
break;
default:
args = 0;
printf("Error: projection invalid\n");
return 0;
}
for (i=0;i<args;i++) {
if (IS_MISSING(proj_args[i])) {
printf("Error: proj_args(%d) invalid\n", i+1);
return 0;
}
}
switch (*vertical) {
case 0:
/* WLH 31 Oct 96 - just fall through
args = 4;
break;
*/
case 1:
args = 0;
if (IS_MISSING(vert_args[0])) {
printf("Error: bottomhgt (vert_args(1)) invalid\n");
return 0;
}
if (IS_MISSING(vert_args[1])) {
printf("Error: hgtinc (vert_args(2)) invalid\n");
return 0;
}
break;
case 2:
case 3:
args = maxnl;
break;
default:
args = 0;
printf("Error: vertical invalid\n");
return 0;
}
for (i=0;i<args;i++) {
if (IS_MISSING(vert_args[i])) {
printf("Error: vert_args(%d) invalid\n", i+1);
return 0;
}
}
return v5dCreate( filename, *numtimes, *numvars, *nr, *nc, nl,
(const char(*)[10]) names, timestamp, datestamp,
*compressmode,
*projection, proj_args, *vertical, vert_args );
}
/*
* Create a simple v5d file. See README file for argument descriptions.
* Return: 1 = ok, 0 = error.
*/
#ifdef UNDERSCORE
int v5dcreatesimple_
#else
# ifdef _CRAY
int V5DCREATESIMPLE
# else
int v5dcreatesimple
# endif
#endif
( const char *name, const int *numtimes, const int *numvars,
const int *nr, const int *nc, const int *nl,
const char varname[][10],
const int timestamp[], const int datestamp[],
const float *northlat, const float *latinc,
const float *westlon, const float *loninc,
const float *bottomhgt, const float *hgtinc )
{
int compressmode, projection, vertical;
float projarg[100], vertarg[MAXLEVELS];
int varnl[MAXVARS];
int i;
for (i=0;i<MAXVARS;i++) {
varnl[i] = *nl;
}
compressmode = 1;
projection = 1;
projarg[0] = *northlat;
projarg[1] = *westlon;
projarg[2] = *latinc;
projarg[3] = *loninc;
vertical = 1;
vertarg[0] = *bottomhgt;
vertarg[1] = *hgtinc;
#ifdef UNDERSCORE
return v5dcreate_
#else
# ifdef _CRAY
return V5DCREATE
# else
return v5dcreate
# endif
#endif
( name, numtimes, numvars, nr, nc, varnl,
varname, timestamp, datestamp, &compressmode,
&projection, projarg, &vertical, vertarg );
}
/*
* Set lowest levels for each variable (other than default of 0).
* Input: lowlev - array [NumVars] of ints
* Return: 1 = ok, 0 = error
*/
#ifdef UNDERSCORE
int v5dsetlowlev_
#else
# ifdef _CRAY
int V5DSETLOWLEV
# else
int v5dsetlowlev
# endif
#endif
( int *lowlev )
{
return v5dSetLowLev(lowlev);
}
/*
* Set the units for a variable.
* Input: var - variable number in [1,NumVars]
* units - a character string
* Return: 1 = ok, 0 = error
*/
#ifdef UNDERSCORE
int v5dsetunits_
#else
# ifdef _CRAY
int V5DSETUNITS
# else
int v5dsetunits
# endif
#endif
( int *var, char *name )
{
return v5dSetUnits( *var, name );
}
/*
* Write a grid of data to the file.
* Input: time - timestep in [1,NumTimes]
* var - timestep in [1,NumVars]
* data - array [nr*nc*nl] of floats
* Return: 1 = ok, 0 = error
*/
#ifdef UNDERSCORE
int v5dwrite_
#else
# ifdef _CRAY
int V5DWRITE
# else
int v5dwrite
# endif
#endif
( const int *time, const int *var, const float *data )
{
return v5dWrite( *time, *var, data );
}
/*
* Specify the McIDAS GR3D file number and grid number which correspond
* to the grid specified by time and var.
* Input: time, var - timestep and variable of grid (starting at 1)
* mcfile, mcgrid - McIDAS grid file number and grid number
* Return: 1 = ok, 0 = errror (bad time or var)
*/
#ifdef UNDERSCORE
int v5dmcfile_
#else
# ifdef _CRAY
int V5DMCFILE
# else
int v5dmcfile
# endif
#endif
( const int *time, const int *var,
const int *mcfile, const int *mcgrid )
{
if (*time<1 || *time>Simple->NumTimes) {
printf("Bad time argument to v5dSetMcIDASgrid: %d\n", *time );
return 0;
}
if (*var<1 || *var>Simple->NumVars) {
printf("Bad var argument to v5dSetMcIDASgrid: %d\n", *var );
return 0;
}
Simple->McFile[*time-1][*var-1] = (short) *mcfile;
Simple->McGrid[*time-1][*var-1] = (short) *mcgrid;
return 1;
}
/*
* Close a simple v5d file.
*/
#ifdef UNDERSCORE
int v5dclose_( void )
#else
# ifdef _CRAY
int V5DCLOSE( void )
# else
int v5dclose( void )
# endif
#endif
{
return v5dClose();
}