/* projlist.c */
/*
* 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"
/*
* Manage a list of map projections and vertical coordinate systems (vcs).
* Basically each grid_info struct (see grid.h) will point to a
* map projection and vcs structure kept in a list here.
*/
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "grid_i.h"
#include "memory.h"
#include "projlist_i.h"
#include "v5d.h"
#ifndef M_PI
# define M_PI 3.1415926
#endif
#define DEG2RAD (M_PI/180.0)
#define EARTH_RADIUS 6371.23
/*
* Test if two floats are nearly equal, i.e. within some epsilon of
* each other.
*/
static int equal( float x, float y )
{
float diff = x - y;
if (diff<0.001 && diff>-0.001) {
return 1;
}
else {
return 0;
}
}
/* Return the sign of x */
static float sign( float x )
{
if (x<0.0) {
return -1.0;
}
else if (x>0.0) {
return 1.0;
}
else {
return 0.0;
}
}
/*
* When doing (row,column) to (lat,lon) conversions it's sometimes useful
* to have some auxillary values around to simplify the calculations. This
* function computes those auxillary values, if any, for a projection and
* stores them in the projection's auxargs array.
* Input: proj - the projection
*/
static void compute_aux_proj_args( struct projection *proj )
{
float lat1, lat2;
switch (proj->Kind) {
case PROJ_LAMBERT: /* Lambert conformal */
#define Lat1 proj->Args[0]
#define Lat2 proj->Args[1]
#define PoleRow proj->Args[2]
#define PoleCol proj->Args[3]
#define CentralLon proj->Args[4]
#define ColInc proj->Args[5]
#define Hemisphere proj->AuxArgs[0]
#define ConeFactor proj->AuxArgs[1]
#define Cone proj->AuxArgs[2]
proj->AuxArgs = (float *) MALLOC( 3 * sizeof(float) );
if (Lat1==Lat2) {
/* polar stereographic? */
if (Lat1>0.0) {
lat1 = (90.0 - Lat1) * DEG2RAD;
/* Cone = 1.0; */
}
else {
lat1 = (90.0 + Lat1) * DEG2RAD;
/* Cone = -1.0; */
}
Cone = cos(lat1); /* WLH 9-27-96 */
Hemisphere = 1.0;
}
else {
/* general Lambert conformal */
float a, b;
if (sign(Lat1) != sign(Lat2)) {
printf("Error: standard latitudes must have the same sign.\n");
exit(1);
}
if (Lat1<Lat2) {
printf("Error: Lat1 must be >= Lat2\n");
exit(1);
}
Hemisphere = 1.0;
lat1 = (90.0 - Lat1) * DEG2RAD;
lat2 = (90.0 - Lat2) * DEG2RAD;
a = log( sin(lat1) ) - log( sin(lat2) );
b = log( tan(lat1/2.0) ) - log( tan(lat2/2.0) );
Cone = a / b;
}
/* Cone is in [-1,1] */
ConeFactor = EARTH_RADIUS * sin(lat1)
/ (ColInc * Cone * pow(tan(lat1/2.0), Cone) );
#undef Lat1
#undef Lat2
#undef PoleRow
#undef PoleCol
#undef CentralLon
#undef ColInc
#undef Hemisphere
#undef ConeFactor
#undef Cone
break;
default:
/* No auxillary args */
proj->AuxArgs = NULL;
}
}
/*
* Return a pointer to a map projection structure given the arguments
* below. If this is a new projection we'll allocate a new projection
* struct and return a pointer to it. If this projection is already in
* the list, just return a pointer to it.
* Input: kind - the type of projection, one of PROJ_*
* nr, nc - number of rows and columns of data
* args - array of projection parameters
* Return: pointer to a projection struct.
*/
struct projection *new_projection( struct grid_db *db, int kind,
int nr, int nc, float *args )
{
int p, i, nargs;
/* determine how many arguments are in the args array */
switch (kind) {
case PROJ_GENERIC: nargs = 4; break;
case PROJ_LINEAR: nargs = 4; break;
case PROJ_LAMBERT: nargs = 6; break;
case PROJ_STEREO: nargs = 5; break;
case PROJ_ROTATED: nargs = 7; break;
case PROJ_EPA: nargs = nr*nc*2; break;
case PROJ_CYLINDRICAL: nargs = 4; break;
case PROJ_SPHERICAL: nargs = 4; break;
case PROJ_MERCATOR: nargs = 4; break;
default:
printf("Fatal error in new_projection!\n");
exit(-1);
}
/* Search projection list for a possible match */
for (p=0; p<db->NumProj; p++) {
if ( db->ProjList[p]->Kind==kind
&& db->ProjList[p]->Nr==nr
&& db->ProjList[p]->Nc==nc) {
int same = 1;
for (i=0;i<nargs;i++) {
if ( !equal(args[i], db->ProjList[p]->Args[i]) ) {
same = 0;
break;
}
}
if (same) {
return db->ProjList[p];
}
}
}
/* if we get here, the projection is not in the list, make a new one */
if (db->NumProj<IMAXPROJ) {
struct projection *newp;
newp = (struct projection *) calloc( 1, sizeof(struct projection) );
newp->Kind = kind;
newp->Nr = nr;
newp->Nc = nc;
newp->Args = (float *) MALLOC( nargs * sizeof(float) );
for (i=0;i<nargs;i++) {
newp->Args[i] = args[i];
}
/* compute extra, optional proj args */
compute_aux_proj_args( newp );
/* add to end of list */
db->ProjList[db->NumProj] = newp;
db->NumProj++;
return newp;
}
else {
printf("Error: too many map projections, %d is limit\n", IMAXPROJ );
return NULL;
}
}
/*
* Deallocate a map projection structure.
*/
void free_projection( struct grid_db *db, struct projection *proj )
{
int i, j;
assert( db );
assert( proj );
for (i=0; i<db->NumProj; i++) {
if (db->ProjList[i]==proj) {
/* found it, remove from list */
for (j=i; j<db->NumProj-1; j++) {
db->ProjList[j] = db->ProjList[j+1];
}
db->NumProj--;
break;
}
}
free( proj->Args );
free( proj );
}
/*
* Given a pointer to a map projection struct, return the numeric position of
* it in the linked list. The first projection being number 1.
* Input: proj - pointer to a map projection struct
* Return: position of proj in list starting at one, else 0 if not in list
*/
int lookup_proj( struct grid_db *db, struct projection *proj )
{
int i;
for (i=0;i<db->NumProj;i++) {
if (db->ProjList[i]==proj) {
return i+1;
}
}
return 0;
}
#define MAX_PROJ_CHARS 1000
char **sprint_projection_list( struct grid_db *db )
{
struct projection *p;
int i;
char **list;
/* construct array of pointers to strings */
if (db->NumProj==0) {
return NULL;
}
list = (char **) MALLOC( db->NumProj * sizeof(char *) );
for (i=0; i<db->NumProj; i++) {
p = db->ProjList[i];
list[i] = (char*) MALLOC( MAX_PROJ_CHARS );
switch (p->Kind) {
case PROJ_GENERIC:
sprintf( list[i], "%3d Generic Linear %4d %4d %g %g %g %g",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
case PROJ_LINEAR:
sprintf( list[i], "%3d Cyl. Equidistant %4d %4d %g %g %g %g",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
case PROJ_LAMBERT:
sprintf( list[i],
"%3d Lambert Conformal %4d %4d %g %g %g %g %g %g",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3], p->Args[4], p->Args[5] );
break;
case PROJ_MERCATOR:
sprintf( list[i],
"%3d Mercator %4d %4d %g %g %g %g",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
case PROJ_STEREO:
sprintf( list[i], "%3d Stereographic %4d %4d %g %g %g %g %g",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3], p->Args[4] );
break;
case PROJ_ROTATED:
sprintf( list[i],
"%3d Rotated %4d %4d %g %g %g %g %g %g %g",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3], p->Args[4], p->Args[5], p->Args[6] );
break;
case PROJ_EPA:
sprintf( list[i], "%3d EPA %4d %4d",
i+1, p->Nr, p->Nc );
break;
case PROJ_CYLINDRICAL:
sprintf( list[i], "%3d Cylindrical projection %4d %4d %g %g %g %g",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
case PROJ_SPHERICAL:
sprintf( list[i], "%3d Spherical projection %4d %4d %g %g %g %g",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
default:
sprintf( list[i], "Error!" );
}
}
return list;
}
void print_projection_list( struct grid_db *db )
{
struct projection *p;
int i;
/* construct array of pointers to strings */
for (i=0; i<db->NumProj; i++) {
p = db->ProjList[i];
if (db->ProjSelected[i]) {
printf("* ");
}
else {
printf(" ");
}
switch (p->Kind) {
case PROJ_GENERIC:
printf( "%3d Generic Linear %4d %4d %g %g %g %g\n",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
case PROJ_LINEAR:
printf( "%3d Cyl. Equidistant %4d %4d %g %g %g %g\n",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
case PROJ_LAMBERT:
printf( "%3d Lambert Conformal %4d %4d %g %g %g %g %g %g\n",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3], p->Args[4], p->Args[5] );
break;
case PROJ_MERCATOR:
printf( "%3d Mercator %4d %4d %g %g %g %g\n",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
case PROJ_STEREO:
printf( "%3d Stereographic %4d %4d %g %g %g %g %g\n",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3], p->Args[4] );
break;
case PROJ_ROTATED:
printf( "%3d Rotated %4d %4d %g %g %g %g %g %g %g\n",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3], p->Args[4], p->Args[5], p->Args[6] );
break;
case PROJ_EPA:
printf( "%3d EPA %4d %4d\n",
i+1, p->Nr, p->Nc );
break;
case PROJ_CYLINDRICAL:
printf( "%3d Cylindrical projection %4d %4d %g %g %g %g\n",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
case PROJ_SPHERICAL:
printf( "%3d Spherical projection %4d %4d %g %g %g %g\n",
i+1, p->Nr, p->Nc, p->Args[0], p->Args[1], p->Args[2],
p->Args[3] );
break;
default:
assert(1==0);
}
}
}
#ifdef LEAVEOUT
void print_projection_list( db )
struct grid_db *db;
{
int i;
struct projection *p;
printf("Projection list:\n");
for (i=0;i<db->NumProj;i++) {
p = db->ProjList[i];
switch (p->Kind) {
case PROJ_GENERIC: printf("Generic "); break;
case PROJ_LINEAR: printf("Linear "); break;
case PROJ_LAMBERT: printf("Lambert "); break;
case PROJ_STEREO: printf("Stereo "); break;
case PROJ_EPA: printf("Epa "); break;
case PROJ_CYLINDRICAL: printf("Cylindrical"); break;
case PROJ_SPHERICAL: printf("Spherical "); break;
case PROJ_MERCATOR: printf("Mercator "); break;
default: printf("Error! ");
}
printf("%d x %d ", p->Nr, p->Nc );
printf("%g %g %g\n", p->Args[0], p->Args[1], p->Args[2] );
}
}
#endif
/*
* Return a pointer to a vcs structure given the arguments below. If
* this is a new projection we'll allocate a new vcs struct and return a
* pointer to it. If this vertical coordinate system is already in the
* list, just return a pointer to it.
* Input: kind - the type of vcs, one of VERT_*
* nl - total number of levels
* lowlev - location of first grid level, below is missing data
* Note: nl-lowlev = actual number of grid levels
* args - array of vcs parameters
* Return: pointer to a vcs struct.
*/
struct vcs *new_vcs( struct grid_db *db, int kind, int nl, int lowlev,
float *args )
{
int v, i, nargs;
assert(db);
assert(args);
/* determine how many arguments are in the args array */
switch (kind) {
case VERT_GENERIC: nargs = 2; break;
case VERT_EQUAL_KM: nargs = 2; break;
case VERT_UNEQUAL_KM: nargs = nl+lowlev; break;
#if V5D_VERSION >= 42
case VERT_UNEQUAL_MB: nargs = nl+lowlev; break;
#endif
case VERT_EPA: nargs = nl; break;
default:
printf("Fatal error in new_vcs!\n");
exit(-1);
}
/* verify arguments are ok */
if (kind==VERT_UNEQUAL_KM) {
for (i=1;i<nargs;i++) {
if (args[i]<=args[i-1]) {
printf("Error in VCS, heights should increase:");
printf(" hgt[%d]=%g hgt[%d]=%g\n", i-1, args[i-1], i, args[i] );
return NULL;
}
}
}
if (kind==VERT_UNEQUAL_MB) {
for (i=1;i<nargs;i++) {
if (args[i]<=args[i-1]) {
printf("Error in VCS, pressures should decrease:");
printf(" hgt[%d]=%g hgt[%d]=%g\n", i-1,
height_to_pressure(args[i-1]), i,
height_to_pressure(args[i]) );
return NULL;
}
}
}
else if (kind==VERT_EQUAL_KM) {
if (args[1]<0.0) {
printf("Error in VCS, increment can't be negative: %g\n", args[1]);
return NULL;
}
}
/* Search list for this vcs */
for (v=0; v<db->NumVcs; v++) {
if ( db->VcsList[v]->Kind==kind
&& db->VcsList[v]->Nl==nl
&& db->VcsList[v]->LowLev==lowlev) {
int same = 1;
for (i=0;i<nargs;i++) {
if (!equal(args[i],db->VcsList[v]->Args[i])) {
same = 0;
break;
}
}
if (same) {
return db->VcsList[v];
}
}
}
/* if we get here, the vcs is not in the list, make a new one */
if (db->NumVcs<IMAXPROJ) {
struct vcs *newv;
newv = (struct vcs *) calloc( 1, sizeof(struct vcs) );
newv->Kind = kind;
newv->Nl = nl+lowlev;
newv->LowLev = lowlev;
newv->Args = (float *) MALLOC( nargs * sizeof(float) );
for (i=0;i<nargs;i++) {
newv->Args[i] = args[i];
}
/* add to end of list */
db->VcsList[db->NumVcs] = newv;
db->NumVcs++;
return newv;
}
else {
printf("Error: too many vertical coordinate systems, %d is limit\n",
IMAXPROJ );
return NULL;
}
}
/*
* Deallocate a VCS struct.
*/
void free_vcs( struct grid_db *db, struct vcs *vcs )
{
int i, j;
assert( db );
assert( vcs );
for (i=0; i<db->NumVcs; i++) {
if (db->VcsList[i]==vcs) {
/* found it, remove from list */
for (j=i; j<db->NumVcs-1; j++) {
db->VcsList[j] = db->VcsList[j+1];
}
db->NumVcs--;
break;
}
}
free( vcs->Args );
free( vcs );
}
char **sprint_vcs_list( struct grid_db *db )
{
struct vcs *v;
int i;
char **list;
/* construct array of pointers to strings */
if (db->NumVcs==0) {
return NULL;
}
list = (char **) MALLOC( db->NumVcs * sizeof(char *) );
for (i=0; i<db->NumVcs; i++) {
v = db->VcsList[i];
list[i] = (char*) MALLOC( MAX_PROJ_CHARS );
switch (v->Kind) {
case VERT_GENERIC:
sprintf( list[i], "%3d Generic Linear %4d %g %g",
i+1, v->Nl, v->Args[0], v->Args[1] );
break;
case VERT_EQUAL_KM:
sprintf( list[i], "%3d Equally-spaced Linear km %4d %g %g",
i+1, v->Nl, v->Args[0], v->Args[1] );
break;
case VERT_UNEQUAL_KM:
sprintf( list[i], "%3d Unequally-spaced Linear km %4d %g %g",
i+1, v->Nl, v->Args[0], v->Args[1] );
break;
case VERT_UNEQUAL_MB:
sprintf( list[i], "%3d Unequally-spaced Pressure mb %4d %g %g",
i+1, v->Nl, height_to_pressure(v->Args[0]),
height_to_pressure(v->Args[1]) );
break;
default:
sprintf( list[i], "Error!" );
}
}
return list;
}
void print_vcs_list( struct grid_db *db )
{
struct vcs *v;
int i,j;
/* construct array of pointers to strings */
for (i=0; i<db->NumVcs; i++) {
v = db->VcsList[i];
if (db->VcsSelected[i]) {
printf("* ");
}
else {
printf(" ");
}
switch (v->Kind) {
case VERT_GENERIC:
printf( "%3d Generic Linear %4d %g %g\n",
i+1, v->Nl, v->Args[0], v->Args[1] );
break;
case VERT_EQUAL_KM:
printf( "%3d Equally-spaced Linear km %4d %g %g\n",
i+1, v->Nl, v->Args[0], v->Args[1] );
break;
case VERT_UNEQUAL_KM:
printf( "%3d Unequally-spaced Linear km %4d %g %g\n",
i+1, v->Nl, v->Args[0], v->Args[1] );
break;
case VERT_UNEQUAL_MB:
printf( "%3d Unequally-spaced Pressure mb %4d \n",
i+1, v->Nl);
for(j=0;j<v->Nl;j++){
printf(" %3d %6g mb\n",j+1,height_to_pressure(v->Args[j]));
}
break;
case VERT_EPA:
printf( "%3d EPA %4d\n",
i+1, v->Nl );
break;
default:
assert(1==0);
}
}
}
/*
* Given a pointer to a vcs struct, return the numeric position of
* it in the linked list. The first vcs being number 1.
* Input: db - the grid data base
* proj - pointer to a vcs struct
* Return: position of vcs in list starting at one, else 0 if not in list
*/
int lookup_vcs( struct grid_db *db, struct vcs *vcs )
{
int i;
for (i=0; i<db->NumVcs; i++) {
if (db->VcsList[i]==vcs) {
return i+1;
}
}
return 0;
}
/*#define MAXLEVELS 100*/
/*
* Combine all 1-level VCSs into one new coordinate system
*/
struct vcs *combine_vcs( struct grid_db *db, int kind )
{
struct vcs *v;
int count;
float height[MAXLEVELS];
struct vcs *vcsentry[MAXLEVELS];
int i, j;
count = 0;
for (i=0; i<db->NumVcs; i++) {
v = db->VcsList[i];
if (v->Kind==kind && v->Nl==1) {
height[count] = v->Args[0];
vcsentry[count] = v;
count++;
}
}
printf("level height\n");
for (i=0;i<count;i++) {
printf("%3d %7g\n", i, height[i] );
}
/* bubble sort the heights */
for (i=0;i<count-1;i++) {
for (j=i+1;j<count;j++) {
if (height[i]>height[j]) {
float tmp;
struct vcs *u;
tmp = height[i];
height[i] = height[j];
height[j] = tmp;
u = vcsentry[i];
vcsentry[i] = vcsentry[j];
vcsentry[j] = u;
}
}
}
printf("level height\n");
for (i=0;i<count;i++) {
printf("%3d %7g\n", i, height[i] );
}
if (kind==VERT_GENERIC) {
return new_vcs( db, VERT_GENERIC, count, 0, height );
}
else if (kind==VERT_EQUAL_KM || kind==VERT_UNEQUAL_KM) {
return new_vcs( db, VERT_UNEQUAL_KM, count, 0, height );
}
#if V5D_VERSION >= 42
else if (kind==VERT_UNEQUAL_MB) {
return new_vcs( db, VERT_UNEQUAL_MB, count, 0, height );
}
#endif
else {
printf("problem in combine_vcs()!\n");
return NULL;
}
}
syntax highlighted by Code2HTML, v. 0.9.1