IO/NetCDF/vtkNetCDFCFReader.cxx

// -*- c++ -*-
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkNetCDFCFReader.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

/*-------------------------------------------------------------------------
  Copyright 2008 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
  the U.S. Government retains certain rights in this software.
-------------------------------------------------------------------------*/

#include "vtkNetCDFCFReader.h"

#include "vtkCellArray.h"
#include "vtkDataArraySelection.h"
#include "vtkDoubleArray.h"
#include "vtkExtentTranslator.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
#include "vtkMath.h"
#include "vtkMergePoints.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPoints.h"
#include "vtkRectilinearGrid.h"
#include "vtkStdString.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStringArray.h"
#include "vtkStructuredGrid.h"
#include "vtkUnstructuredGrid.h"

#include "vtkSmartPointer.h"
#define VTK_CREATE(type, name) vtkSmartPointer<type> name = vtkSmartPointer<type>::New()

#include <set>
#include <vtksys/RegularExpression.hxx>
#include <vtksys/SystemTools.hxx>

#include <cstring>

#include "vtk_netcdf.h"

#define CALL_NETCDF_GENERIC(call, on_error)                                                        \
  {                                                                                                \
    int errorcode = call;                                                                          \
    if (errorcode != NC_NOERR)                                                                     \
    {                                                                                              \
      const char* errorstring = nc_strerror(errorcode);                                            \
      on_error;                                                                                    \
    }                                                                                              \
  }

#define CALL_NETCDF(call)                                                                          \
  CALL_NETCDF_GENERIC(call, vtkErrorMacro(<< "netCDF Error: " << errorstring); return 0;)

#define CALL_NETCDF_GW(call)                                                                       \
  CALL_NETCDF_GENERIC(call, vtkGenericWarningMacro(<< "netCDF Error: " << errorstring); return 0;)

#include <algorithm>

#include <cmath>

//=============================================================================
// Convenience function for getting the text attribute on a variable.  Returns
// true if the attribute exists, false otherwise.
static bool ReadTextAttribute(int ncFD, int varId, const char* name, vtkStdString& result)
{
  size_t length;
  if (nc_inq_attlen(ncFD, varId, name, &length) != NC_NOERR)
  {
    return false;
  }

  result.resize(length);
  if (length > 0)
  {
    if (nc_get_att_text(ncFD, varId, name, &result.at(0)) != NC_NOERR)
    {
      return false;
    }
  }
  else
  {
    // If length == 0, then there really is nothing to read. Do nothing
  }

  // The line below seems weird, but it is here for a good reason.  In general,
  // text attributes are not null terminated, so you have to add your own (which
  // the vtkStdString will do for us).  However, sometimes a null terminating
  // character is written in the attribute anyway.  In a C string this is no big
  // deal.  But it means that the vtkStdString has a null character in it and it
  // is technically different than its own C string.  This line corrects that
  // regardless of whether the null string was written we will get the right
  // string.
  result = result.c_str();

  return true;
}

//------------------------------------------------------------------------------
// Convenience function for getting the range of all values in all components of
// a vtkDoubleArray.
static void GetRangeOfAllComponents(vtkDoubleArray* array, double range[2])
{
  range[0] = VTK_DOUBLE_MAX;
  range[1] = VTK_DOUBLE_MIN;
  for (int component = 0; component < array->GetNumberOfComponents(); component++)
  {
    double componentRange[2];
    array->GetRange(componentRange, component);
    if (componentRange[0] < range[0])
      range[0] = componentRange[0];
    if (componentRange[1] > range[1])
      range[1] = componentRange[1];
  }
}

//=============================================================================
vtkNetCDFCFReader::vtkDimensionInfo::vtkDimensionInfo(int ncFD, int id)
{
  this->DimId = id;

  this->Units = UNDEFINED_UNITS;
  this->HasRegularSpacing = true;
  this->Origin = 0.0;
  this->Spacing = 1.0;

  this->LoadMetaData(ncFD);
}

int vtkNetCDFCFReader::vtkDimensionInfo::LoadMetaData(int ncFD)
{
  this->Units = UNDEFINED_UNITS;

  char name[NC_MAX_NAME + 1];
  CALL_NETCDF_GW(nc_inq_dimname(ncFD, this->DimId, name));
  this->Name = name;

  size_t dimLen;
  CALL_NETCDF_GW(nc_inq_dimlen(ncFD, this->DimId, &dimLen));
  this->Coordinates = vtkSmartPointer<vtkDoubleArray>::New();
  this->Coordinates->SetName((this->Name + "_Coordinates").c_str());
  this->Coordinates->SetNumberOfComponents(1);
  this->Coordinates->SetNumberOfTuples(static_cast<vtkIdType>(dimLen));

  this->Bounds = vtkSmartPointer<vtkDoubleArray>::New();
  this->Bounds->SetName((this->Name + "_Bounds").c_str());
  this->Bounds->SetNumberOfComponents(1);
  this->Bounds->SetNumberOfTuples(static_cast<vtkIdType>(dimLen + 1));

  this->SpecialVariables = vtkSmartPointer<vtkStringArray>::New();

  int varId;
  int varNumDims;
  int varDim;
  // By convention if there is a single dimension variable with the same name as
  // its dimension, then the data contains the coordinates for the dimension.
  if ((nc_inq_varid(ncFD, name, &varId) == NC_NOERR) &&
    (nc_inq_varndims(ncFD, varId, &varNumDims) == NC_NOERR) && (varNumDims == 1) &&
    (nc_inq_vardimid(ncFD, varId, &varDim) == NC_NOERR) && (varDim == this->DimId))
  {
    this->SpecialVariables->InsertNextValue(name);

    // Read coordinates
    CALL_NETCDF_GW(nc_get_var_double(ncFD, varId, this->Coordinates->GetPointer(0)));

    // Check to see if the spacing is regular.
    this->Origin = this->Coordinates->GetValue(0);
    this->Spacing =
      (this->Coordinates->GetValue(static_cast<vtkIdType>(dimLen - 1)) - this->Origin) /
      (dimLen - 1);
    this->HasRegularSpacing = true; // Then check to see if it is false.
    double tolerance = 0.01 * this->Spacing;
    for (size_t i = 1; i < dimLen; i++)
    {
      double expectedValue = this->Origin + i * this->Spacing;
      double actualValue = this->Coordinates->GetValue(static_cast<vtkIdType>(i));
      if ((actualValue < expectedValue - tolerance) || (actualValue > expectedValue + tolerance))
      {
        this->HasRegularSpacing = false;
        break;
      }
    }

    // Check units.
    vtkStdString units;
    if (ReadTextAttribute(ncFD, varId, "units", units))
    {
      units = vtksys::SystemTools::LowerCase(units);
      // Time, latitude, and longitude dimensions are those with units that
      // correspond to strings formatted with the Unidata udunits package.  I'm
      // not sure if these checks are complete, but they matches all of the
      // examples I have seen.
      if ((units.find(" since ") != vtkStdString::npos) ||
        (units.find(" after ") != vtkStdString::npos) || (units == "second") ||
        (units == "seconds") || (units == "day") || (units == "days") || (units == "hour") ||
        (units == "hours") || (units == "minute") || (units == "minutes") || (units == "s") ||
        (units == "sec") || (units == "secs") || (units == "shake") || (units == "shakes") ||
        (units == "sidereal_day") || (units == "sidereal_days") || (units == "sidereal_hour") ||
        (units == "sidereal_hours") || (units == "sidereal_minute") ||
        (units == "sidereal_minutes") || (units == "sidereal_second") ||
        (units == "sidereal_seconds") || (units == "sidereal_year") ||
        (units == "sidereal_years") || (units == "tropical_year") || (units == "tropical_years") ||
        (units == "lunar_month") || (units == "lunar_months") || (units == "common_year") ||
        (units == "common_years") || (units == "leap_year") || (units == "leap_years") ||
        (units == "Julian_year") || (units == "Julian_years") || (units == "Gregorian_year") ||
        (units == "Gregorian_years") || (units == "sidereal_month") ||
        (units == "sidereal_months") || (units == "tropical_month") ||
        (units == "tropical_months") || (units == "d") || (units == "min") || (units == "mins") ||
        (units == "hrs") || (units == "h") || (units == "fortnight") || (units == "fortnights") ||
        (units == "week") || (units == "jiffy") || (units == "jiffies") || (units == "year") ||
        (units == "years") || (units == "yr") || (units == "yrs") || (units == "a") ||
        (units == "eon") || (units == "eons") || (units == "month") || (units == "months"))
      {
        this->Units = TIME_UNITS;
      }
      else if (vtksys::RegularExpression("degrees?_?n").find(units))
      {
        this->Units = LATITUDE_UNITS;
      }
      else if (vtksys::RegularExpression("degrees?_?e").find(units))
      {
        this->Units = LONGITUDE_UNITS;
      }
    }

    // Check axis.
    vtkStdString axis;
    if (ReadTextAttribute(ncFD, varId, "axis", axis))
    {
      // The axis attribute is an alternate way of defining the coordinate type.
      // The string can be "X", "Y", "Z", or "T" which mean longitude, latitude,
      // vertical, and time, respectively.
      if (axis == "X")
      {
        this->Units = LONGITUDE_UNITS;
      }
      else if (axis == "Y")
      {
        this->Units = LATITUDE_UNITS;
      }
      else if (axis == "Z")
      {
        this->Units = VERTICAL_UNITS;
      }
      else if (axis == "T")
      {
        this->Units = TIME_UNITS;
      }
    }

    // Check positive.
    vtkStdString positive;
    if (ReadTextAttribute(ncFD, varId, "positive", positive))
    {
      positive = vtksys::SystemTools::LowerCase(positive);
      if (positive.find("down") != vtkStdString::npos)
      {
        // Flip the values of the coordinates.
        for (vtkIdType i = 0; i < this->Coordinates->GetNumberOfTuples(); i++)
        {
          this->Coordinates->SetValue(i, -(this->Coordinates->GetValue(i)));
        }
        this->Spacing = -this->Spacing;
      }
    }

    // Create the bounds array, which is used in place of the coordinates when
    // loading as cell data.  We will look for the bounds attribute on the
    // description variable to see if cell bounds have been written out.  This
    // code assumes that if such an attribute exists, it is the name of another
    // variable that is of dimensions of size dimLen X 2.  There are no checks
    // for this (other than the existence of the attribute), so if this is not
    // the case then the code could fail.
    vtkStdString boundsName;
    if (ReadTextAttribute(ncFD, varId, "bounds", boundsName))
    {
      this->SpecialVariables->InsertNextValue(boundsName);

      int boundsVarId;
      CALL_NETCDF_GW(nc_inq_varid(ncFD, boundsName.c_str(), &boundsVarId));

      // Read in the first bound value for each entry as a point bound.  If the
      // cells are connected, the second bound value should equal the first
      // bound value of the next entry anyway.
      size_t start[2];
      start[0] = start[1] = 0;
      size_t count[2];
      count[0] = dimLen;
      count[1] = 1;
      CALL_NETCDF_GW(
        nc_get_vars_double(ncFD, boundsVarId, start, count, nullptr, this->Bounds->GetPointer(0)));

      // Read in the last value for the bounds array.  It will be the second
      // bound in the last entry.  This will not be replicated unless the
      // dimension is a longitudinal one that wraps all the way around.
      start[0] = dimLen - 1;
      start[1] = 1;
      count[0] = 1;
      count[1] = 1;
      CALL_NETCDF_GW(nc_get_vars_double(ncFD, boundsVarId, start, count, nullptr,
        this->Bounds->GetPointer(static_cast<vtkIdType>(dimLen))));
    }
    else
    {
      // Bounds not given.  Set them based on the coordinates.
      this->Bounds->SetValue(0, this->Coordinates->GetValue(0) - 0.5 * this->Spacing);
      for (vtkIdType i = 1; i < static_cast<vtkIdType>(dimLen); i++)
      {
        double v0 = this->Coordinates->GetValue(i - 1);
        double v1 = this->Coordinates->GetValue(i);
        this->Bounds->SetValue(i, 0.5 * (v0 + v1));
      }
      this->Bounds->SetValue(static_cast<vtkIdType>(dimLen),
        this->Coordinates->GetValue(static_cast<vtkIdType>(dimLen - 1)) + 0.5 * this->Spacing);
    }
  }
  else
  {
    // Fake coordinates
    for (size_t i = 0; i < dimLen; i++)
    {
      this->Coordinates->SetValue(static_cast<vtkIdType>(i), static_cast<double>(i));
      this->Bounds->SetValue(static_cast<vtkIdType>(i), static_cast<double>(i) - 0.5);
    }
    this->Bounds->SetValue(static_cast<vtkIdType>(dimLen), static_cast<double>(dimLen) - 0.5);
    this->HasRegularSpacing = true;
    this->Origin = 0.0;
    this->Spacing = 1.0;
  }

  return 1;
}

