grid.c

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <assert.h>
#include <string.h>
#include <float.h>  /* FLT_EPSILON */
#include <limits.h> /* INT_MAX     */

#include "dmemory.h"
#include "cdi.h"
#include "cdi_cksum.h"
#include "cdi_int.h"
#include "cdi_uuid.h"
#include "grid.h"
#include "gaussgrid.h"
#include "resource_handle.h"
#include "resource_unpack.h"
#include "namespace.h"
#include "serialize.h"
#include "vlist.h"


/* the value in the second pair of brackets must match the length of
 * the longest string (including terminating NUL) */
static const char Grids[][17] = {
  /*  0 */  "undefined",
  /*  1 */  "generic",
  /*  2 */  "gaussian",
  /*  3 */  "gaussian reduced",
  /*  4 */  "lonlat",
  /*  5 */  "spectral",
  /*  6 */  "fourier",
  /*  7 */  "gme",
  /*  8 */  "trajectory",
  /*  9 */  "unstructured",
  /* 10 */  "curvilinear",
  /* 11 */  "lcc",
  /* 12 */  "projection",
};

/* must match table below */
enum xystdname_idx {
  grid_xystdname_grid_latlon,
  grid_xystdname_latlon,
  grid_xystdname_projection,
};
static const char xystdname_tab[][2][24] = {
  [grid_xystdname_grid_latlon] = { "grid_longitude",
                                   "grid_latitude" },
  [grid_xystdname_latlon] = { "longitude",
                              "latitude" },
  [grid_xystdname_projection] = { "projection_x_coordinate",
                                  "projection_y_coordinate" },

};



static int    gridCompareP    ( void * gridptr1, void * gridptr2 );
static void   gridDestroyP    ( void * gridptr );
static void   gridPrintP      ( void * gridptr, FILE * fp );
static int    gridGetPackSize ( void * gridptr, void *context);
static void   gridPack        ( void * gridptr, void * buff, int size, int *position, void *context);
static int    gridTxCode      ( void );

static const resOps gridOps = {
  gridCompareP,
  gridDestroyP,
  gridPrintP,
  gridGetPackSize,
  gridPack,
  gridTxCode
};

static int  GRID_Debug = 0;   /* If set to 1, debugging */


grid_t *grid_to_pointer(int gridID)
{
  return (grid_t *)reshGetVal(gridID, &gridOps);
}

#define gridMark4Update(gridID) reshSetStatus(gridID, &gridOps, RESH_DESYNC_IN_USE)


void grid_init(grid_t *gridptr)
{
  gridptr->self          = CDI_UNDEFID;
  gridptr->type          = CDI_UNDEFID;
  gridptr->proj          = CDI_UNDEFID;
  gridptr->projtype      = CDI_UNDEFID;
  gridptr->mask          = NULL;
  gridptr->mask_gme      = NULL;
  gridptr->x.vals        = NULL;
  gridptr->y.vals        = NULL;
  gridptr->x.bounds      = NULL;
  gridptr->y.bounds      = NULL;
  gridptr->area          = NULL;
  gridptr->rowlon        = NULL;
  gridptr->nrowlon       = 0;

  gridptr->x.first       = 0.0;
  gridptr->x.last        = 0.0;
  gridptr->x.inc         = 0.0;
  gridptr->y.first       = 0.0;
  gridptr->y.last        = 0.0;
  gridptr->y.inc         = 0.0;

  gridptr->lcc.originLon = 0.0;
  gridptr->lcc.originLat = 0.0;
  gridptr->lcc.lonParY   = 0.0;
  gridptr->lcc.lat1      = 0.0;
  gridptr->lcc.lat2      = 0.0;
  gridptr->lcc.xinc      = 0.0;
  gridptr->lcc.yinc      = 0.0;
  gridptr->lcc.projflag  = 0;
  gridptr->lcc.scanflag  = 0;
  gridptr->lcc.defined   = FALSE;

  gridptr->gme.nd        = 0;
  gridptr->gme.ni        = 0;
  gridptr->gme.ni2       = 0;
  gridptr->gme.ni3       = 0;

  gridptr->trunc         = 0;
  gridptr->nvertex       = 0;
  gridptr->number        = 0;
  gridptr->position      = 0;
  gridptr->reference     = NULL;
  gridptr->prec          = 0;
  gridptr->size          = 0;
  gridptr->x.size        = 0;
  gridptr->y.size        = 0;
  gridptr->np            = 0;
  gridptr->x.flag        = 0;
  gridptr->y.flag        = 0;
  gridptr->isCyclic      = CDI_UNDEFID;

  gridptr->lcomplex      = false;
  gridptr->hasdims       = true;
  gridptr->x.dimname[0]  = 0;
  gridptr->y.dimname[0]  = 0;
  gridptr->x.name[0]     = 0;
  gridptr->y.name[0]     = 0;
  gridptr->x.longname[0] = 0;
  gridptr->y.longname[0] = 0;
  gridptr->x.units[0]    = 0;
  gridptr->y.units[0]    = 0;
  gridptr->x.stdname     = NULL;
  gridptr->y.stdname     = NULL;
  gridptr->vdimname[0]   = 0;
  gridptr->mapname[0]    = 0;
  gridptr->mapping[0]    = 0;
  memset(gridptr->uuid, 0, CDI_UUID_SIZE);
  gridptr->name          = NULL;
  gridptr->vtable        = &cdiGridVtable;
  gridptr->atts.nalloc   = MAX_ATTRIBUTES;
  gridptr->atts.nelems   = 0;
}


