/*
* 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 2 of the License, or
* (at your option) any later version.
*
* As a special exception to the terms of the GNU General Public
* License, you are permitted to link Vis5D with (and distribute the
* resulting source and executables) the LUI library (copyright by
* Stellar Computer Inc. and licensed for distribution with Vis5D),
* the McIDAS library, and/or the NetCDF library, where those
* libraries are governed by the terms of their own licenses.
*
* 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include "../config.h"
/* read_epa.c */
/*
* Functions for working with EPA files
*/
#include <assert.h>
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "file_i.h"
#include "grid_i.h"
#include "misc_i.h"
#include "model_i.h"
#include "proj_i.h"
#include "projlist_i.h"
#include "v5d.h"
#ifdef EPA
#define LATLON_FILE "latlon.dat"
/* This is needed by mm4.c: */
float mm4Sigma[MAXLEVELS];
int model_type;
/* Sigma value for each level in the 3-D grid for
15 level RADM files (how general?) */
static float RadmSigma15[15] = {
0.995,
0.985,
0.970,
0.945,
0.910,
0.865,
0.810,
0.740,
0.650,
0.550,
0.450,
0.350,
0.250,
0.150,
0.050
};
/* Sigma value for each level in the 3-D grid for
6 level RADM files (how general?) */
static float RadmSigma6[6] = {
0.99,
0.955,
0.885,
0.720,
0.450,
0.150,
};
/*
* Convert a date from the EPA date/time string format, which is
* "mo/day/year hour:min:sec", to Vis5d format, which is two
* integers giving the number of days since 1 Jan 1900 and the
* number of seconds since midnight (usually GMT).
*/
static void convert_date( char *datestr, int *date, int *time )
{
int mo, day, year, hour, min, sec, i;
int days_per_month[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
sscanf(datestr, "%d/%d/%d %d:%d:%d",
&mo, &day, &year, &hour, &min, &sec);
/*
printf("mo: %d day: %d year: %d hour: %d min: %d sec: %d\n",
mo, day, year, hour, min, sec);
*/
if ((!(year % 4)) && (year % 100)) days_per_month[1] = 29;
for (i=0; i<mo-1; i++) {
day += days_per_month[i];
}
/*
printf("day: %d year: %d\n", day, year);
*/
*date = (year * 1461) / 4 + day;
*time = hour * 60 * 60 + min * 60 + sec;
return;
}
/*
* Get the lat/lon/height information for an EPA dataset.
* (This information exists in various formats depending on the
* model. For example, in some cases it is on a piece of paper
* somewhere in the top right hand desk drawer in the modeler's
* office.) Also get spatial bunds for the model domain.
*
* Input: nr, nc, nl - grid size
* numvars - number of variables
* varname - array of variable names
* projargs - buffer of 2*nr*nc floats
* vertargs - buffer of nl floats
* Output: proj_out - map projection struct pointer
* vcs_out - vertical coord system struct pointer
* Return: 1 = ok, 0 = error
*/
static int get_epa_projection( struct grid_db *db,
int radm_fd,
int nr, int nc, int nl,
int numvars,
char varname[MAXVARS][10],
struct projection **proj_out,
struct vcs **vcs_out )
{
float pargs[MAXROWS * MAXCOLUMNS * 2];
float vargs[MAXLEVELS * 2];
#define Latitude(ROW,COL) pargs[ ((ROW)*nc+(COL)) * 2 + 0 ]
#define Longitude(ROW,COL) pargs[ ((ROW)*nc+(COL)) * 2 + 1 ]
#define Sigma(LEVEL) vargs[ (LEVEL) ]
/*
float Latitude[MAXROWS][MAXCOLUMNS], Longitude[MAXROWS][MAXCOLUMNS];
float Sigma[MAXLEVELS];
*/
FILE *f;
int i, j, k, kk;
float lat, lon;
char message[80];
char species[80];
char units[80];
int guard;
if (model_type == RADM) {
/*
* Read the (row,column) to (lat,lon) conversion information from
* the LATLON_FILE:
*/
if (nr == 38 && nc == 35) {
f = fopen( LATLON_FILE, "r" );
if (!f) {
printf("Error: unable to open %s\n", LATLON_FILE );
return 0;
}
while (!feof(f)) {
if (fscanf( f, "%d %d %f %f\n", &j, &i, &lat, &lon )<4) {
break;
}
else {
if (i>nr || j>nc) {
printf("BAD I,J: % %d\n", i, j );
}
/* note that we flip the rows from North to South -
with Vis5d the first row is at the North, and
rows increase toward the South -
we also flip the data North to South after we
get it from model_get_data -
of course we could probably skip both of these flips
and let the resampling take care of it */
Latitude(nr-i,j-1) = lat;
Longitude(nr-i,j-1) = lon;
}
}
fclose(f);
}
else if (nr == 61 && nc == 69) {
float lowLat[MAXROWS][MAXCOLUMNS], lowLon[MAXROWS][MAXCOLUMNS];
f = fopen( LATLON_FILE,"r");
if (!f) {
printf("Error: unable to open %s\n", LATLON_FILE );
exit(1);
}
while (!feof(f)) {
if (fscanf( f, "%d %d %f %f\n", &j, &i, &lat, &lon )<4) {
break;
}
else {
if (i>38 || j>35) {
printf("BAD I,J: % %d\n", i, j );
}
lowLat[i-1][j-1] = lat;
lowLon[i-1][j-1] = lon;
}
}
for (i=0; i<nr; i++) {
for (j=0; j<nc; j++) {
float ra, ca;
int ir, jc;
ra = i / 4.0; /* covert 20 km res to 80 km res */
ca = j / 4.0;
ir = ra; /* get integer parts of 80 km grid coord */
jc = ca;
ra = ra - ir; /* get fractional part2 */
ca = ca - jc;
ir = 16 + ir; /* add offset for (1,1) -> (17,15) */
jc = 14 + jc;
/* bilinear interp of lat and lon */
lat = (1.0 - ra) * ((1.0 - ca) * lowLat[ir][jc] +
ca * lowLat[ir][jc+1]) +
ra * ((1.0 - ca) * lowLat[ir+1][jc] +
ca * lowLat[ir+1][jc+1]);
lon = (1.0 - ra) * ((1.0 - ca) * lowLon[ir][jc] +
ca * lowLon[ir][jc+1]) +
ra * ((1.0 - ca) * lowLon[ir+1][jc] +
ca * lowLon[ir+1][jc+1]);
/* note that we flip the rows from North to South -
with Vis5d the first row is at the North, and
rows increase toward the South -
we also flip the data North to South after we
get it from model_get_data -
of course we could probably skip both of these flips
and let the resampling take care of it */
Latitude(nr-1-i,j) = lat;
Longitude(nr-1-i,j) = lon;
}
}
fclose(f);
}
else { /* not "low res" or "hi res", so put it in the ocean */
for (i=0; i<nr; i++) {
for (j=0; j<nc; j++) {
Latitude(i,j) = 5.0 - 0.01 * i;
Longitude(i,j) = 140.0 - 0.01 * j;
}
}
}
if(nl == 15) {
for (k=0; k<nl; k++) Sigma(k) = RadmSigma15[k];
}
else if(nl == 6) {
for (k=0; k<nl; k++) Sigma(k) = RadmSigma6[k];
}
else {
printf("RADM number of sigma levels must be 6 or 15 (is %d)\n", nl);
exit(1);
}
}
else if (model_type == MM4) {
/* Read (row,col) -> (lat,lon) mappings from the data file itself. */
float *gridlat, *gridlon;
int iplat = -9, iplon = -9, ip, status;
for (ip=0;ip<numvars;ip++) {
if (strcmp(varname[ip],"LATC")==0) iplat = ip;
if (strcmp(varname[ip],"LONC")==0) iplon = ip;
}
/* printf("iplat: %d iplon: %d\n", iplat, iplon); */
if (iplat < 0 || iplon < 0) {
printf("bad iplat (%d) or iplon (%d)\n", iplat, iplon);
exit(1);
}
gridlat = (float *) malloc( nr * nc * nl * sizeof(float) );
gridlon = (float *) malloc( nr * nc * nl * sizeof(float) );
strcpy(species, varname[iplat]);
strcpy(message, "");
if (!(status = model_get_data(radm_fd, 1, iplat+1, gridlat,
species, units, message))) {
printf("Error while reading mm4 lat grid: %s\n", message);
exit(1);
}
/*
printf("mm4_get_data gridlat. status: %d, message: %s\n",
status, message);
*/
strcpy(species, varname[iplon]);
strcpy(message, "");
if (!(status = model_get_data(radm_fd, 1, iplon+1, gridlon,
species, units, message))) {
printf("Error while reading mm4 lat grid: %s\n", message);
exit(1);
}
/*
printf("mm4_get_data gridlon. status: %d, message: %s\n",
status, message);
*/
/*
for (i=0; i<nr*nc; i += 200) {
printf("gridlat[%d] = %f, gridlon[%d] = %f\n",
i, gridlat[i], i, gridlon[i]);
}
*/
for (i=0; i<nr; i++) {
for (j=0; j<nc; j++) {
/* note that we flip the rows from North to South -
with Vis5d the first row is at the North, and
rows increase toward the South -
we also flip the data North to South after we
get it from model_get_data -
of course we could probably skip both of these flips
and let the resampling take care of it */
Latitude((nr-1)-i, j) = gridlat[(i*nc) + j];
Longitude((nr-1)-i, j) = -gridlon[(i*nc) + j];
}
}
free(gridlat);
free(gridlon);
for (k=0; k<nl; k++) {
Sigma(k) = (mm4Sigma[k] + mm4Sigma[k+1]) / 2.0;
}
}
else {
printf("lat/lon not yet supported for model_type %d\n", model_type);
exit(1);
}
#undef Latitude
#undef Longitude
#undef Sigma
/* return projection and vcs struct pointers */
*proj_out = new_projection( db, PROJ_EPA, nr, nc, pargs );
*vcs_out = new_vcs( db, VERT_EPA, nl, 0, vargs );
return 1;
}
#endif /*EPA*/
/*
* Scan the named file, createing a grid_info struct for each grid. Store
* the grid_info structs in the grid data base.
* Input: name - name of file to scan
* db - the grid data base
* Return: number of grids found.
*/
int get_epa_info( char *name, struct grid_db *db )
{
#ifdef EPA
int fd;
char message[1000];
int nr, nc, nl, numtimes, numvars;
char specstr[1000], datestr[1000];
char varname[MAXVARS][10];
int varnumber[MAXVARS];
int i, j, k;
int var, time;
int datestamp[MAXTIMES], timestamp[MAXTIMES];
int grids;
struct projection *proj;
struct vcs *vcs;
if (!model_open( &fd, name, message)) {
printf("%s\n", message);
return 0;
}
/* printf("get mm4 info\n");*/
if (!model_inquire( fd, message, &nc, &nr, &nl, &numtimes, &numvars,
specstr )) {
printf("%s\n", message);
close(fd);
return 0;
}
/* printf("nr=%d nc=%d nl=%d times=%d vars=%d\n", nr, nc, nl, numtimes,
numvars );
printf("varnames: %s\n", specstr );
*/
/* Get variable names */
k = 0;
for (var=0;var<numvars;var++) {
/* skip garbage chars */
while(!(isalnum(specstr[k]) || specstr[k] == ' ')) k++;
/* copy good chars */
i = 0;
while(isalnum(specstr[k]) || specstr[k] == ' ') {
varname[var][i++] = specstr[k++];
}
varname[var][i] = 0;
/* remove trailing spaces */
/* CAN'T REMOVE TRAILING SPACES: IT CAUSES PROBLEMS IN THE mm4.c FUNCTIONS!
if (i>0) {
i--;
while (varname[var][i]==' ') {
varname[var][i] = 0;
i--;
}
}
*/
/* printf("Varname[%d] = %s.\n", var, varname[var] );*/
varnumber[var] = var;
}
/* Get time/date stamps */
for (time=0;time<numtimes;time++) {
if (!model_get_date( fd, time+1, datestr, message)) {
printf("%s\n", message );
return 0;
}
convert_date( datestr, &datestamp[time], ×tamp[time] );
timestamp[time] = v5dSecondsToHHMMSS( timestamp[time] );
datestamp[time] = v5dDaysToYYDDD( datestamp[time] );
}
/* Get map projection and vertical coord system */
if (!get_epa_projection( db, fd, nr, nc, nl, numvars, varname,
&proj, &vcs )) {
close(fd);
return 0;
}
/* Remove variables names which aren't used after first time step */
/* This is done because the first timestep has some special "meta" */
/* variables which we don't want to visualize. This has to be done */
/* after get_epa_projection() is called! */
#ifdef LEAVEOUT
if (model_type == MM4 && numtimes > 1) {
float *grid3d = (float *) malloc( nr * nc * nl * sizeof(float) );
for (var=0;var<numvars;var++) {
char species[100], units[100];
strcpy( species, varname[var] );
if (!model_get_data( fd, 2, var+1, grid3d, species, units, message)) {
/*printf("%d: %s %s\n", var, species, message );*/
/* remove this variable */
for (j=var;j<numvars-1;j++) {
strcpy( varname[j], varname[j+1] );
varnumber[j] = varnumber[j+1];
}
numvars--;
var--; /* this is tricky */
}
}
free( grid3d );
}
#endif
#ifdef LEAVEOUT
if (model_type==MM4) {
r->Guard = 2;
}
else if (model_type==RADM) {
r->Guard = 1;
}
else {
r->Guard = 1;
}
#endif
/* printf("numtimes=%d numvars=%d\n", numtimes, numvars );*/
/* Add grid_info nodes to grid list */
grids = 0;
for (time=0;time<numtimes;time++) {
for (var=0;var<numvars;var++) {
struct grid_info *info;
info = alloc_grid_info();
info->FileName = strdup( name );
info->Format = FILE_EPA;
info->TimeStep = time;
info->VarNum = varnumber[var];
info->Nr = nr;
info->Nc = nc;
info->Nl = nl;
info->DateStamp = datestamp[time];
info->TimeStamp = timestamp[time];
info->VarName = strdup( varname[var] );
/* map projection & vcs */
info->Proj = proj;
info->Vcs = vcs;
append_grid( info, db );
grids++;
}
}
/* printf("read %d grids\n", grids );*/
return grids;
#else
printf("Warning: Can't read EPA mm4 files on this system\n");
return 0;
#endif
}
/*
* Reverse the order of the rows in a 2-D array.
*/
static void flip_north_south( float data[], int rows, int columns )
{
int i, nbytes;
float temp[MAXCOLUMNS];
nbytes = columns * sizeof(float);
for (i=0;i<rows/2;i++) {
memcpy( temp, data+i*columns, nbytes );
memcpy( data+i*columns, data+(rows-i-1)*columns, nbytes );
memcpy( data+(rows-i-1)*columns, temp, nbytes );
}
}
/*
* Change 2-D array from row-major to column-major.
*/
static void transpose( float in[], float out[], int rows, int columns )
{
int i, j, p;
static float *temp = NULL;
static int tempsize = 0;
/* allocate temporary buffer if needed or if current one is too small */
if (!temp || tempsize < rows*columns) {
if (temp) free(temp);
temp = (float *) malloc( rows * columns * sizeof(float) );
tempsize = rows * columns;
}
p = 0;
for (j=0;j<columns;j++) {
for (i=0;i<rows;i++) {
temp[ j*rows+i ] = in[i*columns+j];
}
}
memcpy( out, temp, rows*columns*sizeof(float) );
}
/*
* Get the grid data described by g.
* Input: g - pointer to a grid_info structure. It tells us which grid
* in which file is needed.
* Return: pointer to grid data (can be free()'d when finished)
* OR return NULL if there's an error.
*/
float *get_epa_data( struct grid_info *g )
{
#ifdef EPA
float *data;
int count;
char species[80], units[80], message[80];
int fd;
int i;
/* open file */
if (!model_open( &fd, g->FileName, message)) {
printf("%s\n", message);
return NULL;
}
/* allocate buffer */
count = g->Nr * g->Nc * g->Nl;
data = (float *) malloc( count * sizeof(float) );
if (!data) {
printf("Error: out of memory in get_epa_data\n");
close(fd);
return NULL;
}
strcpy( species, g->VarName );
/* printf("EPA get data: t=%d v=%d\n", g->TimeStep+1, g->VarNum+1 );*/
if (!model_get_data( fd, g->TimeStep+1, g->VarNum+1, data,
species, units, message)) {
printf("Error in get_epa_data: %s\n", message );
free(data);
close(fd);
return NULL;
}
/*
printf("begin\n");
printf("spec: %s units: %s\n", species, units );
print_min_max( data, g->Nr * g->Nc * g->Nl );
*/
if (model_type==RADM) {
/* multiply all values by 1000.0 */
for (i=0;i<count;i++) {
data[i] *= 1000.0;
}
}
/* flip north/south */
for (i=0;i<g->Nl;i++) {
flip_north_south( data + i*g->Nr*g->Nc, g->Nr, g->Nc );
}
if (model_type == MM4) {
/* flip top/bottom */
int size = g->Nr * g->Nc * sizeof(float);
float *temp = (float *) malloc( size );
for (i=0; i<g->Nl/2 ;i++) {
memcpy( temp, data + i*g->Nr*g->Nc, size );
memcpy( data + i*g->Nr*g->Nc, data + (g->Nl-1-i)*g->Nr*g->Nc, size );
memcpy( data + (g->Nl-1-i)*g->Nr*g->Nc, temp, size );
}
free(temp);
}
/* change from row-major to column-major order */
for (i=0; i<g->Nl; i++) {
transpose( data + i * g->Nr * g->Nc,
data + i * g->Nr * g->Nc,
g->Nr, g->Nc );
}
/* print_min_max( data, g->Nr * g->Nc * g->Nl );*/
close(fd);
return data;
#else
return NULL;
#endif
}
syntax highlighted by Code2HTML, v. 0.9.1