text: /*=========================================================================
  
  Program:   Visualization Toolkit
  Module:    vtkNetCDFPOPReader.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.
  
  =========================================================================*/
  
  #include "vtkNetCDFPOPReader.h"
  #include "vtkCallbackCommand.h"
  #include "vtkDataArraySelection.h"
  #include "vtkFloatArray.h"
  #include "vtkInformation.h"
  #include "vtkInformationVector.h"
  #include "vtkObjectFactory.h"
  #include "vtkPointData.h"
  #include "vtkRectilinearGrid.h"
  #include "vtkSmartPointer.h"
  #include "vtkStreamingDemandDrivenPipeline.h"
  
  #include "vtk_netcdf.h"
  #include <string>
  #include <vector>
  
  vtkStandardNewMacro(vtkNetCDFPOPReader);
  
  //============================================================================
  #define CALL_NETCDF(call)                                                                          \
    {                                                                                                \
      int errorcode = call;                                                                          \
      if (errorcode != NC_NOERR)                                                                     \
      {                                                                                              \
        vtkErrorMacro(<< "netCDF Error: " << nc_strerror(errorcode));                                \
        return 0;                                                                                    \
      }                                                                                              \
    }
  //============================================================================
  
  class vtkNetCDFPOPReaderInternal
  {
  public:
    vtkSmartPointer<vtkDataArraySelection> VariableArraySelection;
    // a mapping from the list of all variables to the list of available
    // point-based variables
    std::vector<int> VariableMap;
    vtkNetCDFPOPReaderInternal()
    {
      this->VariableArraySelection = vtkSmartPointer<vtkDataArraySelection>::New();
    }
  };
  
  //------------------------------------------------------------------------------
  // set default values
  vtkNetCDFPOPReader::vtkNetCDFPOPReader()
  {
    this->SetNumberOfInputPorts(0);
    this->SetNumberOfOutputPorts(1);
    this->FileName = nullptr;
    this->NCDFFD = 0;
    this->OpenedFileName = nullptr;
    this->Stride[0] = this->Stride[1] = this->Stride[2] = 1;
    this->SelectionObserver = vtkCallbackCommand::New();
    this->SelectionObserver->SetCallback(&vtkNetCDFPOPReader::SelectionModifiedCallback);
    this->SelectionObserver->SetClientData(this);
    this->Internals = new vtkNetCDFPOPReaderInternal;
    this->Internals->VariableArraySelection->AddObserver(
      vtkCommand::ModifiedEvent, this->SelectionObserver);
  }
  
  //------------------------------------------------------------------------------
  // delete filename and netcdf file descriptor
  vtkNetCDFPOPReader::~vtkNetCDFPOPReader()
  {
    this->SetFileName(nullptr);
    if (this->OpenedFileName)
    {
      nc_close(this->NCDFFD);
      this->SetOpenedFileName(nullptr);
    }
    if (this->SelectionObserver)
    {
      this->SelectionObserver->Delete();
      this->SelectionObserver = nullptr;
    }
    delete this->Internals;
    this->Internals = nullptr;
  }
  
  //------------------------------------------------------------------------------
  void vtkNetCDFPOPReader::PrintSelf(ostream& os, vtkIndent indent)
  {
    this->Superclass::PrintSelf(os, indent);
  
    os << indent << "FileName: " << (this->FileName ? this->FileName : "(nullptr)") << endl;
    os << indent << "OpenedFileName: " << (this->OpenedFileName ? this->OpenedFileName : "(nullptr)")
       << endl;
    os << indent << "Stride: {" << this->Stride[0] << ", " << this->Stride[1] << ", "
       << this->Stride[2] << ", "
       << "}" << endl;
    os << indent << "NCDFFD: " << this->NCDFFD << endl;
  
    this->Internals->VariableArraySelection->PrintSelf(os, indent.GetNextIndent());
  }
  
  //------------------------------------------------------------------------------
  // RequestInformation supplies global meta information
  // This should return the reality of what the reader is going to supply.
  // This retrieve the extents for the rectilinear grid
  // NC_MAX_VAR_DIMS comes from the nc library
  int vtkNetCDFPOPReader::RequestInformation(vtkInformation* vtkNotUsed(request),
    vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
  {
    vtkInformation* outInfo = outputVector->GetInformationObject(0);
    if (this->FileName == nullptr)
    {
      vtkErrorMacro("FileName not set.");
      return 0;
    }
  
    if (this->OpenedFileName == nullptr || strcmp(this->OpenedFileName, this->FileName) != 0)
    {
      if (this->OpenedFileName)
      {
        nc_close(this->NCDFFD);
        this->SetOpenedFileName(nullptr);
      }
      int retval = nc_open(this->FileName, NC_NOWRITE, &this->NCDFFD); // read file
      if (retval != NC_NOERR)                                          // checks if read file error
      {
        // we don't need to close the file if there was an error opening the file
        vtkErrorMacro(<< "Can't read file " << nc_strerror(retval));
        return 0;
      }
      this->SetOpenedFileName(this->FileName);
    }
    // get number of variables from file
    int numberOfVariables;
    nc_inq_nvars(this->NCDFFD, &numberOfVariables);
    int dimidsp[NC_MAX_VAR_DIMS];
    int extent[6];
    int dataDimension;
    size_t dimensions[4]; // dimension value
    this->Internals->VariableMap.resize(numberOfVariables);
    char variableName[NC_MAX_NAME + 1];
    int actualVariableCounter = 0;
    // For every variable in the file
    for (int i = 0; i < numberOfVariables; i++)
    {
      this->Internals->VariableMap[i] = -1;
      // get number of dimensions
      CALL_NETCDF(nc_inq_varndims(this->NCDFFD, i, &dataDimension));
      // Variable Dimension ID's containing x,y,z coords for the rectilinear
      // grid spacing
      CALL_NETCDF(nc_inq_vardimid(this->NCDFFD, i, dimidsp));
      if (dataDimension == 3)
      {
        this->Internals->VariableMap[i] = actualVariableCounter++;
        // get variable name
        CALL_NETCDF(nc_inq_varname(this->NCDFFD, i, variableName));
        this->Internals->VariableArraySelection->AddArray(variableName);
        for (int m = 0; m < dataDimension; m++)
        {
          CALL_NETCDF(nc_inq_dimlen(this->NCDFFD, dimidsp[m], dimensions + m));
          // acquire variable dimensions
        }
        extent[0] = extent[2] = extent[4] = 0; // set extent
        extent[1] = static_cast<int>((dimensions[2] - 1) / this->Stride[0]);
        extent[3] = static_cast<int>((dimensions[1] - 1) / this->Stride[1]);
        extent[5] = static_cast<int>((dimensions[0] - 1) / this->Stride[2]);
      }
    }
    // fill in the extent information
    outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent, 6);
    return 1;
  }
  
  //------------------------------------------------------------------------------
  // Setting extents of the rectilinear grid
  int vtkNetCDFPOPReader::RequestData(vtkInformation* request,
    vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
  {
    this->UpdateProgress(0);
    // the default implementation is to do what the old pipeline did find what
    // output is requesting the data, and pass that into ExecuteData
    // which output port did the request come from
    int outputPort = request->Get(vtkDemandDrivenPipeline::FROM_OUTPUT_PORT());
    // if output port is negative then that means this filter is calling the
    // update directly, in that case just assume port 0
    if (outputPort == -1)
    {
      outputPort = 0;
    }
    // get the data object
    vtkInformation* outInfo = outputVector->GetInformationObject(outputPort);
  
    vtkDataObject* output = outInfo->Get(vtkDataObject::DATA_OBJECT());
    int subext[6];
    // vtkInformation * outInfo = output->GetInformationObject(0);
    outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), subext);
    vtkRectilinearGrid* rgrid = vtkRectilinearGrid::SafeDownCast(output);
    rgrid->SetExtent(subext);
    // setup extents for netcdf library to read the netcdf data file
    size_t start[] = { static_cast<size_t>(subext[4] * this->Stride[2]),
      static_cast<size_t>(subext[2] * this->Stride[1]),
      static_cast<size_t>(subext[0] * this->Stride[0]) };
  
    size_t count[] = { static_cast<size_t>(subext[5] - subext[4] + 1),
      static_cast<size_t>(subext[3] - subext[2] + 1),
      static_cast<size_t>(subext[1] - subext[0] + 1) };
  
    ptrdiff_t rStride[3] = { (ptrdiff_t)this->Stride[2], (ptrdiff_t)this->Stride[1],
      (ptrdiff_t)this->Stride[0] };
  
    /*
      vtkMultiThreaderIDType pid = vtkMultiThreader::GetCurrentThreadID();
      //The getpid() and getppid() functions are always successful, and no
      //return value is reserved to indicate an error. So is this check
      //necessary?
      if ((pid = vtkMultiThreader::GetCurrentThreadID()) == 0)
        {
        cerr<< "unable to get pid" << endl;
        }
      else
        {
        cerr << "The process id is " << pid << endl;
       }
  
       cerr << "subext " << subext[0] << " " << subext[1]
            << " " << subext[2] << " " << subext[3] << " "
             << subext[4] << " " << subext[5] << endl;
    */
    // initialize memory (raw data space, x y z axis space) and rectilinear grid
    bool firstPass = true;
    for (size_t i = 0; i < this->Internals->VariableMap.size(); i++)
    {
      if (this->Internals->VariableMap[i] != -1 &&
        this->Internals->VariableArraySelection->GetArraySetting(this->Internals->VariableMap[i]) !=
          0)
      {
        // varidp is probably i in which case nc_inq_varid isn't needed
        int varidp;
        nc_inq_varid(this->NCDFFD,
          this->Internals->VariableArraySelection->GetArrayName(this->Internals->VariableMap[i]),
          &varidp);
  
        if (firstPass)
        {
          int dimidsp[3];
          nc_inq_vardimid(this->NCDFFD, varidp, dimidsp);
          firstPass = false;
          float* x = new float[count[0]];
          float* y = new float[count[1]];
          float* z = new float[count[2]];
          // gets data from x,y,z axis (spacing)
  #if 0
          nc_get_vara_float(this->NCDFFD, dimidsp[0], start, count, x);
          nc_get_vara_float(this->NCDFFD, dimidsp[1], start+1, count+1, y);
          nc_get_vara_float(this->NCDFFD, dimidsp[2], start+2, count+2, z);
  #else
          nc_get_vars_float(this->NCDFFD, dimidsp[0], start, count, rStride, x);
          nc_get_vars_float(this->NCDFFD, dimidsp[1], start + 1, count + 1, rStride + 1, y);
          nc_get_vars_float(this->NCDFFD, dimidsp[2], start + 2, count + 2, rStride + 2, z);
  #endif
          vtkFloatArray* xCoords = vtkFloatArray::New();
          xCoords->SetArray(z, static_cast<vtkIdType>(count[2]), 0, 1);
          vtkFloatArray* yCoords = vtkFloatArray::New();
          yCoords->SetArray(y, static_cast<vtkIdType>(count[1]), 0, 1);
          for (unsigned int q = 0; q < count[0]; q++)
          {
            x[q] = -x[q];
          }
          vtkFloatArray* zCoords = vtkFloatArray::New();
          zCoords->SetArray(x, static_cast<vtkIdType>(count[0]), 0, 1);
          rgrid->SetXCoordinates(xCoords);
          rgrid->SetYCoordinates(yCoords);
          rgrid->SetZCoordinates(zCoords);
          xCoords->Delete();
          yCoords->Delete();
          zCoords->Delete();
        }
        // create vtkFloatArray and get the scalars into it
        vtkFloatArray* scalars = vtkFloatArray::New();
        vtkIdType numberOfTuples = static_cast<vtkIdType>((count[0]) * (count[1]) * (count[2]));
        float* data = new float[numberOfTuples];
  #if 0
        nc_get_vara_float(this->NCDFFD, varidp, start, count, data);
  #else
        nc_get_vars_float(this->NCDFFD, varidp, start, count, rStride, data);
  #endif
        scalars->SetArray(data, numberOfTuples, 0, 1);
        // set list of variables to display data on rectilinear grid
        scalars->SetName(
          this->Internals->VariableArraySelection->GetArrayName(this->Internals->VariableMap[i]));
        rgrid->GetPointData()->AddArray(scalars);
        scalars->Delete();
      }
      this->UpdateProgress((i + 1.0) / this->Internals->VariableMap.size());
    }
    return 1;
  }
  
  //------------------------------------------------------------------------------
  // following 5 functions are used for paraview user interface
  void vtkNetCDFPOPReader::SelectionModifiedCallback(
    vtkObject*, unsigned long, void* clientdata, void*)
  {
    static_cast<vtkNetCDFPOPReader*>(clientdata)->Modified();
  }
  
  //------------------------------------------------------------------------------
  int vtkNetCDFPOPReader::GetNumberOfVariableArrays()
  {
    return this->Internals->VariableArraySelection->GetNumberOfArrays();
  }
  
  //------------------------------------------------------------------------------
  const char* vtkNetCDFPOPReader::GetVariableArrayName(int index)
  {
    if (index < 0 || index >= this->GetNumberOfVariableArrays())
    {
      return nullptr;
    }
    return this->Internals->VariableArraySelection->GetArrayName(index);
  }
  
  //------------------------------------------------------------------------------
  int vtkNetCDFPOPReader::GetVariableArrayStatus(const char* name)
  {
    return this->Internals->VariableArraySelection->ArrayIsEnabled(name);
  }
  
  //------------------------------------------------------------------------------
  void vtkNetCDFPOPReader::SetVariableArrayStatus(const char* name, int status)
  {
    vtkDebugMacro("Set cell array \"" << name << "\" status to: " << status);
    if (this->Internals->VariableArraySelection->ArrayExists(name) == 0)
    {
      vtkErrorMacro(<< name << " is not available in the file.");
      return;
    }
    int enabled = this->Internals->VariableArraySelection->ArrayIsEnabled(name);
    if (status != 0 && enabled == 0)
    {
      this->Internals->VariableArraySelection->EnableArray(name);
      this->Modified();
    }
    else if (status == 0 && enabled != 0)
    {
      this->Internals->VariableArraySelection->DisableArray(name);
      this->Modified();
    }
  }
  