BUG: push bug fix from Ken Mooreland.
[VTK.git] / IO / vtkNetCDFReader.cxx
1 // -*- c++ -*-
2 /*=========================================================================
3
4   Program:   Visualization Toolkit
5   Module:    vtkNetCDFReader.cxx
6
7   Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
8   All rights reserved.
9   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
10
11      This software is distributed WITHOUT ANY WARRANTY; without even
12      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
13      PURPOSE.  See the above copyright notice for more information.
14
15 =========================================================================*/
16
17 /*-------------------------------------------------------------------------
18   Copyright 2008 Sandia Corporation.
19   Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
20   the U.S. Government retains certain rights in this software.
21 -------------------------------------------------------------------------*/
22
23 #include "vtkNetCDFReader.h"
24
25 #include "vtkCallbackCommand.h"
26 #include "vtkCellData.h"
27 #include "vtkDataArraySelection.h"
28 #include "vtkDoubleArray.h"
29 #include "vtkImageData.h"
30 #include "vtkInformation.h"
31 #include "vtkInformationVector.h"
32 #include "vtkIntArray.h"
33 #include "vtkMath.h"
34 #include "vtkObjectFactory.h"
35 #include "vtkPointData.h"
36 #include "vtkRectilinearGrid.h"
37 #include "vtkStreamingDemandDrivenPipeline.h"
38 #include "vtkStringArray.h"
39 #include "vtkStructuredGrid.h"
40
41 #include "vtkSmartPointer.h"
42 #define VTK_CREATE(type, name) \
43   vtkSmartPointer<type> name = vtkSmartPointer<type>::New()
44
45 #include <vtkstd/algorithm>
46 #include <vtksys/SystemTools.hxx>
47
48 #include <netcdf.h>
49
50 #define CALL_NETCDF(call) \
51   { \
52     int errorcode = call; \
53     if (errorcode != NC_NOERR) \
54       { \
55       vtkErrorMacro(<< "netCDF Error: " << nc_strerror(errorcode)); \
56       return 0; \
57       } \
58   }
59
60 #include <ctype.h>
61
62 //=============================================================================
63 static int NetCDFTypeToVTKType(nc_type type)
64 {
65   switch (type)
66     {
67     case NC_BYTE: return VTK_UNSIGNED_CHAR;
68     case NC_CHAR: return VTK_CHAR;
69     case NC_SHORT: return VTK_SHORT;
70     case NC_INT: return VTK_INT;
71     case NC_FLOAT: return VTK_FLOAT;
72     case NC_DOUBLE: return VTK_DOUBLE;
73     default:
74       vtkGenericWarningMacro(<< "Unknown netCDF variable type "
75                              << type);
76       return -1;
77     }
78 }
79
80 //=============================================================================
81 vtkCxxRevisionMacro(vtkNetCDFReader, "1.7.4.1");
82 vtkStandardNewMacro(vtkNetCDFReader);
83
84 //-----------------------------------------------------------------------------
85 vtkNetCDFReader::vtkNetCDFReader()
86 {
87   this->SetNumberOfInputPorts(0);
88
89   this->FileName = NULL;
90   this->ReplaceFillValueWithNan = 0;
91
92   this->LoadingDimensions = vtkSmartPointer<vtkIntArray>::New();
93
94   this->VariableArraySelection = vtkSmartPointer<vtkDataArraySelection>::New();
95   VTK_CREATE(vtkCallbackCommand, cbc);
96   cbc->SetCallback(&vtkNetCDFReader::SelectionModifiedCallback);
97   cbc->SetClientData(this);
98   this->VariableArraySelection->AddObserver(vtkCommand::ModifiedEvent, cbc);
99
100   this->VariableDimensions = vtkStringArray::New();
101   this->AllDimensions = vtkStringArray::New();
102 }
103
104 vtkNetCDFReader::~vtkNetCDFReader()
105 {
106   this->SetFileName(NULL);
107   this->VariableDimensions->Delete();
108   this->AllDimensions->Delete();
109 }
110
111 void vtkNetCDFReader::PrintSelf(ostream &os, vtkIndent indent)
112 {
113   this->Superclass::PrintSelf(os, indent);
114
115   os << indent << "FileName: "
116      << (this->FileName ? this->FileName : "(NULL)") << endl;
117   os << indent << "ReplaceFillValueWithNan: "
118      << this->ReplaceFillValueWithNan << endl;
119
120   os << indent << "VariableArraySelection:" << endl;
121   this->VariableArraySelection->PrintSelf(os, indent.GetNextIndent());
122   os << indent << "VariableDimensions: " << this->VariableDimensions << endl;
123   os << indent << "AllDimensions: " << this->AllDimensions << endl;
124 }
125
126 //-----------------------------------------------------------------------------
127 int vtkNetCDFReader::RequestDataObject(
128                                  vtkInformation *vtkNotUsed(request),
129                                  vtkInformationVector **vtkNotUsed(inputVector),
130                                  vtkInformationVector *outputVector)
131 {
132   vtkInformation *outInfo = outputVector->GetInformationObject(0);
133   vtkDataObject *output = vtkDataObject::GetData(outInfo);
134
135   if (!output || !output->IsA("vtkImageData"))
136     {
137     output = vtkImageData::New();
138     output->SetPipelineInformation(outInfo);
139     output->Delete();   // Not really deleted.
140     }
141
142   return 1;
143 }
144
145 //-----------------------------------------------------------------------------
146 int vtkNetCDFReader::RequestInformation(
147                                  vtkInformation *vtkNotUsed(request),
148                                  vtkInformationVector **vtkNotUsed(inputVector),
149                                  vtkInformationVector *outputVector)
150 {
151   if (!this->UpdateMetaData()) return 0;
152
153   vtkInformation *outInfo = outputVector->GetInformationObject(0);
154
155   int ncFD;
156   CALL_NETCDF(nc_open(this->FileName, NC_NOWRITE, &ncFD));
157
158   VTK_CREATE(vtkDoubleArray, timeValues);
159   VTK_CREATE(vtkIntArray, currentDimensions);
160   this->LoadingDimensions->Initialize();
161   int numArrays = this->VariableArraySelection->GetNumberOfArrays();
162   int numDims = 0;
163   for (int arrayIndex = 0; arrayIndex < numArrays; arrayIndex++)
164     {
165     if (!this->VariableArraySelection->GetArraySetting(arrayIndex)) continue;
166
167     const char *name = this->VariableArraySelection->GetArrayName(arrayIndex);
168     int varId;
169     CALL_NETCDF(nc_inq_varid(ncFD, name, &varId));
170
171     int currentNumDims;
172     CALL_NETCDF(nc_inq_varndims(ncFD, varId, &currentNumDims));
173     if (currentNumDims < 1) continue;
174     currentDimensions->SetNumberOfComponents(1);
175     currentDimensions->SetNumberOfTuples(currentNumDims);
176     CALL_NETCDF(nc_inq_vardimid(ncFD, varId,
177                                 currentDimensions->GetPointer(0)));
178
179     // Assumption: time dimension is first.
180     int timeDim = currentDimensions->GetValue(0);       // Not determined yet.
181     if (this->IsTimeDimension(ncFD, timeDim))
182       {
183       // Get time step information.
184       VTK_CREATE(vtkDoubleArray, compositeTimeValues);
185       vtkSmartPointer<vtkDoubleArray> currentTimeValues
186         = this->GetTimeValues(ncFD, timeDim);
187       double *oldTime = timeValues->GetPointer(0);
188       double *newTime = currentTimeValues->GetPointer(0);
189       double *oldTimeEnd = oldTime + timeValues->GetNumberOfTuples();
190       double *newTimeEnd = newTime + currentTimeValues->GetNumberOfTuples();
191       compositeTimeValues->Allocate(  timeValues->GetNumberOfTuples()
192                                     + currentTimeValues->GetNumberOfTuples());
193       compositeTimeValues->SetNumberOfComponents(1);
194       while ((oldTime < oldTimeEnd) || (newTime < newTimeEnd))
195         {
196         if (   (newTime >= newTimeEnd)
197             || ((oldTime < oldTimeEnd) && (*oldTime < *newTime)) )
198           {
199           compositeTimeValues->InsertNextTuple1(*oldTime);
200           oldTime++;
201           }
202         else if ((oldTime >= oldTimeEnd) || (*newTime < *oldTime))
203           {
204           compositeTimeValues->InsertNextTuple1(*newTime);
205           newTime++;
206           }
207         else // *oldTime == *newTime
208           {
209           compositeTimeValues->InsertNextTuple1(*oldTime);
210           oldTime++;  newTime++;
211           }
212         }
213       timeValues = compositeTimeValues;
214
215       // Strip off time dimension from what we load (we will use it to
216       // subset instead).
217       currentDimensions->RemoveTuple(0);
218       currentNumDims--;
219       }
220
221     // Remember the first variable we encounter.  Use it to determine extents
222     // (below).
223     if (numDims == 0)
224       {
225       numDims = currentNumDims;
226       this->LoadingDimensions->DeepCopy(currentDimensions);
227       }
228     }
229
230   // Using the extent information (captured partially in
231   // this->LoadingDimensions) report extents.
232   vtkDataObject *output = vtkDataObject::GetData(outInfo);
233   if (output && (output->GetExtentType() == VTK_3D_EXTENT))
234     {
235     bool pointData = this->DimensionsAreForPointData(
236                                   this->LoadingDimensions->GetPointer(0),
237                                   this->LoadingDimensions->GetNumberOfTuples());
238     int extent[6];
239     for (int i = 0 ; i < 3; i++)
240       {
241       extent[2*i] = 0;
242       if (i < this->LoadingDimensions->GetNumberOfTuples())
243         {
244         size_t dimlength;
245         // Remember that netCDF arrays are indexed backward from VTK images.
246         int dim = this->LoadingDimensions->GetValue(numDims-i-1);
247         CALL_NETCDF(nc_inq_dimlen(ncFD, dim, &dimlength));
248         extent[2*i+1] = static_cast<int>(dimlength-1);
249         // For cell data, add one to the extent (which is for points).
250         if (!pointData) extent[2*i+1]++;
251         }
252       else
253         {
254         extent[2*i+1] = 0;
255         }
256       }
257     vtkDebugMacro(<< "Whole extents: "
258                   << extent[0] << ", " << extent[1] << ", "
259                   << extent[2] << ", " << extent[3] << ", "
260                   << extent[4] << ", " << extent[5]);
261     outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent, 6);
262     }
263
264   // If we have time, report that.
265   if (timeValues && (timeValues->GetNumberOfTuples() > 0))
266     {
267     outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
268                  timeValues->GetPointer(0),
269                  timeValues->GetNumberOfTuples());
270     double timeRange[2];
271     timeRange[0] = timeValues->GetValue(0);
272     timeRange[1] = timeValues->GetValue(timeValues->GetNumberOfTuples()-1);
273     outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2);
274     }
275   else
276     {
277     outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
278     outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
279     }
280
281   CALL_NETCDF(nc_close(ncFD));
282
283   return 1;
284 }
285
286 //-----------------------------------------------------------------------------
287 int vtkNetCDFReader::RequestData(vtkInformation *vtkNotUsed(request),
288                                  vtkInformationVector **vtkNotUsed(inputVector),
289                                  vtkInformationVector *outputVector)
290 {
291   vtkInformation *outInfo = outputVector->GetInformationObject(0);
292   // If the output is not a vtkDataSet, then the subclass needs to override
293   // this method.
294   vtkDataSet *output = vtkDataSet::GetData(outInfo);
295   if (!output)
296     {
297     vtkErrorMacro(<< "Bad output type.");
298     return 0;
299     }
300
301   // Set up the extent for regular-grid type data sets.
302   vtkImageData *imageOutput = vtkImageData::SafeDownCast(output);
303   vtkRectilinearGrid *rectOutput = vtkRectilinearGrid::SafeDownCast(output);
304   vtkStructuredGrid *structOutput = vtkStructuredGrid::SafeDownCast(output);
305   if (imageOutput)
306     {
307     imageOutput->SetExtent(imageOutput->GetUpdateExtent());
308     }
309   else if (rectOutput)
310     {
311     rectOutput->SetExtent(rectOutput->GetUpdateExtent());
312     }
313   else if (structOutput)
314     {
315     structOutput->SetExtent(structOutput->GetUpdateExtent());
316     }
317   else
318     {
319     // Superclass should handle extent setup if necessary.
320     }
321
322   // Get requested time step.
323   double time = 0.0;
324   if (outInfo->Has(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS()))
325     {
326     time
327       = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEPS())[0];
328     }
329
330   int ncFD;
331   CALL_NETCDF(nc_open(this->FileName, NC_NOWRITE, &ncFD));
332
333   // Iterate over arrays and load selected ones.
334   int numArrays = this->VariableArraySelection->GetNumberOfArrays();
335   for (int arrayIndex = 0; arrayIndex < numArrays; arrayIndex++)
336     {
337     if (!this->VariableArraySelection->GetArraySetting(arrayIndex)) continue;
338
339     const char *name = this->VariableArraySelection->GetArrayName(arrayIndex);
340
341     if (!this->LoadVariable(ncFD, name, time, output)) return 0;
342     }
343
344   CALL_NETCDF(nc_close(ncFD));
345
346   return 1;
347 }
348
349 //-----------------------------------------------------------------------------
350 void vtkNetCDFReader::SetFileName(const char *filename)
351 {
352   if (this->FileName == filename) return;
353   if (this->FileName && filename && (strcmp(this->FileName, filename) == 0))
354     {
355     return;
356     }
357
358   if (this->FileName)
359     {
360     delete[] this->FileName;
361     this->FileName = NULL;
362     }
363
364   if (filename)
365     {
366     this->FileName = new char[strlen(filename)+1];
367     strcpy(this->FileName, filename);
368     }
369
370   this->Modified();
371   this->FileNameMTime.Modified();
372 }
373
374 //----------------------------------------------------------------------------
375 void vtkNetCDFReader::SelectionModifiedCallback(vtkObject*, unsigned long,
376                                                 void* clientdata, void*)
377 {
378   static_cast<vtkNetCDFReader*>(clientdata)->Modified();
379 }
380
381 //-----------------------------------------------------------------------------
382 int vtkNetCDFReader::GetNumberOfVariableArrays()
383 {
384   return this->VariableArraySelection->GetNumberOfArrays();
385 }
386
387 //-----------------------------------------------------------------------------
388 const char* vtkNetCDFReader::GetVariableArrayName(int index)
389 {
390   return this->VariableArraySelection->GetArrayName(index);
391 }
392
393 //-----------------------------------------------------------------------------
394 int vtkNetCDFReader::GetVariableArrayStatus(const char* name)
395 {
396   return this->VariableArraySelection->ArrayIsEnabled(name);
397 }
398
399 //-----------------------------------------------------------------------------
400 void vtkNetCDFReader::SetVariableArrayStatus(const char* name, int status)
401 {
402   vtkDebugMacro("Set cell array \"" << name << "\" status to: " << status);
403   if(status)
404     {
405     this->VariableArraySelection->EnableArray(name);
406     }
407   else
408     {
409     this->VariableArraySelection->DisableArray(name);
410     }
411 }
412
413 //-----------------------------------------------------------------------------
414 void vtkNetCDFReader::SetDimensions(const char *dimensions)
415 {
416   this->VariableArraySelection->DisableAllArrays();
417
418   for (vtkIdType i = 0; i < this->VariableDimensions->GetNumberOfValues(); i++)
419     {
420     if (this->VariableDimensions->GetValue(i) == dimensions)
421       {
422       const char *variableName = this->VariableArraySelection->GetArrayName(i);
423       this->VariableArraySelection->EnableArray(variableName);
424       }
425     }
426 }
427
428 //-----------------------------------------------------------------------------
429 int vtkNetCDFReader::UpdateMetaData()
430 {
431   if (this->MetaDataMTime < this->FileNameMTime)
432     {
433     if (!this->FileName)
434       {
435       vtkErrorMacro(<< "FileName not set.");
436       return 0;
437       }
438
439     int ncFD;
440     CALL_NETCDF(nc_open(this->FileName, NC_NOWRITE, &ncFD));
441
442     int retval = this->ReadMetaData(ncFD);
443
444     if (retval) retval = this->FillVariableDimensions(ncFD);
445
446     if (retval) this->MetaDataMTime.Modified();
447
448     CALL_NETCDF(nc_close(ncFD));
449
450     return retval;
451     }
452   else
453     {
454     return 1;
455     }
456 }
457
458 //-----------------------------------------------------------------------------
459 vtkStdString vtkNetCDFReader::DescribeDimensions(int ncFD,
460                                                  const int *dimIds, int numDims)
461 {
462   vtkStdString description;
463   for (int i = 0; i < numDims; i++)
464     {
465     char name[NC_MAX_NAME+1];
466     CALL_NETCDF(nc_inq_dimname(ncFD, dimIds[i], name));
467     if (i > 0) description += " ";
468     description += name;
469     }
470   return description;
471 }
472
473 //-----------------------------------------------------------------------------
474 int vtkNetCDFReader::ReadMetaData(int ncFD)
475 {
476   int i;
477
478   vtkDebugMacro("ReadMetaData");
479
480   // Look at all variables and record them so that the user can select
481   // which ones he wants.
482   this->VariableArraySelection->RemoveAllArrays();
483
484   int numVariables;
485   CALL_NETCDF(nc_inq_nvars(ncFD, &numVariables));
486
487   for (i = 0; i < numVariables; i++)
488     {
489     char name[NC_MAX_NAME+1];
490     CALL_NETCDF(nc_inq_varname(ncFD, i, name));
491     this->VariableArraySelection->AddArray(name);
492     }
493
494   return 1;
495 }
496
497 //-----------------------------------------------------------------------------
498 int vtkNetCDFReader::FillVariableDimensions(int ncFD)
499 {
500   int numVar = this->GetNumberOfVariableArrays();
501   this->VariableDimensions->SetNumberOfValues(numVar);
502   this->AllDimensions->SetNumberOfValues(0);
503
504   for (int i = 0; i < numVar; i++)
505     {
506     // Get the dimensions of this variable and encode them in a string.
507     const char *varName = this->GetVariableArrayName(i);
508     int varId, numDim, dimIds[NC_MAX_VAR_DIMS];
509     CALL_NETCDF(nc_inq_varid(ncFD, varName, &varId));
510     CALL_NETCDF(nc_inq_varndims(ncFD, varId, &numDim));
511     CALL_NETCDF(nc_inq_vardimid(ncFD, varId, dimIds));
512     vtkStdString dimEncoding("(");
513     for (int j = 0; j < numDim; j++)
514       {
515       if ((j == 0) && (this->IsTimeDimension(ncFD, dimIds[j]))) continue;
516       char dimName[NC_MAX_NAME+1];
517       CALL_NETCDF(nc_inq_dimname(ncFD, dimIds[j], dimName));
518       if (dimEncoding.size() > 1) dimEncoding += ", ";
519       dimEncoding += dimName;
520       }
521     dimEncoding += ")";
522     // Store all dimensions in the VariableDimensions array.
523     this->VariableDimensions->SetValue(i, dimEncoding.c_str());
524     // Store unique dimensions in the AllDimensions array.
525     bool unique = true;
526     for (vtkIdType j = 0; j < this->AllDimensions->GetNumberOfValues(); j++)
527       {
528       if (this->AllDimensions->GetValue(j) == dimEncoding)
529         {
530         unique = false;
531         break;
532         }
533       }
534     if (unique) this->AllDimensions->InsertNextValue(dimEncoding);
535     }
536
537   return 1;
538 }
539
540 //-----------------------------------------------------------------------------
541 int vtkNetCDFReader::IsTimeDimension(int ncFD, int dimId)
542 {
543   char name[NC_MAX_NAME+1];
544   CALL_NETCDF(nc_inq_dimname(ncFD, dimId, name));
545   name[4] = '\0';       // Truncate to 4 characters.
546   return (vtksys::SystemTools::Strucmp(name, "time") == 0);
547 }
548
549 //-----------------------------------------------------------------------------
550 vtkSmartPointer<vtkDoubleArray> vtkNetCDFReader::GetTimeValues(int ncFD,
551                                                                int dimId)
552 {
553   VTK_CREATE(vtkDoubleArray, timeValues);
554   size_t dimLength;
555   CALL_NETCDF(nc_inq_dimlen(ncFD, dimId, &dimLength));
556   timeValues->SetNumberOfComponents(1);
557   timeValues->SetNumberOfTuples(dimLength);
558   for (size_t j = 0; j < dimLength; j++)
559     {
560     timeValues->SetValue(j, static_cast<double>(j));
561     }
562   return timeValues;
563 }
564
565 //-----------------------------------------------------------------------------
566 int vtkNetCDFReader::LoadVariable(int ncFD, const char *varName, double time,
567                                   vtkDataSet *output)
568 {
569   // Get the variable id.
570   int varId;
571   CALL_NETCDF(nc_inq_varid(ncFD, varName, &varId));
572
573   // Get dimension info.
574   int numDims;
575   CALL_NETCDF(nc_inq_varndims(ncFD, varId, &numDims));
576   if (numDims > 4)
577     {
578     vtkErrorMacro(<< "More than 3 dims + time not supported in variable "
579                   << varName);
580     return 0;
581     }
582   int dimIds[4];
583   CALL_NETCDF(nc_inq_vardimid(ncFD, varId, dimIds));
584
585   // Number of values to read.
586   vtkIdType arraySize = 1;
587
588   // Indices to read from.
589   size_t start[4], count[4];
590
591   // Are we using time?
592   int timeIndexOffset = 0;
593   if ((numDims > 0) && this->IsTimeDimension(ncFD, dimIds[0]))
594     {
595     vtkSmartPointer<vtkDoubleArray> timeValues
596       = this->GetTimeValues(ncFD, dimIds[0]);
597     timeIndexOffset = 1;
598     // Find the index for the given time.
599     // We could speed this up with a binary search or something.
600     for (start[0] = 0;
601          start[0] < static_cast<size_t>(timeValues->GetNumberOfTuples());
602          start[0]++)
603       {
604       if (timeValues->GetValue(start[0]) >= time) break;
605       }
606     count[0] = 1;
607     numDims--;
608     }
609
610   if (numDims > 3)
611     {
612     vtkErrorMacro(<< "More than 3 dims without time not supported in variable "
613                   << varName);
614     return 0;
615     }
616
617   bool loadingPointData = this->DimensionsAreForPointData(
618                                   this->LoadingDimensions->GetPointer(0),
619                                   this->LoadingDimensions->GetNumberOfTuples());
620
621   // Set up read indices.  Also check to make sure the dimensions are consistent
622   // with other loaded variables.
623   int extent[6];
624   output->GetUpdateExtent(extent);
625   if (numDims != this->LoadingDimensions->GetNumberOfTuples())
626     {
627     vtkWarningMacro(<< "Variable " << varName << " dimensions ("
628                     << this->DescribeDimensions(ncFD, dimIds+timeIndexOffset,
629                                                 numDims).c_str()
630                     << ") are different than the other variable dimensions ("
631                     << this->DescribeDimensions(ncFD,
632                            this->LoadingDimensions->GetPointer(0),
633                            this->LoadingDimensions->GetNumberOfTuples()).c_str()
634                     << ").  Skipping");
635     return 1;
636     }
637   for (int i = 0; i < numDims; i++)
638     {
639     if (dimIds[i+timeIndexOffset] != this->LoadingDimensions->GetValue(i))
640       {
641       vtkWarningMacro(<< "Variable " << varName << " dimensions ("
642                       << this->DescribeDimensions(ncFD, dimIds+timeIndexOffset,
643                                                   numDims).c_str()
644                       << ") are different than the other variable dimensions ("
645                       << this->DescribeDimensions(ncFD,
646                            this->LoadingDimensions->GetPointer(0),
647                            this->LoadingDimensions->GetNumberOfTuples()).c_str()
648                       << ").  Skipping");
649       return 1;
650       }
651     // Remember that netCDF arrays are indexed backward from VTK images.
652     start[i+timeIndexOffset] = extent[2*(numDims-i-1)];
653     count[i+timeIndexOffset]
654       = extent[2*(numDims-i-1)+1]-extent[2*(numDims-i-1)]+1;
655
656     // If loading cell data, subtract one from the data being loaded.
657     if (!loadingPointData) count[i+timeIndexOffset]--;
658
659     arraySize *= count[i+timeIndexOffset];
660     }
661
662   // Allocate an array of the right type.
663   nc_type ncType;
664   CALL_NETCDF(nc_inq_vartype(ncFD, varId, &ncType));
665   int vtkType = NetCDFTypeToVTKType(ncType);
666   if (vtkType < 1) return 0;
667   vtkSmartPointer<vtkDataArray> dataArray;
668   dataArray.TakeReference(vtkDataArray::CreateDataArray(vtkType));
669   dataArray->SetNumberOfComponents(1);
670   dataArray->SetNumberOfTuples(arraySize);
671
672   // Read the array from the file.
673   CALL_NETCDF(nc_get_vars(ncFD, varId, start, count, NULL,
674                           dataArray->GetVoidPointer(0)));
675
676   // Check for a fill value.
677   size_t attribLength;
678   if (   (nc_inq_attlen(ncFD, varId, "_FillValue", &attribLength) == NC_NOERR)
679       && (attribLength == 1) )
680     {
681     if (this->ReplaceFillValueWithNan)
682       {
683       // NaN only available with float and double.
684       if (dataArray->GetDataType() == VTK_FLOAT)
685         {
686         float fillValue;
687         nc_get_att_float(ncFD, varId, "_FillValue", &fillValue);
688         vtkstd::replace(reinterpret_cast<float*>(dataArray->GetVoidPointer(0)),
689                         reinterpret_cast<float*>(dataArray->GetVoidPointer(
690                                                dataArray->GetNumberOfTuples())),
691                         fillValue, static_cast<float>(vtkMath::Nan()));
692         }
693       else if (dataArray->GetDataType() == VTK_DOUBLE)
694         {
695         double fillValue;
696         nc_get_att_double(ncFD, varId, "_FillValue", &fillValue);
697         vtkstd::replace(reinterpret_cast<double*>(dataArray->GetVoidPointer(0)),
698                         reinterpret_cast<double*>(dataArray->GetVoidPointer(
699                                                dataArray->GetNumberOfTuples())),
700                         fillValue, vtkMath::Nan());
701         }
702       else
703         {
704         vtkDebugMacro(<< "No NaN available for data of type "
705                       << dataArray->GetDataType());
706         }
707       }
708     }
709
710   // Check to see if there is a scale or offset.
711   double scale = 1.0;
712   double offset = 0.0;
713   if (   (nc_inq_attlen(ncFD, varId, "scale_factor", &attribLength) == NC_NOERR)
714       && (attribLength == 1) )
715     {
716     CALL_NETCDF(nc_get_att_double(ncFD, varId, "scale_factor", &scale));
717     }
718   if (   (nc_inq_attlen(ncFD, varId, "add_offset", &attribLength) == NC_NOERR)
719       && (attribLength == 1) )
720     {
721     CALL_NETCDF(nc_get_att_double(ncFD, varId, "add_offset", &offset));
722     }
723
724   if ((scale != 1.0) || (offset != 0.0))
725     {
726     VTK_CREATE(vtkDoubleArray, adjustedArray);
727     adjustedArray->SetNumberOfComponents(1);
728     adjustedArray->SetNumberOfTuples(arraySize);
729     for (vtkIdType i = 0; i < arraySize; i++)
730       {
731       adjustedArray->SetValue(i, dataArray->GetTuple1(i)*scale + offset);
732       }
733     dataArray = adjustedArray;
734     }
735
736   // Add data to the output.
737   dataArray->SetName(varName);
738   if (loadingPointData)
739     {
740     output->GetPointData()->AddArray(dataArray);
741     }
742   else
743     {
744     output->GetCellData()->AddArray(dataArray);
745     }
746
747   return 1;
748 }
749