static
void grid_free_components(grid_t *gridptr)
{
  void *p2free[] = { gridptr->mask, gridptr->mask_gme,
                     gridptr->x.vals, gridptr->y.vals,
                     gridptr->x.bounds, gridptr->y.bounds,
                     gridptr->rowlon, gridptr->area,
                     gridptr->reference, gridptr->name};

  for ( size_t i = 0; i < sizeof(p2free)/sizeof(p2free[0]); ++i )
    if ( p2free[i] ) Free(p2free[i]);
}

void grid_free(grid_t *gridptr)
{
  grid_free_components(gridptr);
  grid_init(gridptr);
}

static
grid_t *gridNewEntry(cdiResH resH)
{
  grid_t *gridptr = (grid_t*) Malloc(sizeof(grid_t));
  grid_init(gridptr);

  if ( resH == CDI_UNDEFID )
    gridptr->self = reshPut(gridptr, &gridOps);
  else
    {
      gridptr->self = resH;
      reshReplace(resH, gridptr, &gridOps);
    }

  return gridptr;
}

static
void gridInit(void)
{
  static bool gridInitialized = false;
  if ( gridInitialized ) return;
  gridInitialized = true;

  const char *env = getenv("GRID_DEBUG");
  if ( env ) GRID_Debug = atoi(env);
}

static
void grid_copy_base_scalar_fields(grid_t *gridptrOrig, grid_t *gridptrDup)
{
  memcpy(gridptrDup, gridptrOrig, sizeof(grid_t));
  gridptrDup->self = CDI_UNDEFID;
  if ( gridptrOrig->reference )
    gridptrDup->reference = strdupx(gridptrOrig->reference);
}


static
grid_t *grid_copy_base(grid_t *gridptrOrig)
{
  grid_t *gridptrDup = (grid_t *)Malloc(sizeof (*gridptrDup));
  gridptrOrig->vtable->copyScalarFields(gridptrOrig, gridptrDup);
  gridptrOrig->vtable->copyArrayFields(gridptrOrig, gridptrDup);
  return gridptrDup;
}


unsigned cdiGridCount(void)
{
  return reshCountType(&gridOps);
}

static inline
void gridSetString(char *gridstrname, const char *name, size_t len)
{
  if ( len > CDI_MAX_NAME ) len = CDI_MAX_NAME;
  strncpy(gridstrname, name, len);
  gridstrname[len - 1] = 0;
}

static inline
void gridGetString(char *name, const char *gridstrname, size_t len)
{
  if ( len > CDI_MAX_NAME ) len = CDI_MAX_NAME;
  strncpy(name, gridstrname, len);
  name[len - 1] = 0;
}

static inline
void gridSetName(char *gridstrname, const char *name)
{
  strncpy(gridstrname, name, CDI_MAX_NAME);
  gridstrname[CDI_MAX_NAME - 1] = 0;
}