//------------------------------------------------------------------------------
class vtkNetCDFCFReader::vtkDimensionInfoVector
{
public:
  std::vector<vtkDimensionInfo> v;
};

//=============================================================================
vtkNetCDFCFReader::vtkDependentDimensionInfo::vtkDependentDimensionInfo(
  int ncFD, int varId, vtkNetCDFCFReader* parent)
{
  if (this->LoadMetaData(ncFD, varId, parent))
  {
    this->Valid = true;
  }
  else
  {
    this->Valid = false;
  }
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::vtkDependentDimensionInfo::LoadMetaData(
  int ncFD, int varId, vtkNetCDFCFReader* parent)
{
  int longitudeCoordVarId, latitudeCoordVarId;
  int longitudeBoundsVarId, latitudeBoundsVarId;
  longitudeCoordVarId = latitudeCoordVarId = -1;
  longitudeBoundsVarId = latitudeBoundsVarId = -1;
  this->GridDimensions = vtkSmartPointer<vtkIntArray>::New();
  this->SpecialVariables = vtkSmartPointer<vtkStringArray>::New();

  // The grid dimensions are the dimensions on the variable.  Since multiple
  // variables can be put on the same grid and this class identifies grids by
  // their variables, I group all of the dimension combinations together for 2D
  // coordinate lookup.  Technically, the CF specification allows you to specify
  // different coordinate variables, but we do not support that because there is
  // no easy way to differentiate grids with the same dimensions.  If different
  // coordinates are needed, then duplicate dimensions should be created.
  // Anyone who disagrees should write their own reader class.
  int numGridDimensions;
  CALL_NETCDF_GW(nc_inq_varndims(ncFD, varId, &numGridDimensions));

  if (numGridDimensions == 0)
  {
    // If a variable has no dimensions, there is no reason to have dependent
    // dimension variables. Just exit here for safety.
    return 0;
  }

  this->GridDimensions->SetNumberOfTuples(numGridDimensions);
  CALL_NETCDF_GW(nc_inq_vardimid(ncFD, varId, this->GridDimensions->GetPointer(0)));

  // Remove initial time dimension, which has no effect on data type.
  if (parent->IsTimeDimension(ncFD, this->GridDimensions->GetValue(0)))
  {
    this->GridDimensions->RemoveTuple(0);
    numGridDimensions--;
    if (numGridDimensions == 0)
    {
      // If a variable has no dimensions (ignoring time), there is no reason
      // to have dependent dimension variables. Just exit here for safety.
      return 0;
    }
  }

  // Most coordinate variables are defined by a variable the same name as the
  // dimension they describe.  Those are handled elsewhere.  This class handles
  // dependent variables that define coordinates that are not the same name as
  // any dimension.  This is only done when the coordinates cannot be expressed
  // as a 1D table lookup from dimension index.  This occurs in only two places
  // in the CF convention.  First, it happens for 2D coordinate variables with
  // 4-sided cells.  This is basically when the grid is a 2D curvilinear grid.
  // Each i,j topological point can be placed anywhere in space.  Second, it
  // happens for multi-dimensional coordinate variables with p-sided cells.
  // These are unstructured collections of polygons.

  vtkStdString coordinates;
  if (!ReadTextAttribute(ncFD, varId, "coordinates", coordinates))
    return 0;

  std::vector<std::string> coordName;
  vtksys::SystemTools::Split(coordinates, coordName, ' ');

  int numAuxCoordDims = -1;

  for (std::vector<std::string>::iterator iter = coordName.begin(); iter != coordName.end(); ++iter)
  {
    int auxCoordVarId;
    if (nc_inq_varid(ncFD, iter->c_str(), &auxCoordVarId) != NC_NOERR)
      continue;

    // Make sure that the coordinate variables have the same dimensions and that
    // those dimensions are the same as the last two dimensions on the grid.
    // Not sure if that is enforced by the specification, but I am going to make
    // that assumption.
    int numDims;
    CALL_NETCDF_GW(nc_inq_varndims(ncFD, auxCoordVarId, &numDims));
    // I am only supporting either 1 or 2 dimensions in the coordinate
    // variables.  See the comment below regarding identifying the
    // CellsUnstructured flag.
    if (numDims > 2)
      continue;

    int auxCoordDims[2];
    CALL_NETCDF_GW(nc_inq_vardimid(ncFD, auxCoordVarId, auxCoordDims));
    int* gridDims = this->GridDimensions->GetPointer(numGridDimensions - numDims);
    bool auxCoordDimsValid = true;
    for (int dimId = 0; dimId < numDims; dimId++)
    {
      if (auxCoordDims[dimId] != gridDims[dimId])
      {
        auxCoordDimsValid = false;
        break;
      }
    }
    if (!auxCoordDimsValid)
      continue;

    // The variable is no use to me unless it is identified as either longitude
    // or latitude.
    vtkStdString units;
    if (!ReadTextAttribute(ncFD, auxCoordVarId, "units", units))
      continue;
    units = vtksys::SystemTools::LowerCase(units);
    if (vtksys::RegularExpression("degrees?_?n").find(units))
    {
      latitudeCoordVarId = auxCoordVarId;
    }
    else if (vtksys::RegularExpression("degrees?_?e").find(units))
    {
      longitudeCoordVarId = auxCoordVarId;
    }
    else
    {
      continue;
    }
    this->SpecialVariables->InsertNextValue(*iter);

    if (numAuxCoordDims < 0)
    {
      // Record the number of dimensions in the coordinate arrays.
      numAuxCoordDims = numDims;
    }
    else if (numAuxCoordDims != numDims)
    {
      // Number of dimensions in different coordinate arrays do not match.
      return 0;
    }
  }

  if ((longitudeCoordVarId == -1) || (latitudeCoordVarId == -1))
  {
    // Did not find all coordinate variables.
    return 0;
  }

  // Technically, p-sided cells can be accessed with any number of dimensions.
  // However, it makes little sense to have more than a 1D array of cell ids.
  // Supporting more than this is a pain because it exceeds the dimensionality
  // supported by vtkDataArray.  It also makes it exceedingly difficult to
  // determine whether the topology is implicitly defined by 2D 4-sided cells or
  // explicitly with p-sided cells.  Thus, I am only implementing p-sided cells
  // with 1D coordinate variables.
  if (numAuxCoordDims == 1)
  {
    this->CellsUnstructured = true;
  }
  else if (numAuxCoordDims == 2)
  {
    this->CellsUnstructured = false;
  }
  else
  {
    // Not supporting this.
    return 0;
  }

  vtkStdString bounds;
  if (ReadTextAttribute(ncFD, longitudeCoordVarId, "bounds", bounds))
  {
    // The bounds is supposed to point to an array with numAuxCoordDims+1
    // dimensions.  The first numAuxCoordDims should be the same as the coord
    // arrays.  The last dimension has the number of vertices in each cell.
    // Maybe I should check this, but I'm not.
    CALL_NETCDF_GW(nc_inq_varid(ncFD, bounds.c_str(), &longitudeBoundsVarId));
    this->SpecialVariables->InsertNextValue(bounds);
  }
  if (ReadTextAttribute(ncFD, latitudeCoordVarId, "bounds", bounds))
  {
    // The bounds is supposed to point to an array with numAuxCoordDims+1
    // dimensions.  The first numAuxCoordDims should be the same as the coord
    // arrays.  The last dimension has the number of vertices in each cell.
    // Maybe I should check this, but I'm not.
    CALL_NETCDF_GW(nc_inq_varid(ncFD, bounds.c_str(), &latitudeBoundsVarId));
    this->SpecialVariables->InsertNextValue(bounds);
  }

  this->HasBounds = ((longitudeBoundsVarId != -1) && (latitudeBoundsVarId != -1));

  // Load in all the longitude and latitude coordinates.  Maybe not the most
  // efficient thing to do for large data, but it is just a 2D surface, so it
  // should be OK for most things.
  this->LongitudeCoordinates = vtkSmartPointer<vtkDoubleArray>::New();
  this->LatitudeCoordinates = vtkSmartPointer<vtkDoubleArray>::New();
  if (this->CellsUnstructured)
  {
    if (this->HasBounds)
    {
      int ok;
      ok = this->LoadUnstructuredBoundsVariable(
        ncFD, longitudeBoundsVarId, this->LongitudeCoordinates);
      if (!ok)
        return 0;
      ok =
        this->LoadUnstructuredBoundsVariable(ncFD, latitudeBoundsVarId, this->LatitudeCoordinates);
      if (!ok)
        return 0;
    }
    else
    {
      // Unstructured cells need bounds to define the topology.
      return 0;
    }
  }
  else
  {
    if (this->HasBounds)
    {
      int ok;
      ok = this->LoadBoundsVariable(ncFD, longitudeBoundsVarId, this->LongitudeCoordinates);
      if (!ok)
        return 0;
      ok = this->LoadBoundsVariable(ncFD, latitudeBoundsVarId, this->LatitudeCoordinates);
      if (!ok)
        return 0;
    }
    else
    {
      int ok;
      ok = this->LoadCoordinateVariable(ncFD, longitudeCoordVarId, this->LongitudeCoordinates);
      if (!ok)
        return 0;
      ok = this->LoadCoordinateVariable(ncFD, latitudeCoordVarId, this->LatitudeCoordinates);
      if (!ok)
        return 0;
    }
  }

  return 1;
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::vtkDependentDimensionInfo::LoadCoordinateVariable(
  int ncFD, int varId, vtkDoubleArray* coords)
{
  int dimIds[2];
  CALL_NETCDF_GW(nc_inq_vardimid(ncFD, varId, dimIds));

  size_t dimSizes[2];
  for (int i = 0; i < 2; i++)
  {
    CALL_NETCDF_GW(nc_inq_dimlen(ncFD, dimIds[i], &dimSizes[i]));
  }

  coords->SetNumberOfComponents(static_cast<int>(dimSizes[1]));
  coords->SetNumberOfTuples(static_cast<vtkIdType>(dimSizes[0]));
  CALL_NETCDF_GW(nc_get_var_double(ncFD, varId, coords->GetPointer(0)));

  return 1;
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::vtkDependentDimensionInfo::LoadBoundsVariable(
  int ncFD, int varId, vtkDoubleArray* coords)
{
  int dimIds[3];
  CALL_NETCDF_GW(nc_inq_vardimid(ncFD, varId, dimIds));

  size_t dimSizes[3];
  for (int i = 0; i < 3; i++)
  {
    CALL_NETCDF_GW(nc_inq_dimlen(ncFD, dimIds[i], &dimSizes[i]));
  }

  if (dimSizes[2] != 4)
  {
    vtkGenericWarningMacro(<< "Expected 2D dependent coordinate bounds to have"
                           << " 4 entries in final dimension.  Instead has " << dimSizes[2]);
    return 0;
  }

  // Bounds are stored as 4-tuples for every cell.  Tuple entries 0 and 1
  // connect to the cell in the -i topological direction.  Tuple entries 0 and 3
  // connect to the cell in the -j topological direction.
  std::vector<double> boundsData(dimSizes[0] * dimSizes[1] * 4);
  if (!boundsData.empty())
  {
    CALL_NETCDF_GW(nc_get_var_double(ncFD, varId, &boundsData.at(0)));
  }

  // The coords array are the coords at the points.  There is one more point
  // than cell in each topological direction.
  int numComponents = static_cast<int>(dimSizes[1]);
  vtkIdType numTuples = static_cast<vtkIdType>(dimSizes[0]);
  coords->SetNumberOfComponents(numComponents + 1);
  coords->SetNumberOfTuples(numTuples + 1);

  // Copy from the bounds data to the coordinates data.  Most values will
  // be copied from the bound's 0'th tuple entry.  Values at the extremes
  // will be copied from other entries.
  for (vtkIdType j = 0; j < numTuples; j++)
  {
    for (int i = 0; i < numComponents; i++)
    {
      coords->SetComponent(j, i, boundsData[(j * numComponents + i) * 4 + 0]);
    }
    coords->SetComponent(j, numComponents, boundsData[((j + 1) * numComponents - 1) * 4 + 1]);
  }
  for (int i = 0; i < numComponents; i++)
  {
    coords->SetComponent(numTuples, i, boundsData[((numTuples - 1) * numComponents) * 4 + 2]);
  }
  coords->SetComponent(
    numTuples, numComponents, boundsData[(numTuples * numComponents - 1) * 4 + 3]);

  return 1;
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::vtkDependentDimensionInfo::LoadUnstructuredBoundsVariable(
  int ncFD, int varId, vtkDoubleArray* coords)
{
  int dimIds[2];
  CALL_NETCDF_GW(nc_inq_vardimid(ncFD, varId, dimIds));

  size_t dimSizes[2];
  for (int i = 0; i < 2; i++)
  {
    CALL_NETCDF_GW(nc_inq_dimlen(ncFD, dimIds[i], &dimSizes[i]));
  }

  int numVertPerCell = static_cast<int>(dimSizes[1]);
  vtkIdType numCells = static_cast<vtkIdType>(dimSizes[0]);

  coords->SetNumberOfComponents(numVertPerCell);
  coords->SetNumberOfTuples(numCells);
  CALL_NETCDF_GW(nc_get_var_double(ncFD, varId, coords->GetPointer(0)));

  return 1;
}

//------------------------------------------------------------------------------
class vtkNetCDFCFReader::vtkDependentDimensionInfoVector
{
public:
  std::vector<vtkNetCDFCFReader::vtkDependentDimensionInfo> v;
};

//=============================================================================
vtkStandardNewMacro(vtkNetCDFCFReader);

//------------------------------------------------------------------------------
vtkNetCDFCFReader::vtkNetCDFCFReader()
{
  this->SphericalCoordinates = 1;
  this->VerticalScale = 1.0;
  this->VerticalBias = 0.0;
  this->OutputType = -1;

  this->DimensionInfo = new vtkDimensionInfoVector;
  this->DependentDimensionInfo = new vtkDependentDimensionInfoVector;
}

vtkNetCDFCFReader::~vtkNetCDFCFReader()
{
  delete this->DimensionInfo;
  delete this->DependentDimensionInfo;
}

void vtkNetCDFCFReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "SphericalCoordinates: " << this->SphericalCoordinates << endl;
  os << indent << "VerticalScale: " << this->VerticalScale << endl;
  os << indent << "VerticalBias: " << this->VerticalBias << endl;
  os << indent << "OutputType: " << this->OutputType << endl;
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::CanReadFile(const char* filename)
{
  // We really just read basic arrays from netCDF files.  If the netCDF library
  // says we can read it, then we can read it.
  int ncFD;
  int errorcode = nc_open(filename, NC_NOWRITE, &ncFD);
  if (errorcode == NC_NOERR)
  {
    nc_close(ncFD);
    return 1;
  }
  else
  {
    return 0;
  }
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::SetOutputType(int type)
{
  vtkDebugMacro(<< this->GetClassName() << " (" << this
                << "):"
                   " setting OutputType to "
                << type);
  if (this->OutputType != type)
  {
    bool typeValid = ((type == -1) || (type == VTK_IMAGE_DATA) || (type == VTK_RECTILINEAR_GRID) ||
      (type == VTK_STRUCTURED_GRID) || (type == VTK_UNSTRUCTURED_GRID));
    if (typeValid)
    {
      this->OutputType = type;
      this->Modified();
    }
    else
    {
      vtkErrorMacro(<< "Invalid OutputType: " << type);
    }
  }
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::RequestDataObject(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
{
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkDataObject* output = vtkDataObject::GetData(outInfo);

  // This is really too early to know the appropriate data type.  We need to
  // have meta data and let the user select arrays.  We have to do part
  // of the RequestInformation to get the appropriate meta data.
  if (!this->UpdateMetaData())
    return 0;

  // Check that the dataType is correct or automatically set it if it is set to
  // -1.
  int dataType = this->OutputType;

  int ncFD;
  CALL_NETCDF(nc_open(this->FileName, NC_NOWRITE, &ncFD));

  int numArrays = this->VariableArraySelection->GetNumberOfArrays();
  for (int arrayIndex = 0; arrayIndex < numArrays; arrayIndex++)
  {
    if (!this->VariableArraySelection->GetArraySetting(arrayIndex))
      continue;

    const char* name = this->VariableArraySelection->GetArrayName(arrayIndex);
    int varId;
    CALL_NETCDF(nc_inq_varid(ncFD, name, &varId));

    int currentNumDims;
    CALL_NETCDF(nc_inq_varndims(ncFD, varId, &currentNumDims));
    if (currentNumDims < 1)
      continue;
    VTK_CREATE(vtkIntArray, currentDimensions);
    currentDimensions->SetNumberOfComponents(1);
    currentDimensions->SetNumberOfTuples(currentNumDims);
    CALL_NETCDF(nc_inq_vardimid(ncFD, varId, currentDimensions->GetPointer(0)));

    // Remove initial time dimension, which has no effect on data type.
    if (this->IsTimeDimension(ncFD, currentDimensions->GetValue(0)))
    {
      currentDimensions->RemoveTuple(0);
      currentNumDims--;
      if (currentNumDims < 1)
        continue;
    }

    CoordinateTypesEnum coordType = this->CoordinateType(currentDimensions);

    int preferredDataType;
    switch (coordType)
    {
      case COORDS_UNIFORM_RECTILINEAR:
        preferredDataType = VTK_IMAGE_DATA;
        break;
      case COORDS_NONUNIFORM_RECTILINEAR:
        preferredDataType = VTK_RECTILINEAR_GRID;
        break;
      case COORDS_REGULAR_SPHERICAL:
      case COORDS_2D_EUCLIDEAN:
      case COORDS_2D_SPHERICAL:
      case COORDS_EUCLIDEAN_4SIDED_CELLS:
      case COORDS_SPHERICAL_4SIDED_CELLS:
        preferredDataType = VTK_STRUCTURED_GRID;
        break;
      case COORDS_EUCLIDEAN_PSIDED_CELLS:
      case COORDS_SPHERICAL_PSIDED_CELLS:
        preferredDataType = VTK_UNSTRUCTURED_GRID;
        break;
      default:
        vtkErrorMacro(<< "Internal error: unknown coordinate type.");
        return 0;
    }

    // Check the data type.
    if (dataType == -1)
    {
      dataType = preferredDataType;
    }
    else
    {
      switch (dataType)
      {
        case VTK_IMAGE_DATA:
          if (preferredDataType != VTK_IMAGE_DATA)
          {
            vtkWarningMacro("You have set the OutputType to a data type that"
                            " cannot fully represent the topology of the data."
                            " Some of the topology will be ignored.");
          }
          break;
        case VTK_RECTILINEAR_GRID:
          if ((preferredDataType != VTK_IMAGE_DATA) && (preferredDataType != VTK_RECTILINEAR_GRID))
          {
            vtkWarningMacro("You have set the OutputType to a data type that"
                            " cannot fully represent the topology of the data."
                            " Some of the topology will be ignored.");
          }
          break;
        case VTK_STRUCTURED_GRID:
          if ((preferredDataType != VTK_IMAGE_DATA) &&
            (preferredDataType != VTK_RECTILINEAR_GRID) &&
            (preferredDataType != VTK_STRUCTURED_GRID))
          {
            vtkWarningMacro("You have set the OutputType to a data type that"
                            " cannot fully represent the topology of the data."
                            " Some of the topology will be ignored.");
          }
          break;
        case VTK_UNSTRUCTURED_GRID:
          // Unstructured grid supports all topologies.
          break;
        default:
          vtkErrorMacro(<< "Sanity check failed: bad internal type.");
          return 0;
      }
    }

    // That's right, break.  We really only want to look at the dimensions of
    // the first valid loaded variable.  If we got here, we found a variable.
    // The loop is really only for the continues at the top.
    break;
  }

  CALL_NETCDF(nc_close(ncFD));

  if (dataType == -1)
  {
    // Either no variables are selected or only variables with no dimensions.
    // In either case, an image data should be fine.
    dataType = VTK_IMAGE_DATA;
  }

  if (dataType == VTK_IMAGE_DATA)
  {
    if (!output || !output->IsA("vtkImageData"))
    {
      output = vtkImageData::New();
      outInfo->Set(vtkDataObject::DATA_OBJECT(), output);
      output->Delete(); // Not really deleted.
    }
  }
  else if (dataType == VTK_RECTILINEAR_GRID)
  {
    if (!output || !output->IsA("vtkRectilinearGrid"))
    {
      output = vtkRectilinearGrid::New();
      outInfo->Set(vtkDataObject::DATA_OBJECT(), output);
      output->Delete(); // Not really deleted.
    }
  }
  else if (dataType == VTK_STRUCTURED_GRID)
  {
    if (!output || !output->IsA("vtkStructuredGrid"))
    {
      output = vtkStructuredGrid::New();
      outInfo->Set(vtkDataObject::DATA_OBJECT(), output);
      output->Delete(); // Not really deleted.
    }
  }
  else if (dataType == VTK_UNSTRUCTURED_GRID)
  {
    if (!output || !output->IsA("vtkUnstructuredGrid"))
    {
      output = vtkUnstructuredGrid::New();
      outInfo->Set(vtkDataObject::DATA_OBJECT(), output);
      output->Delete(); // Not really deleted.
    }
  }
  else
  {
    vtkErrorMacro(<< "Sanity check failed: bad internal type.");
    return 0;
  }

  return 1;
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::RequestInformation(
  vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // Let the superclass do the heavy lifting.
  if (!this->Superclass::RequestInformation(request, inputVector, outputVector))
  {
    return 0;
  }

  // Superclass understands structured data, but we have to handle unstructured
  // "extents" (pieces).
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkDataObject* output = vtkDataObject::GetData(outInfo);
  if (output)
  {
    if (output->GetExtentType() != VTK_3D_EXTENT)
    {
      outInfo->Set(CAN_HANDLE_PIECE_REQUEST(), 1);
    }
    else
    {
      outInfo->Set(CAN_PRODUCE_SUB_EXTENT(), 1);
    }
  }
  else
  {
    return 0;
  }
  return 1;
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::RequestData(
  vtkInformation* request, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // If the output does not directly support 3D extents, then we have to make
  // some from the piece information so the superclass knows what portion of
  // arrays to load.
  vtkDataObject* output = vtkDataObject::GetData(outputVector);
  if (output)
  {
    if (output->GetExtentType() == VTK_3D_EXTENT)
    {
      // Do nothing.  3D extents already set.
    }
    else if (output->GetExtentType() == VTK_PIECES_EXTENT)
    {
      int pieceNumber, numberOfPieces, ghostLevels;
      vtkInformation* outInfo = outputVector->GetInformationObject(0);
      pieceNumber = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
      numberOfPieces = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
      ghostLevels = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());

      int extent[6];
      this->ExtentForDimensionsAndPiece(pieceNumber, numberOfPieces, ghostLevels, extent);

      // Store the update extent in the output's information object to make it
      // easy to find whenever loading data for this object.
      output->GetInformation()->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent, 6);
    }
    else
    {
      vtkWarningMacro(<< "Invalid extent type encountered.  Data arrays may"
                      << " be loaded incorrectly.");
    }
  }
  else // output == nullptr
  {
    vtkErrorMacro(<< "No output object.");
    return 0;
  }

  // Let the superclass do the heavy lifting.
  if (!this->Superclass::RequestData(request, inputVector, outputVector))
  {
    return 0;
  }

  // Add spacing information defined by the COARDS conventions.

  vtkImageData* imageOutput = vtkImageData::GetData(outputVector);
  if (imageOutput)
  {
    this->AddRectilinearCoordinates(imageOutput);
  }

  vtkRectilinearGrid* rectilinearOutput = vtkRectilinearGrid::GetData(outputVector);
  if (rectilinearOutput)
  {
    switch (this->CoordinateType(this->LoadingDimensions))
    {
      case COORDS_EUCLIDEAN_PSIDED_CELLS:
      case COORDS_SPHERICAL_PSIDED_CELLS:
        // There is no sensible way to store p-sided cells in a structured grid.
        // Just fake some coordinates (related to ParaView bug #11543).
        this->FakeRectilinearCoordinates(rectilinearOutput);
        break;
      case COORDS_UNIFORM_RECTILINEAR:
      case COORDS_NONUNIFORM_RECTILINEAR:
      case COORDS_REGULAR_SPHERICAL:
      case COORDS_2D_EUCLIDEAN:
      case COORDS_2D_SPHERICAL:
      case COORDS_EUCLIDEAN_4SIDED_CELLS:
      case COORDS_SPHERICAL_4SIDED_CELLS:
      default:
        this->AddRectilinearCoordinates(rectilinearOutput);
    }
  }

  vtkStructuredGrid* structuredOutput = vtkStructuredGrid::GetData(outputVector);
  if (structuredOutput)
  {
    switch (this->CoordinateType(this->LoadingDimensions))
    {
      case COORDS_UNIFORM_RECTILINEAR:
      case COORDS_NONUNIFORM_RECTILINEAR:
        this->Add1DRectilinearCoordinates(structuredOutput);
        break;
      case COORDS_REGULAR_SPHERICAL:
        this->Add1DSphericalCoordinates(structuredOutput);
        break;
      case COORDS_2D_EUCLIDEAN:
      case COORDS_EUCLIDEAN_4SIDED_CELLS:
        this->Add2DRectilinearCoordinates(structuredOutput);
        break;
      case COORDS_2D_SPHERICAL:
      case COORDS_SPHERICAL_4SIDED_CELLS:
        this->Add2DSphericalCoordinates(structuredOutput);
        break;
      case COORDS_EUCLIDEAN_PSIDED_CELLS:
      case COORDS_SPHERICAL_PSIDED_CELLS:
        // There is no sensible way to store p-sided cells in a structured grid.
        // Just fake some coordinates (ParaView bug #11543).
        this->FakeStructuredCoordinates(structuredOutput);
        break;
      default:
        vtkErrorMacro("Internal error: unknown coordinate type.");
        return 0;
    }
  }

  vtkUnstructuredGrid* unstructuredOutput = vtkUnstructuredGrid::GetData(outputVector);
  if (unstructuredOutput)
  {
    int extent[6];
    this->GetUpdateExtentForOutput(unstructuredOutput, extent);

    switch (this->CoordinateType(this->LoadingDimensions))
    {
      case COORDS_UNIFORM_RECTILINEAR:
      case COORDS_NONUNIFORM_RECTILINEAR:
        this->Add1DRectilinearCoordinates(unstructuredOutput, extent);
        break;
      case COORDS_REGULAR_SPHERICAL:
        this->Add1DSphericalCoordinates(unstructuredOutput, extent);
        break;
      case COORDS_2D_EUCLIDEAN:
      case COORDS_EUCLIDEAN_4SIDED_CELLS:
        this->Add2DRectilinearCoordinates(unstructuredOutput, extent);
        break;
      case COORDS_2D_SPHERICAL:
      case COORDS_SPHERICAL_4SIDED_CELLS:
        this->Add2DSphericalCoordinates(unstructuredOutput, extent);
        break;
      case COORDS_EUCLIDEAN_PSIDED_CELLS:
        this->AddUnstructuredRectilinearCoordinates(unstructuredOutput, extent);
        break;
      case COORDS_SPHERICAL_PSIDED_CELLS:
        this->AddUnstructuredSphericalCoordinates(unstructuredOutput, extent);
        break;
      default:
        vtkErrorMacro("Internal error: unknown coordinate type.");
        return 0;
    }
  }

  return 1;
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::ExtentForDimensionsAndPiece(
  int pieceNumber, int numberOfPieces, int ghostLevels, int extent[6])
{
  VTK_CREATE(vtkExtentTranslator, extentTranslator);

  extentTranslator->SetWholeExtent(this->WholeExtent);
  extentTranslator->SetPiece(pieceNumber);
  extentTranslator->SetNumberOfPieces(numberOfPieces);
  extentTranslator->SetGhostLevel(ghostLevels);

  extentTranslator->PieceToExtent();

  extentTranslator->GetExtent(extent);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::GetUpdateExtentForOutput(vtkDataSet* output, int extent[6])
{
  vtkInformation* info = output->GetInformation();
  if (info->Has(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT()))
  {
    info->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent);
  }
  else
  {
    return this->Superclass::GetUpdateExtentForOutput(output, extent);
  }
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::AddRectilinearCoordinates(vtkImageData* imageOutput)
{
  double origin[3];
  origin[0] = origin[1] = origin[2] = 0.0;
  double spacing[3];
  spacing[0] = spacing[1] = spacing[2] = 1.0;

  int numDim = this->LoadingDimensions->GetNumberOfTuples();
  if (numDim >= 3)
    numDim = 3;

  for (int i = 0; i < numDim; i++)
  {
    // Remember that netCDF dimension ordering is backward from VTK.
    int dim = this->LoadingDimensions->GetValue(numDim - i - 1);
    vtkDimensionInfo* dimInfo = this->GetDimensionInfo(dim);
    origin[i] = dimInfo->GetOrigin();
    spacing[i] = dimInfo->GetSpacing();
  }

  imageOutput->SetOrigin(origin);
  imageOutput->SetSpacing(spacing);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::AddRectilinearCoordinates(vtkRectilinearGrid* rectilinearOutput)
{
  int extent[6];
  rectilinearOutput->GetExtent(extent);

  int numDim = this->LoadingDimensions->GetNumberOfTuples();
  for (int i = 0; i < 3; i++)
  {
    vtkSmartPointer<vtkDoubleArray> coords;
    if (i < numDim)
    {
      // Remember that netCDF dimension ordering is backward from VTK.
      int dim = this->LoadingDimensions->GetValue(numDim - i - 1);
      coords = this->GetDimensionInfo(dim)->GetCoordinates();
      int extLow = extent[2 * i];
      int extHi = extent[2 * i + 1];
      if ((extLow != 0) || (extHi != coords->GetNumberOfTuples() - 1))
      {
        // Getting a subset of this dimension.
        VTK_CREATE(vtkDoubleArray, newcoords);
        newcoords->SetNumberOfComponents(1);
        newcoords->SetNumberOfTuples(extHi - extLow + 1);
        memcpy(newcoords->GetPointer(0), coords->GetPointer(extLow),
          (extHi - extLow + 1) * sizeof(double));
        coords = newcoords;
      }
    }
    else
    {
      coords = vtkSmartPointer<vtkDoubleArray>::New();
      coords->SetNumberOfTuples(1);
      coords->SetComponent(0, 0, 0.0);
    }
    switch (i)
    {
      case 0:
        rectilinearOutput->SetXCoordinates(coords);
        break;
      case 1:
        rectilinearOutput->SetYCoordinates(coords);
        break;
      case 2:
        rectilinearOutput->SetZCoordinates(coords);
        break;
      default:
        vtkErrorMacro("Sanity check failed!");
        break;
    }
  }
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::FakeRectilinearCoordinates(vtkRectilinearGrid* rectilinearOutput)
{
  int dimensionSizes[3];
  rectilinearOutput->GetDimensions(dimensionSizes);

  for (int dim = 0; dim < 3; dim++)
  {
    vtkNew<vtkDoubleArray> coordinate;
    coordinate->SetNumberOfComponents(1);
    coordinate->SetNumberOfTuples(dimensionSizes[dim]);

    for (int index = 0; index < dimensionSizes[dim]; index++)
    {
      coordinate->SetComponent(index, 0, static_cast<double>(index));
    }
    switch (dim)
    {
      case 0:
        rectilinearOutput->SetXCoordinates(coordinate);
        break;
      case 1:
        rectilinearOutput->SetYCoordinates(coordinate);
        break;
      case 2:
        rectilinearOutput->SetZCoordinates(coordinate);
        break;
      default:
        vtkErrorMacro("Sanity check failed!");
        break;
    }
  }
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add1DRectilinearCoordinates(vtkPoints* points, const int extent[6])
{
  points->SetDataTypeToDouble();
  points->SetNumberOfPoints(
    (extent[1] - extent[0] + 1) * (extent[3] - extent[2] + 1) * (extent[5] - extent[4] + 1));
  vtkDataArray* pointData = points->GetData();

  int numDimNetCDF = this->LoadingDimensions->GetNumberOfTuples();
  for (int dimVTK = 0; dimVTK < 3; dimVTK++)
  {
    vtkSmartPointer<vtkDoubleArray> coords;
    if (dimVTK < numDimNetCDF)
    {
      // Remember that netCDF dimension ordering is backward from VTK.
      int dimNetCDF = this->LoadingDimensions->GetValue(numDimNetCDF - dimVTK - 1);
      coords = this->GetDimensionInfo(dimNetCDF)->GetCoordinates();

      int ijk[3];
      vtkIdType pointIdx = 0;
      for (ijk[2] = extent[2 * 2]; ijk[2] <= extent[2 * 2 + 1]; ijk[2]++)
      {
        for (ijk[1] = extent[1 * 2]; ijk[1] <= extent[1 * 2 + 1]; ijk[1]++)
        {
          for (ijk[0] = extent[0 * 2]; ijk[0] <= extent[0 * 2 + 1]; ijk[0]++)
          {
            pointData->SetComponent(pointIdx, dimVTK, coords->GetValue(ijk[dimVTK]));
            pointIdx++;
          }
        }
      }
    }
    else
    {
      int ijk[3];
      vtkIdType pointIdx = 0;
      for (ijk[2] = extent[2 * 2]; ijk[2] <= extent[2 * 2 + 1]; ijk[2]++)
      {
        for (ijk[1] = extent[1 * 2]; ijk[1] <= extent[1 * 2 + 1]; ijk[1]++)
        {
          for (ijk[0] = extent[0 * 2]; ijk[0] <= extent[0 * 2 + 1]; ijk[0]++)
          {
            pointData->SetComponent(pointIdx, dimVTK, 0.0);
            pointIdx++;
          }
        }
      }
    }
  }
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add2DRectilinearCoordinates(vtkPoints* points, const int extent[6])
{
  points->SetDataTypeToDouble();
  points->Allocate(
    (extent[1] - extent[0] + 1) * (extent[3] - extent[2] + 1) * (extent[5] - extent[4] + 1));

  vtkDependentDimensionInfo* info = this->FindDependentDimensionInfo(this->LoadingDimensions);

  vtkDoubleArray* longitudeCoordinates = info->GetLongitudeCoordinates();
  vtkDoubleArray* latitudeCoordinates = info->GetLatitudeCoordinates();

  vtkDoubleArray* verticalCoordinates = nullptr;
  if (this->LoadingDimensions->GetNumberOfTuples() == 3)
  {
    int vertDim = this->LoadingDimensions->GetValue(0);
    if (info->GetHasBounds())
    {
      verticalCoordinates = this->GetDimensionInfo(vertDim)->GetBounds();
    }
    else
    {
      verticalCoordinates = this->GetDimensionInfo(vertDim)->GetCoordinates();
    }
  }

  for (int k = extent[4]; k <= extent[5]; k++)
  {
    double h;
    if (verticalCoordinates)
    {
      h = verticalCoordinates->GetValue(k);
    }
    else
    {
      h = 0.0;
    }
    for (int j = extent[2]; j <= extent[3]; j++)
    {
      for (int i = extent[0]; i <= extent[1]; i++)
      {
        double lon = longitudeCoordinates->GetComponent(j, i);
        double lat = latitudeCoordinates->GetComponent(j, i);
        points->InsertNextPoint(lon, lat, h);
      }
    }
  }
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add1DRectilinearCoordinates(vtkStructuredGrid* structuredOutput)
{
  int extent[6];
  structuredOutput->GetExtent(extent);

  VTK_CREATE(vtkPoints, points);
  this->Add1DRectilinearCoordinates(points, extent);
  structuredOutput->SetPoints(points);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add2DRectilinearCoordinates(vtkStructuredGrid* structuredOutput)
{
  int extent[6];
  structuredOutput->GetExtent(extent);

  VTK_CREATE(vtkPoints, points);
  this->Add2DRectilinearCoordinates(points, extent);
  structuredOutput->SetPoints(points);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::FakeStructuredCoordinates(vtkStructuredGrid* structuredOutput)
{
  int extent[6];
  structuredOutput->GetExtent(extent);

  vtkNew<vtkPoints> points;
  points->SetDataTypeToDouble();
  points->Allocate(
    (extent[1] - extent[0] + 1) * (extent[3] - extent[2] + 1) * (extent[5] - extent[4] + 1));

  for (int kIndex = extent[4]; kIndex <= extent[5]; kIndex++)
  {
    for (int jIndex = extent[2]; jIndex <= extent[3]; jIndex++)
    {
      for (int iIndex = extent[0]; iIndex <= extent[1]; iIndex++)
      {
        points->InsertNextPoint(iIndex, jIndex, kIndex);
      }
    }
  }

  structuredOutput->SetPoints(points);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add1DRectilinearCoordinates(
  vtkUnstructuredGrid* unstructuredOutput, const int extent[6])
{
  VTK_CREATE(vtkPoints, points);
  this->Add1DRectilinearCoordinates(points, extent);
  unstructuredOutput->SetPoints(points);

  this->AddStructuredCells(unstructuredOutput, extent);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add2DRectilinearCoordinates(
  vtkUnstructuredGrid* unstructuredOutput, const int extent[6])
{
  VTK_CREATE(vtkPoints, points);
  this->Add2DRectilinearCoordinates(points, extent);
  unstructuredOutput->SetPoints(points);

  this->AddStructuredCells(unstructuredOutput, extent);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add1DSphericalCoordinates(vtkPoints* points, const int extent[6])
{
  points->SetDataTypeToDouble();
  points->Allocate(
    (extent[1] - extent[0] + 1) * (extent[3] - extent[2] + 1) * (extent[5] - extent[4] + 1));

  vtkDoubleArray* coordArrays[3];
  for (vtkIdType i = 0; i < this->LoadingDimensions->GetNumberOfTuples(); i++)
  {
    int dim = this->LoadingDimensions->GetValue(i);
    coordArrays[i] = this->GetDimensionInfo(dim)->GetBounds();
  }

  int longitudeDim, latitudeDim, verticalDim;
  this->IdentifySphericalCoordinates(
    this->LoadingDimensions, longitudeDim, latitudeDim, verticalDim);

  if ((longitudeDim < 0) || (latitudeDim < 0))
  {
    vtkErrorMacro(<< "Internal error: treating non spherical coordinates as if"
                  << " they were spherical.");
    return;
  }

  // Check the height scale and bias.
  double vertScale = this->VerticalScale;
  double vertBias = this->VerticalBias;
  if (verticalDim >= 0)
  {
    double* verticalRange = coordArrays[verticalDim]->GetRange();
    if ((verticalRange[0] * vertScale + vertBias < 0) ||
      (verticalRange[1] * vertScale + vertBias < 0))
    {
      vertBias = -std::min(verticalRange[0], verticalRange[1]) * vertScale;
    }
  }
  else
  {
    if (vertScale + vertBias <= 0)
    {
      vertScale = 1.0;
      vertBias = 0.0;
    }
  }

  int ijk[3];
  for (ijk[0] = extent[4]; ijk[0] <= extent[5]; ijk[0]++)
  {
    for (ijk[1] = extent[2]; ijk[1] <= extent[3]; ijk[1]++)
    {
      for (ijk[2] = extent[0]; ijk[2] <= extent[1]; ijk[2]++)
      {
        double lon, lat, h;
        if (verticalDim >= 0)
        {
          lon = coordArrays[longitudeDim]->GetValue(ijk[longitudeDim]);
          lat = coordArrays[latitudeDim]->GetValue(ijk[latitudeDim]);
          h = coordArrays[verticalDim]->GetValue(ijk[verticalDim]);
        }
        else
        {
          lon = coordArrays[longitudeDim]->GetValue(ijk[longitudeDim + 1]);
          lat = coordArrays[latitudeDim]->GetValue(ijk[latitudeDim + 1]);
          h = 1.0;
        }
        lon = vtkMath::RadiansFromDegrees(lon);
        lat = vtkMath::RadiansFromDegrees(lat);
        h = h * vertScale + vertBias;

        double cartesianCoord[3];
        cartesianCoord[0] = h * cos(lon) * cos(lat);
        cartesianCoord[1] = h * sin(lon) * cos(lat);
        cartesianCoord[2] = h * sin(lat);
        points->InsertNextPoint(cartesianCoord);
      }
    }
  }
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add2DSphericalCoordinates(vtkPoints* points, const int extent[6])
{
  points->SetDataTypeToDouble();
  points->Allocate(
    (extent[1] - extent[0] + 1) * (extent[3] - extent[2] + 1) * (extent[5] - extent[4] + 1));

  vtkDependentDimensionInfo* info = this->FindDependentDimensionInfo(this->LoadingDimensions);

  vtkDoubleArray* longitudeCoordinates = info->GetLongitudeCoordinates();
  vtkDoubleArray* latitudeCoordinates = info->GetLatitudeCoordinates();

  vtkDoubleArray* verticalCoordinates = nullptr;
  if (this->LoadingDimensions->GetNumberOfTuples() == 3)
  {
    int vertDim = this->LoadingDimensions->GetValue(0);
    if (info->GetHasBounds())
    {
      verticalCoordinates = this->GetDimensionInfo(vertDim)->GetBounds();
    }
    else
    {
      verticalCoordinates = this->GetDimensionInfo(vertDim)->GetCoordinates();
    }
  }

  // Check the height scale and bias.
  double vertScale = this->VerticalScale;
  double vertBias = this->VerticalBias;
  if (verticalCoordinates)
  {
    double* verticalRange = verticalCoordinates->GetRange();
    if ((verticalRange[0] * vertScale + vertBias < 0) ||
      (verticalRange[1] * vertScale + vertBias < 0))
    {
      vertBias = -std::min(verticalRange[0], verticalRange[1]) * vertScale;
    }
  }
  else
  {
    if (vertScale + vertBias <= 0)
    {
      vertScale = 1.0;
      vertBias = 0.0;
    }
  }

  for (int k = extent[4]; k <= extent[5]; k++)
  {
    double h;
    if (verticalCoordinates)
    {
      h = verticalCoordinates->GetValue(k) * vertScale + vertBias;
    }
    else
    {
      h = vertScale + vertBias;
    }
    for (int j = extent[2]; j <= extent[3]; j++)
    {
      for (int i = extent[0]; i <= extent[1]; i++)
      {
        double lon = longitudeCoordinates->GetComponent(j, i);
        double lat = latitudeCoordinates->GetComponent(j, i);
        lon = vtkMath::RadiansFromDegrees(lon);
        lat = vtkMath::RadiansFromDegrees(lat);

        double cartesianCoord[3];
        cartesianCoord[0] = h * cos(lon) * cos(lat);
        cartesianCoord[1] = h * sin(lon) * cos(lat);
        cartesianCoord[2] = h * sin(lat);
        points->InsertNextPoint(cartesianCoord);
      }
    }
  }
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add1DSphericalCoordinates(vtkStructuredGrid* structuredOutput)
{
  int extent[6];
  structuredOutput->GetExtent(extent);

  VTK_CREATE(vtkPoints, points);
  this->Add1DSphericalCoordinates(points, extent);
  structuredOutput->SetPoints(points);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add2DSphericalCoordinates(vtkStructuredGrid* structuredOutput)
{
  int extent[6];
  structuredOutput->GetExtent(extent);

  VTK_CREATE(vtkPoints, points);
  this->Add2DSphericalCoordinates(points, extent);
  structuredOutput->SetPoints(points);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add1DSphericalCoordinates(
  vtkUnstructuredGrid* unstructuredOutput, const int extent[6])
{
  VTK_CREATE(vtkPoints, points);
  this->Add1DSphericalCoordinates(points, extent);
  unstructuredOutput->SetPoints(points);

  this->AddStructuredCells(unstructuredOutput, extent);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::Add2DSphericalCoordinates(
  vtkUnstructuredGrid* unstructuredOutput, const int extent[6])
{
  VTK_CREATE(vtkPoints, points);
  this->Add2DSphericalCoordinates(points, extent);
  unstructuredOutput->SetPoints(points);

  this->AddStructuredCells(unstructuredOutput, extent);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::AddStructuredCells(
  vtkUnstructuredGrid* unstructuredOutput, const int extent[6])
{
  vtkIdType numPoints[3];
  numPoints[0] = extent[1] - extent[0] + 1;
  numPoints[1] = extent[3] - extent[2] + 1;
  numPoints[2] = extent[5] - extent[4] + 1;

  vtkIdType numCells[3];
  numCells[0] = numPoints[0] - 1;
  numCells[1] = numPoints[1] - 1;
  numCells[2] = numPoints[2] - 1;

  vtkIdType nextPointRow = numPoints[0];
  vtkIdType nextPointSlab = nextPointRow * numPoints[1];

  bool extentIs2D = (numCells[2] < 1);

  if (extentIs2D)
  {
    vtkIdType totalNumCells = numCells[0] * numCells[1];
    unstructuredOutput->Allocate(totalNumCells);
    vtkCellArray* cells = unstructuredOutput->GetCells();
    cells->AllocateEstimate(totalNumCells, 4);

    for (int j = 0; j < numCells[1]; j++)
    {
      vtkIdType rowStart = j * nextPointRow;
      for (int i = 0; i < numCells[0]; i++)
      {
        vtkIdType lowCellPoint = rowStart + i;

        vtkIdType pointIds[4];
        pointIds[0] = lowCellPoint;
        pointIds[1] = lowCellPoint + 1;
        pointIds[2] = lowCellPoint + nextPointRow + 1;
        pointIds[3] = lowCellPoint + nextPointRow;

        unstructuredOutput->InsertNextCell(VTK_QUAD, 4, pointIds);
      }
    }
  }
  else // !extentIs2D
  {
    vtkIdType totalNumCells = numCells[0] * numCells[1] * numCells[2];
    unstructuredOutput->Allocate(totalNumCells);
    vtkCellArray* cells = unstructuredOutput->GetCells();
    cells->AllocateEstimate(totalNumCells, 8);

    for (int k = 0; k < numCells[2]; k++)
    {
      vtkIdType slabStart = k * nextPointSlab;
      for (int j = 0; j < numCells[1]; j++)
      {
        vtkIdType rowStart = slabStart + j * nextPointRow;
        for (int i = 0; i < numCells[0]; i++)
        {
          vtkIdType lowCellPoint = rowStart + i;

          // This code is assuming that all axis are scaling up.  If that is
          // not the case, this will probably make inverted hexahedra.
          vtkIdType pointIds[8];
          pointIds[0] = lowCellPoint;
          pointIds[1] = lowCellPoint + 1;
          pointIds[2] = lowCellPoint + nextPointRow + 1;
          pointIds[3] = lowCellPoint + nextPointRow;
          pointIds[4] = lowCellPoint + nextPointSlab;
          pointIds[5] = lowCellPoint + nextPointSlab + 1;
          pointIds[6] = lowCellPoint + nextPointSlab + nextPointRow + 1;
          pointIds[7] = lowCellPoint + nextPointSlab + nextPointRow;

          unstructuredOutput->InsertNextCell(VTK_HEXAHEDRON, 8, pointIds);
        }
      }
    }
  }
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::AddUnstructuredRectilinearCoordinates(
  vtkUnstructuredGrid* unstructuredOutput, const int extent[6])
{
  vtkDependentDimensionInfo* info = this->FindDependentDimensionInfo(this->LoadingDimensions);

  vtkDoubleArray* longitudeCoordinates = info->GetLongitudeCoordinates();
  vtkDoubleArray* latitudeCoordinates = info->GetLatitudeCoordinates();

  int numPointsPerCell = longitudeCoordinates->GetNumberOfComponents();
  vtkIdType totalNumCells = longitudeCoordinates->GetNumberOfTuples();

  double bounds[6];
  GetRangeOfAllComponents(longitudeCoordinates, bounds + 0);
  GetRangeOfAllComponents(latitudeCoordinates, bounds + 2);
  bounds[4] = bounds[5] = 0.0;

  VTK_CREATE(vtkPoints, points);
  points->SetDataTypeToDouble();
  points->Allocate(totalNumCells);

  VTK_CREATE(vtkMergePoints, locator);
  locator->InitPointInsertion(points, bounds);

  // Make space in output unstructured grid.
  unstructuredOutput->Allocate(extent[1] - extent[0]);
  vtkCellArray* cells = unstructuredOutput->GetCells();
  cells->AllocateEstimate(extent[1] - extent[0], numPointsPerCell);

  std::vector<vtkIdType> cellPoints(numPointsPerCell);

  // This is a rather lame way to break up cells amongst processes.  It will be
  // slow and ghost cells are totally screwed up.
  for (int cellId = extent[0]; cellId < extent[1]; cellId++)
  {
    for (int cellPointId = 0; cellPointId < numPointsPerCell; cellPointId++)
    {
      double coord[3];
      coord[0] = longitudeCoordinates->GetComponent(cellId, cellPointId);
      coord[1] = latitudeCoordinates->GetComponent(cellId, cellPointId);
      coord[2] = 0.0;

      vtkIdType pointId;
      locator->InsertUniquePoint(coord, pointId);

      cellPoints[cellPointId] = pointId;
    }
    unstructuredOutput->InsertNextCell(VTK_POLYGON, numPointsPerCell, &cellPoints.at(0));
  }

  points->Squeeze();
  unstructuredOutput->SetPoints(points);
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::AddUnstructuredSphericalCoordinates(
  vtkUnstructuredGrid* unstructuredOutput, const int extent[6])
{
  // First load the data as rectilinear coordinates, and then convert them
  // to spherical coordinates.  Not only does this reuse code, but it also
  // probably makes the locator more efficient this way.
  this->AddUnstructuredRectilinearCoordinates(unstructuredOutput, extent);

  double height = 1.0 * this->VerticalScale + this->VerticalBias;
  if (height <= 0.0)
  {
    height = 1.0;
  }

  vtkPoints* points = unstructuredOutput->GetPoints();
  vtkIdType numPoints = points->GetNumberOfPoints();
  for (vtkIdType pointId = 0; pointId < numPoints; pointId++)
  {
    double lonLat[3];
    points->GetPoint(pointId, lonLat);
    double lon = vtkMath::RadiansFromDegrees(lonLat[0]);
    double lat = vtkMath::RadiansFromDegrees(lonLat[1]);

    double cartesianCoord[3];
    cartesianCoord[0] = height * cos(lon) * cos(lat);
    cartesianCoord[1] = height * sin(lon) * cos(lat);
    cartesianCoord[2] = height * sin(lat);
    points->SetPoint(pointId, cartesianCoord);
  }
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::ReadMetaData(int ncFD)
{
  vtkDebugMacro("ReadMetaData");

  int numDimensions;
  CALL_NETCDF(nc_inq_ndims(ncFD, &numDimensions));
  this->DimensionInfo->v.resize(numDimensions);

  std::set<vtkStdString> specialVariables;

  for (int i = 0; i < numDimensions; i++)
  {
    this->DimensionInfo->v[i] = vtkDimensionInfo(ncFD, i);

    // Record any special variables for this dimension.
    vtkStringArray* dimensionVariables = this->DimensionInfo->v[i].GetSpecialVariables();
    for (vtkIdType j = 0; j < dimensionVariables->GetNumberOfValues(); j++)
    {
      specialVariables.insert(dimensionVariables->GetValue(j));
    }
  }

  int numVariables;
  CALL_NETCDF(nc_inq_nvars(ncFD, &numVariables));

  // Check all variables for special 2D coordinates.
  for (int i = 0; i < numVariables; i++)
  {
    vtkDependentDimensionInfo info(ncFD, i, this);
    if (!info.GetValid())
      continue;
    if (this->FindDependentDimensionInfo(info.GetGridDimensions()) != nullptr)
    {
      continue;
    }

    this->DependentDimensionInfo->v.push_back(info);

    // Record any special variables.
    vtkStringArray* dimensionVariables = info.GetSpecialVariables();
    for (vtkIdType j = 0; j < dimensionVariables->GetNumberOfValues(); j++)
    {
      specialVariables.insert(dimensionVariables->GetValue(j));
    }
  }

  // Look at all variables and record them so that the user can select which
  // ones he wants.  This oddness of adding and removing from
  // VariableArraySelection is to preserve any current settings for variables.
  typedef std::set<vtkStdString> stringSet;
  stringSet variablesToAdd;
  stringSet variablesToRemove;

  // Initialize variablesToRemove with all the variables.  Then remove them from
  // the list as we find them.
  for (int i = 0; i < this->VariableArraySelection->GetNumberOfArrays(); i++)
  {
    variablesToRemove.insert(this->VariableArraySelection->GetArrayName(i));
  }

  for (int i = 0; i < numVariables; i++)
  {
    char name[NC_MAX_NAME + 1];
    CALL_NETCDF(nc_inq_varname(ncFD, i, name));
    if (specialVariables.find(name) == specialVariables.end())
    {
      if (variablesToRemove.find(name) == variablesToRemove.end())
      {
        // Variable not already here.  Insert it in the variables to add.
        variablesToAdd.insert(name);
      }
      else
      {
        // Variable already exists.  Leave it be.  Remove it from the
        // variablesToRemove list.
        variablesToRemove.erase(name);
      }
    }
  }

  // Add and remove variables.  This will be a no-op if the variables have not
  // changed.
  for (stringSet::iterator removeItr = variablesToRemove.begin();
       removeItr != variablesToRemove.end(); ++removeItr)
  {
    this->VariableArraySelection->RemoveArrayByName(removeItr->c_str());
  }
  for (stringSet::iterator addItr = variablesToAdd.begin(); addItr != variablesToAdd.end();
       ++addItr)
  {
    this->VariableArraySelection->AddArray(addItr->c_str());
  }

  return 1;
}

//------------------------------------------------------------------------------
int vtkNetCDFCFReader::IsTimeDimension(int vtkNotUsed(ncFD), int dimId)
{
  return (this->GetDimensionInfo(dimId)->GetUnits() == vtkDimensionInfo::TIME_UNITS);
}

//------------------------------------------------------------------------------
vtkSmartPointer<vtkDoubleArray> vtkNetCDFCFReader::GetTimeValues(int vtkNotUsed(ncFD), int dimId)
{
  return this->GetDimensionInfo(dimId)->GetCoordinates();
}

//------------------------------------------------------------------------------
inline vtkNetCDFCFReader::vtkDimensionInfo* vtkNetCDFCFReader::GetDimensionInfo(int dimension)
{
  return &(this->DimensionInfo->v.at(dimension));
}

//------------------------------------------------------------------------------
vtkNetCDFCFReader::vtkDependentDimensionInfo* vtkNetCDFCFReader::FindDependentDimensionInfo(
  vtkIntArray* dims)
{
  for (size_t i = 0; i < this->DependentDimensionInfo->v.size(); i++)
  {
    vtkIntArray* dependentDims = this->DependentDimensionInfo->v[i].GetGridDimensions();
    if (dims->GetNumberOfTuples() == dependentDims->GetNumberOfTuples())
    {
      bool same = true;
      for (vtkIdType j = 0; j < dims->GetNumberOfTuples(); j++)
      {
        if (dims->GetValue(j) != dependentDims->GetValue(j))
        {
          same = false;
          break;
        }
      }
      if (same)
        return &(this->DependentDimensionInfo->v[i]);
    }
  }
  return nullptr;
}

//------------------------------------------------------------------------------
void vtkNetCDFCFReader::IdentifySphericalCoordinates(
  vtkIntArray* dimensions, int& longitudeDim, int& latitudeDim, int& verticalDim)
{
  longitudeDim = latitudeDim = verticalDim = -1;
  for (int i = 0; i < dimensions->GetNumberOfTuples(); i++)
  {
    switch (this->GetDimensionInfo(dimensions->GetValue(i))->GetUnits())
    {
      case vtkDimensionInfo::LONGITUDE_UNITS:
        longitudeDim = i;
        break;
      case vtkDimensionInfo::LATITUDE_UNITS:
        latitudeDim = i;
        break;
      case vtkDimensionInfo::UNDEFINED_UNITS:
      case vtkDimensionInfo::TIME_UNITS:
      case vtkDimensionInfo::VERTICAL_UNITS:
      default:
        verticalDim = i;
        break;
    }
  }
}

//------------------------------------------------------------------------------
vtkNetCDFCFReader::CoordinateTypesEnum vtkNetCDFCFReader::CoordinateType(vtkIntArray* dimensions)
{
  vtkDependentDimensionInfo* dependentDimInfo = this->FindDependentDimensionInfo(dimensions);

  // Check to see if using p-sided cells.
  if (dependentDimInfo && dependentDimInfo->GetCellsUnstructured())
  {
    if (this->SphericalCoordinates)
    {
      return COORDS_SPHERICAL_PSIDED_CELLS;
    }
    else
    {
      return COORDS_EUCLIDEAN_PSIDED_CELLS;
    }
  }

  // Check to see if using 4-sided cells.
  if (dependentDimInfo && !dependentDimInfo->GetCellsUnstructured() &&
    dependentDimInfo->GetHasBounds())
  {
    if (this->SphericalCoordinates)
    {
      return COORDS_SPHERICAL_4SIDED_CELLS;
    }
    else
    {
      return COORDS_EUCLIDEAN_4SIDED_CELLS;
    }
  }

  // Check to see if using 2D coordinate lookup.
  if (dependentDimInfo && !dependentDimInfo->GetCellsUnstructured() &&
    !dependentDimInfo->GetHasBounds())
  {
    if (this->SphericalCoordinates)
    {
      return COORDS_2D_SPHERICAL;
    }
    else
    {
      return COORDS_2D_EUCLIDEAN;
    }
  }

  // Check to see if we should (otherwise) be using spherical coordinates.
  if (this->SphericalCoordinates)
  {
    int longitudeDim, latitudeDim, verticalDim;
    this->IdentifySphericalCoordinates(dimensions, longitudeDim, latitudeDim, verticalDim);
    if ((longitudeDim != -1) && (latitudeDim != -1) &&
      ((dimensions->GetNumberOfTuples() == 2) || (verticalDim != -1)))
    {
      return COORDS_REGULAR_SPHERICAL;
    }
  }

  // Check to see if any dimension as irregular spacing.
  for (int i = 0; i < dimensions->GetNumberOfTuples(); i++)
  {
    int dimId = dimensions->GetValue(i);
    if (!this->GetDimensionInfo(dimId)->GetHasRegularSpacing())
    {
      return COORDS_NONUNIFORM_RECTILINEAR;
    }
  }

  // All dimensions appear to be uniform rectilinear.
  return COORDS_UNIFORM_RECTILINEAR;
}

//------------------------------------------------------------------------------
bool vtkNetCDFCFReader::DimensionsAreForPointData(vtkIntArray* dimensions)
{
  switch (this->CoordinateType(dimensions))
  {
    case COORDS_UNIFORM_RECTILINEAR:
      return true;
    case COORDS_NONUNIFORM_RECTILINEAR:
      return true;
    case COORDS_REGULAR_SPHERICAL:
      return false;
    case COORDS_2D_EUCLIDEAN:
      return true;
    case COORDS_2D_SPHERICAL:
      return true;
    case COORDS_EUCLIDEAN_4SIDED_CELLS:
      return false;
    case COORDS_SPHERICAL_4SIDED_CELLS:
      return false;
    case COORDS_EUCLIDEAN_PSIDED_CELLS:
      return false;
    case COORDS_SPHERICAL_PSIDED_CELLS:
      return false;
    default:
      vtkErrorMacro("Internal error: unknown coordinate type.");
      return true;
  }
}