void cdiGridTypeInit(grid_t *gridptr, int gridtype, int size)
{
  gridptr->type = gridtype;
  gridptr->size = size;

  if      ( gridtype == GRID_CURVILINEAR  ) gridptr->nvertex = 4;
  else if ( gridtype == GRID_UNSTRUCTURED ) gridptr->x.size = size;

  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
    case GRID_GAUSSIAN_REDUCED:
    case GRID_TRAJECTORY:
    case GRID_CURVILINEAR:
    case GRID_UNSTRUCTURED:
    case GRID_GME:
      {
        if ( gridtype == GRID_TRAJECTORY )
          {
            if ( gridptr->x.name[0] == 0 ) gridSetName(gridptr->x.name, "tlon");
            if ( gridptr->y.name[0] == 0 ) gridSetName(gridptr->y.name, "tlat");
          }
        else
          {
            if ( gridptr->x.name[0] == 0 ) gridSetName(gridptr->x.name, "lon");
            if ( gridptr->y.name[0] == 0 ) gridSetName(gridptr->y.name, "lat");
          }

        gridSetName(gridptr->x.longname, "longitude");
        gridSetName(gridptr->y.longname, "latitude");

        gridptr->x.stdname = xystdname_tab[grid_xystdname_latlon][0];
        gridptr->y.stdname = xystdname_tab[grid_xystdname_latlon][1];
        gridSetName(gridptr->x.units, "degrees_east");
        gridSetName(gridptr->y.units, "degrees_north");

        break;
      }
    case GRID_GENERIC:
    case GRID_PROJECTION:
      {
        if ( gridptr->x.name[0] == 0 ) gridSetName(gridptr->x.name, "x");
        if ( gridptr->y.name[0] == 0 ) gridSetName(gridptr->y.name, "y");
        if ( gridtype == GRID_PROJECTION )
          {
            gridSetName(gridptr->mapname, "Projection");

            gridptr->x.stdname = xystdname_tab[grid_xystdname_projection][0];
            gridptr->y.stdname = xystdname_tab[grid_xystdname_projection][1];
            gridSetName(gridptr->x.units, "m");
            gridSetName(gridptr->y.units, "m");
          }
        break;
      }
    }
}


// used also in CDO
void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *restrict xvals)
{
  if ( (! (fabs(xinc) > 0)) && xsize > 1 )
    {
      if ( xfirst >= xlast )
        {
          while ( xfirst >= xlast ) xlast += 360;
          xinc = (xlast-xfirst)/(xsize);
        }
      else
        {
          xinc = (xlast-xfirst)/(xsize-1);
        }
    }

  for ( int i = 0; i < xsize; ++i )
    xvals[i] = xfirst + i*xinc;
}

static
void calc_gaussgrid(double *restrict yvals, int ysize, double yfirst, double ylast)
{
  double *restrict yw = (double *) Malloc((size_t)ysize * sizeof(double));
  gaussaw(yvals, yw, (size_t)ysize);
  Free(yw);
  for (int i = 0; i < ysize; i++ )
    yvals[i] = asin(yvals[i])/M_PI*180.0;

  if ( yfirst < ylast && yfirst > -90.0 && ylast < 90.0 )
    {
      int yhsize = ysize/2;
      for (int i = 0; i < yhsize; i++ )
        {
          double ytmp = yvals[i];
          yvals[i] = yvals[ysize-i-1];
          yvals[ysize-i-1] = ytmp;
        }
    }
}

// used also in CDO
void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double yinc, double *restrict yvals)
{
  const double deleps = 0.002;

  if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
    {
      if ( ysize > 2 )
	{
	  calc_gaussgrid(yvals, ysize, yfirst, ylast);

	  if ( ! (IS_EQUAL(yfirst, 0) && IS_EQUAL(ylast, 0)) )
	    if ( fabs(yvals[0] - yfirst) > deleps || fabs(yvals[ysize-1] - ylast) > deleps )
	      {
		double *restrict ytmp = NULL;
		int nstart, lfound = 0;
		int ny = (int) (180./(fabs(ylast-yfirst)/(ysize-1)) + 0.5);
		ny -= ny%2;
		if ( ny > ysize && ny < 4096 )
		  {
		    ytmp = (double *) Malloc((size_t)ny * sizeof (double));
		    calc_gaussgrid(ytmp, ny, yfirst, ylast);
                    {
                      int i;
                      for ( i = 0; i < (ny-ysize); i++ )
                        if ( fabs(ytmp[i] - yfirst) < deleps ) break;
                      nstart = i;
                    }

		    lfound = (nstart+ysize-1) < ny
                      && fabs(ytmp[nstart+ysize-1] - ylast) < deleps;
                    if ( lfound )
                      {
                        for (int i = 0; i < ysize; i++) yvals[i] = ytmp[i+nstart];
                      }
		  }

		if ( !lfound )
		  {
		    Warning("Cannot calculate gaussian latitudes for lat1 = %g latn = %g!", yfirst, ylast);
		    for (int i = 0; i < ysize; i++ ) yvals[i] = 0;
		    yvals[0] = yfirst;
		    yvals[ysize-1] = ylast;
		  }

		if ( ytmp ) Free(ytmp);
	      }
	}
      else
        {
          yvals[0] = yfirst;
          yvals[ysize-1] = ylast;
        }
    }
  /*     else if ( gridtype == GRID_LONLAT || gridtype == GRID_GENERIC ) */
  else
    {
      if ( (! (fabs(yinc) > 0)) && ysize > 1 )
        {
          if ( IS_EQUAL(yfirst, ylast) && IS_NOT_EQUAL(yfirst, 0) ) ylast *= -1;

          if ( yfirst > ylast )
            yinc = (yfirst-ylast)/(ysize-1);
          else if ( yfirst < ylast )
            yinc = (ylast-yfirst)/(ysize-1);
          else
            {
              if ( ysize%2 != 0 )
                {
                  yinc = 180.0/(ysize-1);
                  yfirst = -90;
                }
              else
                {
                  yinc = 180.0/ysize;
                  yfirst = -90 + yinc/2;
                }
            }
        }

      if ( yfirst > ylast && yinc > 0 ) yinc = -yinc;

      for (int i = 0; i < ysize; i++ )
        yvals[i] = yfirst + i*yinc;
    }
  /*
    else
    Error("unable to calculate values for %s grid!", gridNamePtr(gridtype));
  */
}

/*
@Function  gridCreate
@Title     Create a horizontal Grid

@Prototype int gridCreate(int gridtype, int size)
@Parameter
    @Item  gridtype  The type of the grid, one of the set of predefined CDI grid types.
                     The valid CDI grid types are @func{GRID_GENERIC}, @func{GRID_GAUSSIAN},
                     @func{GRID_LONLAT}, @func{GRID_LCC}, @func{GRID_SPECTRAL},
                     @func{GRID_GME}, @func{GRID_CURVILINEAR} and @func{GRID_UNSTRUCTURED}.
    @Item  size      Number of gridpoints.

@Description
The function @func{gridCreate} creates a horizontal Grid.

@Result
@func{gridCreate} returns an identifier to the Grid.

@Example
Here is an example using @func{gridCreate} to create a regular lon/lat Grid:

@Source
#include "cdi.h"
   ...
#define  nlon  12
#define  nlat   6
   ...
double lons[nlon] = {0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330};
double lats[nlat] = {-75, -45, -15, 15, 45, 75};
int gridID;
   ...
gridID = gridCreate(GRID_LONLAT, nlon*nlat);
gridDefXsize(gridID, nlon);
gridDefYsize(gridID, nlat);
gridDefXvals(gridID, lons);
gridDefYvals(gridID, lats);
   ...
@EndSource
@EndFunction
*/
int gridCreate(int gridtype, int size)
{
  if ( CDI_Debug ) Message("gridtype=%s  size=%d", gridNamePtr(gridtype), size);

  if ( size < 0 || size > INT_MAX ) Error("Grid size (%d) out of bounds (0 - %d)!", size, INT_MAX);

  gridInit();

  grid_t *gridptr = gridNewEntry(CDI_UNDEFID);
  if ( ! gridptr ) Error("No memory");

  int gridID = gridptr->self;

  if ( CDI_Debug ) Message("gridID: %d", gridID);

  cdiGridTypeInit(gridptr, gridtype, size);

  return gridID;
}

static
void gridDestroyKernel( grid_t * gridptr )
{
  xassert ( gridptr );

  int id = gridptr->self;

  grid_free_components(gridptr);
  Free( gridptr );

  reshRemove ( id, &gridOps );
}

/*
@Function  gridDestroy
@Title     Destroy a horizontal Grid

@Prototype void gridDestroy(int gridID)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.

@EndFunction
*/
void gridDestroy(int gridID)
{
  grid_t *gridptr = grid_to_pointer(gridID);
  gridptr->vtable->destroy(gridptr);
}

static
void gridDestroyP(void * gridptr)
{
  ((grid_t *)gridptr)->vtable->destroy((grid_t *)gridptr);
}


const char *gridNamePtr(int gridtype)
{
  int size = (int) (sizeof(Grids)/sizeof(Grids[0]));

  const char *name = (gridtype >= 0 && gridtype < size) ? Grids[gridtype] : Grids[GRID_GENERIC];

  return name;
}


void gridName(int gridtype, char *gridname)
{
  strcpy(gridname, gridNamePtr(gridtype));
}

static
void *grid_key_to_ptr(grid_t *gridptr, int key)
{
  void *keyptr = NULL;

  switch (key)
    {
    case CDI_KEY_XNAME:      keyptr = (void*)gridptr->x.name; break;
    case CDI_KEY_XLONGNAME:  keyptr = (void*)gridptr->x.longname; break;
    case CDI_KEY_XUNITS:     keyptr = (void*)gridptr->x.units; break;
    case CDI_KEY_YNAME:      keyptr = (void*)gridptr->y.name; break;
    case CDI_KEY_YLONGNAME:  keyptr = (void*)gridptr->y.longname; break;
    case CDI_KEY_YUNITS:     keyptr = (void*)gridptr->y.units; break;
    case CDI_KEY_XDIMNAME:   keyptr = (void*)gridptr->x.dimname; break;
    case CDI_KEY_YDIMNAME:   keyptr = (void*)gridptr->y.dimname; break;
    case CDI_KEY_VDIMNAME:   keyptr = (void*)gridptr->vdimname; break;
    case CDI_KEY_MAPNAME:    keyptr = (void*)gridptr->mapname; break;
    case CDI_KEY_MAPPING:    keyptr = (void*)gridptr->mapping; break;
    }

  return keyptr;
}

/*
@Function  cdiGridDefKeyStr
@Title     Define a CDI grid string value from a key

@Prototype int cdiGridDefKeyStr(int gridID, int key, int size, const char *mesg)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  key      The key to be searched
    @Item  size     The allocated length of the string on input
    @Item  mesg     The address of a string where the data will be read

@Description
The function @func{cdiGridDefKeyStr} defines a CDI grid string value from a key.

@Result
@func{cdiGridDefKeyStr} returns 0 if OK and integer value on error.

@EndFunction
*/
int cdiGridDefKeyStr(int gridID, int key, int size, const char *mesg)
{
  if ( size < 1 || mesg == NULL || *mesg == 0 ) return -1;

  grid_t *gridptr = grid_to_pointer(gridID);

  char *keyptr = (char*)grid_key_to_ptr(gridptr, key);
  if ( keyptr == NULL )
    {
      Warning("CDI grid string key %d not supported!", key);
      return -1;
    }

  gridSetString(keyptr, mesg, (size_t)size);
  gridMark4Update(gridID);

  return 0;
}

/*
@Function  cdiGridInqKeyStr
@Title     Get a CDI grid string value from a key

@Prototype int cdiGridInqKeyStr(int gridID, int key, int size, char *mesg)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  key      The key to be searched.
    @Item  size     The allocated length of the string on input.
    @Item  mesg     The address of a string where the data will be retrieved.
                    The caller must allocate space for the returned string.
                    The maximum possible length, in characters, of the string
                    is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{cdiGridInqKeyStr} return a CDI grid string value from a key.

@Result
@func{cdiGridInqKeyStr} returns 0 if OK and integer value on error.

@EndFunction
*/
int cdiGridInqKeyStr(int gridID, int key, int size, char *mesg)
{
  if ( size < 1 || mesg == NULL ) return -1;

  grid_t *gridptr = grid_to_pointer(gridID);
  const char *keyptr = (const char*)grid_key_to_ptr(gridptr, key);
  if ( keyptr == NULL)
    {
      Warning("CDI grid string key %d not supported!", key);
      return -1;
    }

  gridGetString(mesg, keyptr, (size_t)size);

  return 0;
}

/*
@Function  gridDefXname
@Title     Define the name of a X-axis

@Prototype void gridDefXname(int gridID, const char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  name     Name of the X-axis.

@Description
The function @func{gridDefXname} defines the name of a X-axis.

@EndFunction
*/
void gridDefXname(int gridID, const char *xname)
{
  (void)cdiGridDefKeyStr(gridID, CDI_KEY_XNAME, CDI_MAX_NAME, xname);
}

/*
@Function  gridDefXlongname
@Title     Define the longname of a X-axis

@Prototype void gridDefXlongname(int gridID, const char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  longname Longname of the X-axis.

@Description
The function @func{gridDefXlongname} defines the longname of a X-axis.

@EndFunction
*/
void gridDefXlongname(int gridID, const char *xlongname)
{
  (void)cdiGridDefKeyStr(gridID, CDI_KEY_XLONGNAME, CDI_MAX_NAME, xlongname);
}

/*
@Function  gridDefXunits
@Title     Define the units of a X-axis

@Prototype void gridDefXunits(int gridID, const char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  units    Units of the X-axis.

@Description
The function @func{gridDefXunits} defines the units of a X-axis.

@EndFunction
*/
void gridDefXunits(int gridID, const char *xunits)
{
  (void)cdiGridDefKeyStr(gridID, CDI_KEY_XUNITS, CDI_MAX_NAME, xunits);
}

/*
@Function  gridDefYname
@Title     Define the name of a Y-axis

@Prototype void gridDefYname(int gridID, const char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  name     Name of the Y-axis.

@Description
The function @func{gridDefYname} defines the name of a Y-axis.

@EndFunction
*/
void gridDefYname(int gridID, const char *yname)
{
  (void)cdiGridDefKeyStr(gridID, CDI_KEY_YNAME, CDI_MAX_NAME, yname);
}

/*
@Function  gridDefYlongname
@Title     Define the longname of a Y-axis

@Prototype void gridDefYlongname(int gridID, const char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  longname Longname of the Y-axis.

@Description
The function @func{gridDefYlongname} defines the longname of a Y-axis.

@EndFunction
*/
void gridDefYlongname(int gridID, const char *ylongname)
{
  (void)cdiGridDefKeyStr(gridID, CDI_KEY_YLONGNAME, CDI_MAX_NAME, ylongname);
}

/*
@Function  gridDefYunits
@Title     Define the units of a Y-axis

@Prototype void gridDefYunits(int gridID, const char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  units    Units of the Y-axis.

@Description
The function @func{gridDefYunits} defines the units of a Y-axis.

@EndFunction
*/
void gridDefYunits(int gridID, const char *yunits)
{
  (void)cdiGridDefKeyStr(gridID, CDI_KEY_YUNITS, CDI_MAX_NAME, yunits);
}

/*
@Function  gridInqXname
@Title     Get the name of a X-axis

@Prototype void gridInqXname(int gridID, char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  name     Name of the X-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqXname} returns the name of a X-axis.

@Result
@func{gridInqXname} returns the name of the X-axis to the parameter name.

@EndFunction
*/
void gridInqXname(int gridID, char *xname)
{
  (void)cdiGridInqKeyStr(gridID, CDI_KEY_XNAME, CDI_MAX_NAME, xname);
}

/*
@Function  gridInqXlongname
@Title     Get the longname of a X-axis

@Prototype void gridInqXlongname(int gridID, char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  longname Longname of the X-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqXlongname} returns the longname of a X-axis.

@Result
@func{gridInqXlongname} returns the longname of the X-axis to the parameter longname.

@EndFunction
*/
void gridInqXlongname(int gridID, char *xlongname)
{
  (void)cdiGridInqKeyStr(gridID, CDI_KEY_XLONGNAME, CDI_MAX_NAME, xlongname);
}

/*
@Function  gridInqXunits
@Title     Get the units of a X-axis

@Prototype void gridInqXunits(int gridID, char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  units    Units of the X-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqXunits} returns the units of a X-axis.

@Result
@func{gridInqXunits} returns the units of the X-axis to the parameter units.

@EndFunction
*/
void gridInqXunits(int gridID, char *xunits)
{
  (void)cdiGridInqKeyStr(gridID, CDI_KEY_XUNITS, CDI_MAX_NAME, xunits);
}


void gridInqXstdname(int gridID, char *xstdname)
{
  if ( xstdname )
    {
      xstdname[0] = 0;
      grid_t *gridptr = grid_to_pointer(gridID);
      if ( gridptr->x.stdname ) strcpy(xstdname, gridptr->x.stdname);
    }
}

/*
@Function  gridInqYname
@Title     Get the name of a Y-axis

@Prototype void gridInqYname(int gridID, char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  name     Name of the Y-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqYname} returns the name of a Y-axis.

@Result
@func{gridInqYname} returns the name of the Y-axis to the parameter name.

@EndFunction
*/
void gridInqYname(int gridID, char *yname)
{
  (void)cdiGridInqKeyStr(gridID, CDI_KEY_YNAME, CDI_MAX_NAME, yname);
}

/*
@Function  gridInqYlongname
@Title     Get the longname of a Y-axis

@Prototype void gridInqYlongname(int gridID, char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  longname Longname of the Y-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqYlongname} returns the longname of a Y-axis.

@Result
@func{gridInqYlongname} returns the longname of the Y-axis to the parameter longname.

@EndFunction
*/
void gridInqYlongname(int gridID, char *ylongname)
{
  (void)cdiGridInqKeyStr(gridID, CDI_KEY_YLONGNAME, CDI_MAX_NAME, ylongname);
}

/*
@Function  gridInqYunits
@Title     Get the units of a Y-axis

@Prototype void gridInqYunits(int gridID, char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
    @Item  units    Units of the Y-axis. The caller must allocate space for the
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqYunits} returns the units of a Y-axis.

@Result
@func{gridInqYunits} returns the units of the Y-axis to the parameter units.

@EndFunction
*/
void gridInqYunits(int gridID, char *yunits)
{
  (void)cdiGridInqKeyStr(gridID, CDI_KEY_YUNITS, CDI_MAX_NAME, yunits);
}


void gridInqYstdname(int gridID, char *ystdname)
{
  if ( ystdname )
    {
      ystdname[0] = 0;
      grid_t *gridptr = grid_to_pointer(gridID);
      if ( gridptr->y.stdname ) strcpy(ystdname, gridptr->y.stdname);
    }
}


void gridDefProj(int gridID, int projID)
{
  grid_t *gridptr = grid_to_pointer(gridID);
  gridptr->proj = projID;

  if ( gridptr->type == GRID_CURVILINEAR )
    {
      grid_t *projptr = grid_to_pointer(projID);
      if ( projptr->x.name[0] ) strcpy(gridptr->x.dimname, projptr->x.name);
      if ( projptr->y.name[0] ) strcpy(gridptr->y.dimname, projptr->y.name);
    }
}


int gridInqProj(int gridID)
{
  grid_t *gridptr = grid_to_pointer(gridID);
  return gridptr->proj;
}


int gridInqProjType(int gridID)
{
  grid_t *gridptr = grid_to_pointer(gridID);

  int projtype = gridptr->projtype;

  if ( projtype == -1 )
    {
      char mapping[CDI_MAX_NAME]; mapping[0] = 0;
      cdiGridInqKeyStr(gridID, CDI_KEY_MAPPING, CDI_MAX_NAME, mapping);
      if ( mapping[0] )
        {
          if      ( strcmp(mapping, "rotated_latitude_longitude") == 0 )   projtype = CDI_PROJ_RLL;
          else if ( strcmp(mapping, "lambert_azimuthal_equal_area") == 0 ) projtype = CDI_PROJ_LAEA;
          else if ( strcmp(mapping, "lambert_conformal_conic") == 0 )      projtype = CDI_PROJ_LCC;
          else if ( strcmp(mapping, "sinusoidal") == 0 )                   projtype = CDI_PROJ_SINU;

          gridptr->projtype = projtype;
        }
    }

  return projtype;
}