1 // -*- c++ -*-
2 /*=========================================================================
3
4 Program: Visualization Toolkit
5 Module: vtkNetCDFCFReader.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 "vtkNetCDFCFReader.h"
24
25 #include "vtkCellArray.h"
26 #include "vtkDataArraySelection.h"
27 #include "vtkDoubleArray.h"
28 #include "vtkExtentTranslator.h"
29 #include "vtkImageData.h"
30 #include "vtkInformation.h"
31 #include "vtkInformationVector.h"
32 #include "vtkIntArray.h"
33 #include "vtkMath.h"
34 #include "vtkMergePoints.h"
35 #include "vtkNew.h"
36 #include "vtkObjectFactory.h"
37 #include "vtkPoints.h"
38 #include "vtkRectilinearGrid.h"
39 #include "vtkStdString.h"
40 #include "vtkStringArray.h"
41 #include "vtkStreamingDemandDrivenPipeline.h"
42 #include "vtkStructuredGrid.h"
43 #include "vtkUnstructuredGrid.h"
44
45 #include "vtkSmartPointer.h"
46 #define VTK_CREATE(type, name) \
47 vtkSmartPointer<type> name = vtkSmartPointer<type>::New()
48
49 #include <set>
50 #include <vtksys/RegularExpression.hxx>
51 #include <vtksys/SystemTools.hxx>
52
53 #include <cstring>
54
55 #include "vtk_netcdf.h"
56
57 #define CALL_NETCDF_GENERIC(call, on_error) \
58 { \
59 int errorcode = call; \
60 if (errorcode != NC_NOERR) \
61 { \
62 const char * errorstring = nc_strerror(errorcode); \
63 on_error; \
64 } \
65 }
66
67 #define CALL_NETCDF(call) \
68 CALL_NETCDF_GENERIC(call, vtkErrorMacro(<< "netCDF Error: " << errorstring); return 0;)
69
70 #define CALL_NETCDF_GW(call) \
71 CALL_NETCDF_GENERIC(call, vtkGenericWarningMacro(<< "netCDF Error: " << errorstring); return 0;)
72
73 #include <algorithm>
74
75 #include <cmath>
76
77 //=============================================================================
78 // Convenience function for getting the text attribute on a variable. Returns
79 // true if the attribute exists, false otherwise.
ReadTextAttribute(int ncFD,int varId,const char * name,vtkStdString & result)80 static bool ReadTextAttribute(int ncFD, int varId, const char *name,
81 vtkStdString &result)
82 {
83 size_t length;
84 if (nc_inq_attlen(ncFD, varId, name, &length) != NC_NOERR) { return false; }
85
86 result.resize(length);
87 if (length > 0)
88 {
89 if (nc_get_att_text(ncFD, varId, name, &result.at(0)) != NC_NOERR)
90 {
91 return false;
92 }
93 }
94 else
95 {
96 // If length == 0, then there really is nothing to read. Do nothing
97 }
98
99 // The line below seems weird, but it is here for a good reason. In general,
100 // text attributes are not null terminated, so you have to add your own (which
101 // the vtkStdString will do for us). However, sometimes a null terminating
102 // character is written in the attribute anyway. In a C string this is no big
103 // deal. But it means that the vtkStdString has a null character in it and it
104 // is technically different than its own C string. This line corrects that
105 // regardless of whether the null string was written we will get the right
106 // string.
107 result = result.c_str();
108
109 return true;
110 }
111
112 //-----------------------------------------------------------------------------
113 // Convenience function for getting the range of all values in all components of
114 // a vtkDoubleArray.
GetRangeOfAllComponents(vtkDoubleArray * array,double range[2])115 static void GetRangeOfAllComponents(vtkDoubleArray *array, double range[2])
116 {
117 range[0] = VTK_DOUBLE_MAX;
118 range[1] = VTK_DOUBLE_MIN;
119 for (int component = 0; component<array->GetNumberOfComponents(); component++)
120 {
121 double componentRange[2];
122 array->GetRange(componentRange, component);
123 if (componentRange[0] < range[0]) range[0] = componentRange[0];
124 if (componentRange[1] > range[1]) range[1] = componentRange[1];
125 }
126 }
127
128 //=============================================================================
vtkDimensionInfo(int ncFD,int id)129 vtkNetCDFCFReader::vtkDimensionInfo::vtkDimensionInfo(int ncFD, int id)
130 {
131 this->DimId = id;
132
133 this->Units = UNDEFINED_UNITS;
134 this->HasRegularSpacing = true;
135 this->Origin = 0.0;
136 this->Spacing = 1.0;
137
138 this->LoadMetaData(ncFD);
139 }
140
LoadMetaData(int ncFD)141 int vtkNetCDFCFReader::vtkDimensionInfo::LoadMetaData(int ncFD)
142 {
143 this->Units = UNDEFINED_UNITS;
144
145 char name[NC_MAX_NAME+1];
146 CALL_NETCDF_GW(nc_inq_dimname(ncFD, this->DimId, name));
147 this->Name = name;
148
149 size_t dimLen;
150 CALL_NETCDF_GW(nc_inq_dimlen(ncFD, this->DimId, &dimLen));
151 this->Coordinates = vtkSmartPointer<vtkDoubleArray>::New();
152 this->Coordinates->SetName((this->Name + "_Coordinates").c_str());
153 this->Coordinates->SetNumberOfComponents(1);
154 this->Coordinates->SetNumberOfTuples(static_cast<vtkIdType>(dimLen));
155
156 this->Bounds = vtkSmartPointer<vtkDoubleArray>::New();
157 this->Bounds->SetName((this->Name + "_Bounds").c_str());
158 this->Bounds->SetNumberOfComponents(1);
159 this->Bounds->SetNumberOfTuples(static_cast<vtkIdType>(dimLen+1));
160
161 this->SpecialVariables = vtkSmartPointer<vtkStringArray>::New();
162
163 int varId;
164 int varNumDims;
165 int varDim;
166 // By convention if there is a single dimension variable with the same name as
167 // its dimension, then the data contains the coordinates for the dimension.
168 if ( (nc_inq_varid(ncFD, name, &varId) == NC_NOERR)
169 && (nc_inq_varndims(ncFD, varId, &varNumDims) == NC_NOERR)
170 && (varNumDims == 1)
171 && (nc_inq_vardimid(ncFD, varId, &varDim) == NC_NOERR)
172 && (varDim == this->DimId) )
173 {
174 this->SpecialVariables->InsertNextValue(name);
175
176 // Read coordinates
177 CALL_NETCDF_GW(nc_get_var_double(ncFD, varId,
178 this->Coordinates->GetPointer(0)));
179
180 // Check to see if the spacing is regular.
181 this->Origin = this->Coordinates->GetValue(0);
182 this->Spacing
183 = (this->Coordinates->GetValue(static_cast<vtkIdType>(dimLen-1)) - this->Origin)/(dimLen-1);
184 this->HasRegularSpacing = true; // Then check to see if it is false.
185 double tolerance = 0.01*this->Spacing;
186 for (size_t i = 1; i < dimLen; i++)
187 {
188 double expectedValue = this->Origin + i*this->Spacing;
189 double actualValue = this->Coordinates->GetValue(static_cast<vtkIdType>(i));
190 if ( (actualValue < expectedValue-tolerance)
191 || (actualValue > expectedValue+tolerance) )
192 {
193 this->HasRegularSpacing = false;
194 break;
195 }
196 }
197
198 // Check units.
199 vtkStdString units;
200 if (ReadTextAttribute(ncFD, varId, "units", units))
201 {
202 units = vtksys::SystemTools::LowerCase(units);
203 // Time, latitude, and longitude dimensions are those with units that
204 // correspond to strings formatted with the Unidata udunits package. I'm
205 // not sure if these checks are complete, but they matches all of the
206 // examples I have seen.
207 if ((units.find(" since ") != vtkStdString::npos)
208 || (units.find(" after ") != vtkStdString::npos)
209 || (units == "second")
210 || (units == "seconds")
211 || (units == "day")
212 || (units == "days")
213 || (units == "hour")
214 || (units == "hours")
215 || (units == "minute")
216 || (units == "minutes")
217 || (units == "s")
218 || (units == "sec")
219 || (units == "secs")
220 || (units == "shake")
221 || (units == "shakes")
222 || (units == "sidereal_day")
223 || (units == "sidereal_days")
224 || (units == "sidereal_hour")
225 || (units == "sidereal_hours")
226 || (units == "sidereal_minute")
227 || (units == "sidereal_minutes")
228 || (units == "sidereal_second")
229 || (units == "sidereal_seconds")
230 || (units == "sidereal_year")
231 || (units == "sidereal_years")
232 || (units == "tropical_year")
233 || (units == "tropical_years")
234 || (units == "lunar_month")
235 || (units == "lunar_months")
236 || (units == "common_year")
237 || (units == "common_years")
238 || (units == "leap_year")
239 || (units == "leap_years")
240 || (units == "Julian_year")
241 || (units == "Julian_years")
242 || (units == "Gregorian_year")
243 || (units == "Gregorian_years")
244 || (units == "sidereal_month")
245 || (units == "sidereal_months")
246 || (units == "tropical_month")
247 || (units == "tropical_months")
248 || (units == "d")
249 || (units == "min")
250 || (units == "mins")
251 || (units == "hrs")
252 || (units == "h")
253 || (units == "fortnight")
254 || (units == "fortnights")
255 || (units == "week")
256 || (units == "jiffy")
257 || (units == "jiffies")
258 || (units == "year")
259 || (units == "years")
260 || (units == "yr")
261 || (units == "yrs")
262 || (units == "a")
263 || (units == "eon")
264 || (units == "eons")
265 || (units == "month")
266 || (units == "months") )
267 {
268 this->Units = TIME_UNITS;
269 }
270 else if (vtksys::RegularExpression("degrees?_?n").find(units))
271 {
272 this->Units = LATITUDE_UNITS;
273 }
274 else if (vtksys::RegularExpression("degrees?_?e").find(units))
275 {
276 this->Units = LONGITUDE_UNITS;
277 }
278 }
279
280 // Check axis.
281 vtkStdString axis;
282 if (ReadTextAttribute(ncFD, varId, "axis", axis))
283 {
284 // The axis attribute is an alternate way of defining the coordinate type.
285 // The string can be "X", "Y", "Z", or "T" which mean longitude, latitude,
286 // vertical, and time, respectively.
287 if (axis == "X")
288 {
289 this->Units = LONGITUDE_UNITS;
290 }
291 else if (axis == "Y")
292 {
293 this->Units = LATITUDE_UNITS;
294 }
295 else if (axis == "Z")
296 {
297 this->Units = VERTICAL_UNITS;
298 }
299 else if (axis == "T")
300 {
301 this->Units = TIME_UNITS;
302 }
303 }
304
305 // Check positive.
306 vtkStdString positive;
307 if (ReadTextAttribute(ncFD, varId, "positive", positive))
308 {
309 positive = vtksys::SystemTools::LowerCase(positive);
310 if (positive.find("down") != vtkStdString::npos)
311 {
312 // Flip the values of the coordinates.
313 for (vtkIdType i = 0; i < this->Coordinates->GetNumberOfTuples(); i++)
314 {
315 this->Coordinates->SetValue(i, -(this->Coordinates->GetValue(i)));
316 }
317 this->Spacing = -this->Spacing;
318 }
319 }
320
321 // Create the bounds array, which is used in place of the coordinates when
322 // loading as cell data. We will look for the bounds attribute on the
323 // description variable to see if cell bounds have been written out. This
324 // code assumes that if such an attribute exists, it is the name of another
325 // variable that is of dimensions of size dimLen X 2. There are no checks
326 // for this (other than the existence of the attribute), so if this is not
327 // the case then the code could fail.
328 vtkStdString boundsName;
329 if (ReadTextAttribute(ncFD, varId, "bounds", boundsName))
330 {
331 this->SpecialVariables->InsertNextValue(boundsName);
332
333 int boundsVarId;
334 CALL_NETCDF_GW(nc_inq_varid(ncFD, boundsName.c_str(), &boundsVarId));
335
336 // Read in the first bound value for each entry as a point bound. If the
337 // cells are connected, the second bound value should equal the first
338 // bound value of the next entry anyway.
339 size_t start[2]; start[0] = start[1] = 0;
340 size_t count[2]; count[0] = dimLen; count[1] = 1;
341 CALL_NETCDF_GW(nc_get_vars_double(ncFD, boundsVarId, start, count, nullptr,
342 this->Bounds->GetPointer(0)));
343
344 // Read in the last value for the bounds array. It will be the second
345 // bound in the last entry. This will not be replicated unless the
346 // dimension is a longitudinal one that wraps all the way around.
347 start[0] = dimLen-1; start[1] = 1;
348 count[0] = 1; count[1] = 1;
349 CALL_NETCDF_GW(nc_get_vars_double(ncFD, boundsVarId, start, count, nullptr,
350 this->Bounds->GetPointer(static_cast<vtkIdType>(dimLen))));
351 }
352 else
353 {
354 // Bounds not given. Set them based on the coordinates.
355 this->Bounds->SetValue(
356 0, this->Coordinates->GetValue(0) - 0.5*this->Spacing);
357 for (vtkIdType i = 1; i < static_cast<vtkIdType>(dimLen); i++)
358 {
359 double v0 = this->Coordinates->GetValue(i-1);
360 double v1 = this->Coordinates->GetValue(i);
361 this->Bounds->SetValue(i, 0.5*(v0+v1));
362 }
363 this->Bounds->SetValue(static_cast<vtkIdType>(dimLen),
364 this->Coordinates->GetValue(static_cast<vtkIdType>(dimLen-1))+0.5*this->Spacing);
365 }
366 }
367 else
368 {
369 // Fake coordinates
370 for (size_t i = 0; i < dimLen; i++)
371 {
372 this->Coordinates->SetValue(static_cast<vtkIdType>(i),
373 static_cast<double>(i));
374 this->Bounds->SetValue(static_cast<vtkIdType>(i),
375 static_cast<double>(i) - 0.5);
376 }
377 this->Bounds->SetValue(static_cast<vtkIdType>(dimLen),
378 static_cast<double>(dimLen) - 0.5);
379 this->HasRegularSpacing = true;
380 this->Origin = 0.0;
381 this->Spacing = 1.0;
382 }
383
384 return 1;
385 }
386
387 //-----------------------------------------------------------------------------
388 class vtkNetCDFCFReader::vtkDimensionInfoVector
389 {
390 public:
391 std::vector<vtkDimensionInfo> v;
392 };
393
394 //=============================================================================
vtkDependentDimensionInfo(int ncFD,int varId,vtkNetCDFCFReader * parent)395 vtkNetCDFCFReader::vtkDependentDimensionInfo::vtkDependentDimensionInfo(
396 int ncFD, int varId,
397 vtkNetCDFCFReader *parent)
398 {
399 if (this->LoadMetaData(ncFD, varId, parent))
400 {
401 this->Valid = true;
402 }
403 else
404 {
405 this->Valid = false;
406 }
407 }
408
409 //-----------------------------------------------------------------------------
LoadMetaData(int ncFD,int varId,vtkNetCDFCFReader * parent)410 int vtkNetCDFCFReader::vtkDependentDimensionInfo::LoadMetaData(
411 int ncFD, int varId,
412 vtkNetCDFCFReader *parent)
413 {
414 int longitudeCoordVarId, latitudeCoordVarId;
415 int longitudeBoundsVarId, latitudeBoundsVarId;
416 longitudeCoordVarId = latitudeCoordVarId = -1;
417 longitudeBoundsVarId = latitudeBoundsVarId = -1;
418 this->GridDimensions = vtkSmartPointer<vtkIntArray>::New();
419 this->SpecialVariables = vtkSmartPointer<vtkStringArray>::New();
420
421 // The grid dimensions are the dimensions on the variable. Since multiple
422 // variables can be put on the same grid and this class identifies grids by
423 // their variables, I group all of the dimension combinations together for 2D
424 // coordinate lookup. Technically, the CF specification allows you to specify
425 // different coordinate variables, but we do not support that because there is
426 // no easy way to differentiate grids with the same dimensions. If different
427 // coordinates are needed, then duplicate dimensions should be created.
428 // Anyone who disagrees should write their own reader class.
429 int numGridDimensions;
430 CALL_NETCDF_GW(nc_inq_varndims(ncFD, varId, &numGridDimensions));
431
432 if(numGridDimensions == 0)
433 {
434 // If a variable has no dimensions, there is no reason to have dependent
435 // dimension variables. Just exit here for safety.
436 return 0;
437 }
438
439 this->GridDimensions->SetNumberOfTuples(numGridDimensions);
440 CALL_NETCDF_GW(nc_inq_vardimid(ncFD, varId,
441 this->GridDimensions->GetPointer(0)));
442
443 // Remove initial time dimension, which has no effect on data type.
444 if (parent->IsTimeDimension(ncFD, this->GridDimensions->GetValue(0)))
445 {
446 this->GridDimensions->RemoveTuple(0);
447 numGridDimensions--;
448 if(numGridDimensions == 0)
449 {
450 // If a variable has no dimensions (ignoring time), there is no reason
451 // to have dependent dimension variables. Just exit here for safety.
452 return 0;
453 }
454 }
455
456 // Most coordinate variables are defined by a variable the same name as the
457 // dimension they describe. Those are handled elsewhere. This class handles
458 // dependent variables that define coordinates that are not the same name as
459 // any dimension. This is only done when the coordinates cannot be expressed
460 // as a 1D table lookup from dimension index. This occurs in only two places
461 // in the CF convention. First, it happens for 2D coordinate variables with
462 // 4-sided cells. This is basically when the grid is a 2D curvilinear grid.
463 // Each i,j topological point can be placed anywhere in space. Second, it
464 // happens for multi-dimensional coordinate variables with p-sided cells.
465 // These are unstructured collections of polygons.
466
467 vtkStdString coordinates;
468 if (!ReadTextAttribute(ncFD, varId, "coordinates", coordinates)) return 0;
469
470 std::vector<std::string> coordName;
471 vtksys::SystemTools::Split(coordinates, coordName, ' ');
472
473 int numAuxCoordDims = -1;
474
475 for (std::vector<std::string>::iterator iter = coordName.begin();
476 iter != coordName.end(); ++iter)
477 {
478 int auxCoordVarId;
479 if (nc_inq_varid(ncFD, iter->c_str(), &auxCoordVarId) != NC_NOERR) continue;
480
481 // Make sure that the coordinate variables have the same dimensions and that
482 // those dimensions are the same as the last two dimensions on the grid.
483 // Not sure if that is enforced by the specification, but I am going to make
484 // that assumption.
485 int numDims;
486 CALL_NETCDF_GW(nc_inq_varndims(ncFD, auxCoordVarId, &numDims));
487 // I am only supporting either 1 or 2 dimensions in the coordinate
488 // variables. See the comment below regarding identifying the
489 // CellsUnstructured flag.
490 if (numDims > 2) continue;
491
492 int auxCoordDims[2];
493 CALL_NETCDF_GW(nc_inq_vardimid(ncFD, auxCoordVarId, auxCoordDims));
494 int *gridDims = this->GridDimensions->GetPointer(numGridDimensions-numDims);
495 bool auxCoordDimsValid = true;
496 for (int dimId = 0; dimId < numDims; dimId++)
497 {
498 if (auxCoordDims[dimId] != gridDims[dimId])
499 {
500 auxCoordDimsValid = false;
501 break;
502 }
503 }
504 if (!auxCoordDimsValid) continue;
505
506 // The variable is no use to me unless it is identified as either longitude
507 // or latitude.
508 vtkStdString units;
509 if (!ReadTextAttribute(ncFD, auxCoordVarId, "units", units)) continue;
510 units = vtksys::SystemTools::LowerCase(units);
511 if (vtksys::RegularExpression("degrees?_?n").find(units))
512 {
513 latitudeCoordVarId = auxCoordVarId;
514 }
515 else if (vtksys::RegularExpression("degrees?_?e").find(units))
516 {
517 longitudeCoordVarId = auxCoordVarId;
518 }
519 else
520 {
521 continue;
522 }
523 this->SpecialVariables->InsertNextValue(*iter);
524
525 if (numAuxCoordDims < 0)
526 {
527 // Record the number of dimensions in the coordinate arrays.
528 numAuxCoordDims = numDims;
529 }
530 else if (numAuxCoordDims != numDims)
531 {
532 // Number of dimensions in different coordinate arrays do not match.
533 return 0;
534 }
535 }
536
537 if ((longitudeCoordVarId == -1) || (latitudeCoordVarId == -1))
538 {
539 // Did not find all coordinate variables.
540 return 0;
541 }
542
543 // Technically, p-sided cells can be accessed with any number of dimensions.
544 // However, it makes little sense to have more than a 1D array of cell ids.
545 // Supporting more than this is a pain because it exceeds the dimensionality
546 // supported by vtkDataArray. It also makes it exceedingly difficult to
547 // determine whether the topology is implicitly defined by 2D 4-sided cells or
548 // explicitly with p-sided cells. Thus, I am only implementing p-sided cells
549 // with 1D coordinate variables.
550 if (numAuxCoordDims == 1)
551 {
552 this->CellsUnstructured = true;
553 }
554 else if (numAuxCoordDims == 2)
555 {
556 this->CellsUnstructured = false;
557 }
558 else
559 {
560 // Not supporting this.
561 return 0;
562 }
563
564 vtkStdString bounds;
565 if (ReadTextAttribute(ncFD, longitudeCoordVarId, "bounds", bounds))
566 {
567 // The bounds is supposed to point to an array with numAuxCoordDims+1
568 // dimensions. The first numAuxCoordDims should be the same as the coord
569 // arrays. The last dimension has the number of vertices in each cell.
570 // Maybe I should check this, but I'm not.
571 CALL_NETCDF_GW(nc_inq_varid(ncFD, bounds.c_str(), &longitudeBoundsVarId));
572 this->SpecialVariables->InsertNextValue(bounds);
573 }
574 if (ReadTextAttribute(ncFD, latitudeCoordVarId, "bounds", bounds))
575 {
576 // The bounds is supposed to point to an array with numAuxCoordDims+1
577 // dimensions. The first numAuxCoordDims should be the same as the coord
578 // arrays. The last dimension has the number of vertices in each cell.
579 // Maybe I should check this, but I'm not.
580 CALL_NETCDF_GW(nc_inq_varid(ncFD, bounds.c_str(), &latitudeBoundsVarId));
581 this->SpecialVariables->InsertNextValue(bounds);
582 }
583
584 this->HasBounds = ((longitudeBoundsVarId != -1)&&(latitudeBoundsVarId != -1));
585
586 // Load in all the longitude and latitude coordinates. Maybe not the most
587 // efficient thing to do for large data, but it is just a 2D surface, so it
588 // should be OK for most things.
589 this->LongitudeCoordinates = vtkSmartPointer<vtkDoubleArray>::New();
590 this->LatitudeCoordinates = vtkSmartPointer<vtkDoubleArray>::New();
591 if (this->CellsUnstructured)
592 {
593 if (this->HasBounds)
594 {
595 int ok;
596 ok = this->LoadUnstructuredBoundsVariable(ncFD, longitudeBoundsVarId,
597 this->LongitudeCoordinates);
598 if (!ok) return 0;
599 ok = this->LoadUnstructuredBoundsVariable(ncFD, latitudeBoundsVarId,
600 this->LatitudeCoordinates);
601 if (!ok) return 0;
602 }
603 else
604 {
605 // Unstructured cells need bounds to define the topology.
606 return 0;
607 }
608 }
609 else
610 {
611 if (this->HasBounds)
612 {
613 int ok;
614 ok = this->LoadBoundsVariable(ncFD, longitudeBoundsVarId,
615 this->LongitudeCoordinates);
616 if (!ok) return 0;
617 ok = this->LoadBoundsVariable(ncFD, latitudeBoundsVarId,
618 this->LatitudeCoordinates);
619 if (!ok) return 0;
620 }
621 else
622 {
623 int ok;
624 ok = this->LoadCoordinateVariable(ncFD, longitudeCoordVarId,
625 this->LongitudeCoordinates);
626 if (!ok) return 0;
627 ok = this->LoadCoordinateVariable(ncFD, latitudeCoordVarId,
628 this->LatitudeCoordinates);
629 if (!ok) return 0;
630 }
631 }
632
633 return 1;
634 }
635
636 //-----------------------------------------------------------------------------
LoadCoordinateVariable(int ncFD,int varId,vtkDoubleArray * coords)637 int vtkNetCDFCFReader::vtkDependentDimensionInfo::LoadCoordinateVariable(
638 int ncFD, int varId,
639 vtkDoubleArray *coords)
640 {
641 int dimIds[2];
642 CALL_NETCDF_GW(nc_inq_vardimid(ncFD, varId, dimIds));
643
644 size_t dimSizes[2];
645 for (int i = 0; i < 2; i++)
646 {
647 CALL_NETCDF_GW(nc_inq_dimlen(ncFD, dimIds[i], &dimSizes[i]));
648 }
649
650 coords->SetNumberOfComponents(static_cast<int>(dimSizes[1]));
651 coords->SetNumberOfTuples(static_cast<vtkIdType>(dimSizes[0]));
652 CALL_NETCDF_GW(nc_get_var_double(ncFD, varId, coords->GetPointer(0)));
653
654 return 1;
655 }
656
657 //-----------------------------------------------------------------------------
LoadBoundsVariable(int ncFD,int varId,vtkDoubleArray * coords)658 int vtkNetCDFCFReader::vtkDependentDimensionInfo::LoadBoundsVariable(
659 int ncFD, int varId,
660 vtkDoubleArray *coords)
661 {
662 int dimIds[3];
663 CALL_NETCDF_GW(nc_inq_vardimid(ncFD, varId, dimIds));
664
665 size_t dimSizes[3];
666 for (int i = 0; i < 3; i++)
667 {
668 CALL_NETCDF_GW(nc_inq_dimlen(ncFD, dimIds[i], &dimSizes[i]));
669 }
670
671 if (dimSizes[2] != 4)
672 {
673 vtkGenericWarningMacro(<< "Expected 2D dependent coordinate bounds to have"
674 << " 4 entries in final dimension. Instead has "
675 << dimSizes[2]);
676 return 0;
677 }
678
679 // Bounds are stored as 4-tuples for every cell. Tuple entries 0 and 1
680 // connect to the cell in the -i topological direction. Tuple entries 0 and 3
681 // connect to the cell in the -j topological direction.
682 std::vector<double> boundsData(dimSizes[0]*dimSizes[1]*4);
683 if (!boundsData.empty())
684 {
685 CALL_NETCDF_GW(nc_get_var_double(ncFD, varId, &boundsData.at(0)));
686 }
687
688 // The coords array are the coords at the points. There is one more point
689 // than cell in each topological direction.
690 int numComponents = static_cast<int>(dimSizes[1]);
691 vtkIdType numTuples = static_cast<vtkIdType>(dimSizes[0]);
692 coords->SetNumberOfComponents(numComponents+1);
693 coords->SetNumberOfTuples(numTuples+1);
694
695 // Copy from the bounds data to the coordinates data. Most values will
696 // be copied from the bound's 0'th tuple entry. Values at the extremes
697 // will be copied from other entries.
698 for (vtkIdType j = 0; j < numTuples; j++)
699 {
700 for (int i = 0; i < numComponents; i++)
701 {
702 coords->SetComponent(j, i, boundsData[(j*numComponents + i)*4 + 0]);
703 }
704 coords->SetComponent(j, numComponents,
705 boundsData[((j+1)*numComponents-1)*4 + 1]);
706 }
707 for (int i = 0; i < numComponents; i++)
708 {
709 coords->SetComponent(numTuples, i,
710 boundsData[((numTuples-1)*numComponents)*4 + 2]);
711 }
712 coords->SetComponent(numTuples, numComponents,
713 boundsData[(numTuples*numComponents-1)*4 + 3]);
714
715 return 1;
716 }
717
718 //-----------------------------------------------------------------------------
LoadUnstructuredBoundsVariable(int ncFD,int varId,vtkDoubleArray * coords)719 int vtkNetCDFCFReader::vtkDependentDimensionInfo::LoadUnstructuredBoundsVariable(
720 int ncFD, int varId, vtkDoubleArray *coords)
721 {
722 int dimIds[2];
723 CALL_NETCDF_GW(nc_inq_vardimid(ncFD, varId, dimIds));
724
725 size_t dimSizes[2];
726 for (int i = 0; i < 2; i++)
727 {
728 CALL_NETCDF_GW(nc_inq_dimlen(ncFD, dimIds[i], &dimSizes[i]));
729 }
730
731 int numVertPerCell = static_cast<int>(dimSizes[1]);
732 vtkIdType numCells = static_cast<vtkIdType>(dimSizes[0]);
733
734 coords->SetNumberOfComponents(numVertPerCell);
735 coords->SetNumberOfTuples(numCells);
736 CALL_NETCDF_GW(nc_get_var_double(ncFD, varId, coords->GetPointer(0)));
737
738 return 1;
739 }
740
741 //-----------------------------------------------------------------------------
742 class vtkNetCDFCFReader::vtkDependentDimensionInfoVector
743 {
744 public:
745 std::vector<vtkNetCDFCFReader::vtkDependentDimensionInfo> v;
746 };
747
748 //=============================================================================
749 vtkStandardNewMacro(vtkNetCDFCFReader);
750
751 //-----------------------------------------------------------------------------
vtkNetCDFCFReader()752 vtkNetCDFCFReader::vtkNetCDFCFReader()
753 {
754 this->SphericalCoordinates = 1;
755 this->VerticalScale = 1.0;
756 this->VerticalBias = 0.0;
757 this->OutputType = -1;
758
759 this->DimensionInfo = new vtkDimensionInfoVector;
760 this->DependentDimensionInfo = new vtkDependentDimensionInfoVector;
761 }
762
~vtkNetCDFCFReader()763 vtkNetCDFCFReader::~vtkNetCDFCFReader()
764 {
765 delete this->DimensionInfo;
766 delete this->DependentDimensionInfo;
767 }
768
PrintSelf(ostream & os,vtkIndent indent)769 void vtkNetCDFCFReader::PrintSelf(ostream &os, vtkIndent indent)
770 {
771 this->Superclass::PrintSelf(os, indent);
772
773 os << indent << "SphericalCoordinates: " << this->SphericalCoordinates <<endl;
774 os << indent << "VerticalScale: " << this->VerticalScale <<endl;
775 os << indent << "VerticalBias: " << this->VerticalBias <<endl;
776 os << indent << "OutputType: " << this->OutputType <<endl;
777 }
778
779 //-----------------------------------------------------------------------------
CanReadFile(const char * filename)780 int vtkNetCDFCFReader::CanReadFile(const char *filename)
781 {
782 // We really just read basic arrays from netCDF files. If the netCDF library
783 // says we can read it, then we can read it.
784 int ncFD;
785 int errorcode = nc_open(filename, NC_NOWRITE, &ncFD);
786 if (errorcode == NC_NOERR)
787 {
788 nc_close(ncFD);
789 return 1;
790 }
791 else
792 {
793 return 0;
794 }
795 }
796
797 //-----------------------------------------------------------------------------
SetOutputType(int type)798 void vtkNetCDFCFReader::SetOutputType(int type)
799 {
800 vtkDebugMacro(<< this->GetClassName() << " (" << this << "):"
801 " setting OutputType to " << type);
802 if (this->OutputType != type)
803 {
804 bool typeValid = ( (type == -1)
805 || (type == VTK_IMAGE_DATA)
806 || (type == VTK_RECTILINEAR_GRID)
807 || (type == VTK_STRUCTURED_GRID)
808 || (type == VTK_UNSTRUCTURED_GRID) );
809 if (typeValid)
810 {
811 this->OutputType = type;
812 this->Modified();
813 }
814 else
815 {
816 vtkErrorMacro(<< "Invalid OutputType: " << type);
817 }
818 }
819 }
820
821 //-----------------------------------------------------------------------------
RequestDataObject(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)822 int vtkNetCDFCFReader::RequestDataObject(
823 vtkInformation *vtkNotUsed(request),
824 vtkInformationVector **vtkNotUsed(inputVector),
825 vtkInformationVector *outputVector)
826 {
827 vtkInformation *outInfo = outputVector->GetInformationObject(0);
828 vtkDataObject *output = vtkDataObject::GetData(outInfo);
829
830 // This is really too early to know the appropriate data type. We need to
831 // have meta data and let the user select arrays. We have to do part
832 // of the RequestInformation to get the appropriate meta data.
833 if (!this->UpdateMetaData()) return 0;
834
835 // Check that the dataType is correct or automatically set it if it is set to
836 // -1.
837 int dataType = this->OutputType;
838
839 int ncFD;
840 CALL_NETCDF(nc_open(this->FileName, NC_NOWRITE, &ncFD));
841
842 int numArrays = this->VariableArraySelection->GetNumberOfArrays();
843 for (int arrayIndex = 0; arrayIndex < numArrays; arrayIndex++)
844 {
845 if (!this->VariableArraySelection->GetArraySetting(arrayIndex)) continue;
846
847 const char *name = this->VariableArraySelection->GetArrayName(arrayIndex);
848 int varId;
849 CALL_NETCDF(nc_inq_varid(ncFD, name, &varId));
850
851 int currentNumDims;
852 CALL_NETCDF(nc_inq_varndims(ncFD, varId, ¤tNumDims));
853 if (currentNumDims < 1) continue;
854 VTK_CREATE(vtkIntArray, currentDimensions);
855 currentDimensions->SetNumberOfComponents(1);
856 currentDimensions->SetNumberOfTuples(currentNumDims);
857 CALL_NETCDF(nc_inq_vardimid(ncFD, varId,
858 currentDimensions->GetPointer(0)));
859
860 // Remove initial time dimension, which has no effect on data type.
861 if (this->IsTimeDimension(ncFD, currentDimensions->GetValue(0)))
862 {
863 currentDimensions->RemoveTuple(0);
864 currentNumDims--;
865 if (currentNumDims < 1) continue;
866 }
867
868 CoordinateTypesEnum coordType = this->CoordinateType(currentDimensions);
869
870 int preferredDataType;
871 switch(coordType)
872 {
873 case COORDS_UNIFORM_RECTILINEAR:
874 preferredDataType = VTK_IMAGE_DATA;
875 break;
876 case COORDS_NONUNIFORM_RECTILINEAR:
877 preferredDataType = VTK_RECTILINEAR_GRID;
878 break;
879 case COORDS_REGULAR_SPHERICAL:
880 case COORDS_2D_EUCLIDEAN:
881 case COORDS_2D_SPHERICAL:
882 case COORDS_EUCLIDEAN_4SIDED_CELLS:
883 case COORDS_SPHERICAL_4SIDED_CELLS:
884 preferredDataType = VTK_STRUCTURED_GRID;
885 break;
886 case COORDS_EUCLIDEAN_PSIDED_CELLS:
887 case COORDS_SPHERICAL_PSIDED_CELLS:
888 preferredDataType = VTK_UNSTRUCTURED_GRID;
889 break;
890 default:
891 vtkErrorMacro(<< "Internal error: unknown coordinate type.");
892 return 0;
893 }
894
895 // Check the data type.
896 if (dataType == -1)
897 {
898 dataType = preferredDataType;
899 }
900 else
901 {
902 switch (dataType)
903 {
904 case VTK_IMAGE_DATA:
905 if (preferredDataType != VTK_IMAGE_DATA)
906 {
907 vtkWarningMacro("You have set the OutputType to a data type that"
908 " cannot fully represent the topology of the data."
909 " Some of the topology will be ignored.");
910 }
911 break;
912 case VTK_RECTILINEAR_GRID:
913 if ( (preferredDataType != VTK_IMAGE_DATA)
914 && (preferredDataType != VTK_RECTILINEAR_GRID) )
915 {
916 vtkWarningMacro("You have set the OutputType to a data type that"
917 " cannot fully represent the topology of the data."
918 " Some of the topology will be ignored.");
919 }
920 break;
921 case VTK_STRUCTURED_GRID:
922 if ( (preferredDataType != VTK_IMAGE_DATA)
923 && (preferredDataType != VTK_RECTILINEAR_GRID)
924 && (preferredDataType != VTK_STRUCTURED_GRID) )
925 {
926 vtkWarningMacro("You have set the OutputType to a data type that"
927 " cannot fully represent the topology of the data."
928 " Some of the topology will be ignored.");
929 }
930 break;
931 case VTK_UNSTRUCTURED_GRID:
932 // Unstructured grid supports all topologies.
933 break;
934 default:
935 vtkErrorMacro(<< "Sanity check failed: bad internal type.");
936 return 0;
937 }
938 }
939
940 // That's right, break. We really only want to look at the dimensions of
941 // the first valid loaded variable. If we got here, we found a variable.
942 // The loop is really only for the continues at the top.
943 break;
944 }
945
946 CALL_NETCDF(nc_close(ncFD));
947
948 if (dataType == -1)
949 {
950 // Either no variables are selected or only variables with no dimensions.
951 // In either case, an image data should be fine.
952 dataType = VTK_IMAGE_DATA;
953 }
954
955 if (dataType == VTK_IMAGE_DATA)
956 {
957 if (!output || !output->IsA("vtkImageData"))
958 {
959 output = vtkImageData::New();
960 outInfo->Set(vtkDataObject::DATA_OBJECT(), output);
961 output->Delete(); // Not really deleted.
962 }
963 }
964 else if (dataType == VTK_RECTILINEAR_GRID)
965 {
966 if (!output || !output->IsA("vtkRectilinearGrid"))
967 {
968 output = vtkRectilinearGrid::New();
969 outInfo->Set(vtkDataObject::DATA_OBJECT(), output);
970 output->Delete(); // Not really deleted.
971 }
972 }
973 else if (dataType == VTK_STRUCTURED_GRID)
974 {
975 if (!output || !output->IsA("vtkStructuredGrid"))
976 {
977 output = vtkStructuredGrid::New();
978 outInfo->Set(vtkDataObject::DATA_OBJECT(), output);
979 output->Delete(); // Not really deleted.
980 }
981 }
982 else if (dataType == VTK_UNSTRUCTURED_GRID)
983 {
984 if (!output || !output->IsA("vtkUnstructuredGrid"))
985 {
986 output = vtkUnstructuredGrid::New();
987 outInfo->Set(vtkDataObject::DATA_OBJECT(), output);
988 output->Delete(); // Not really deleted.
989 }
990 }
991 else
992 {
993 vtkErrorMacro(<< "Sanity check failed: bad internal type.");
994 return 0;
995 }
996
997 return 1;
998 }
999
1000 //-----------------------------------------------------------------------------
RequestInformation(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)1001 int vtkNetCDFCFReader::RequestInformation(vtkInformation *request,
1002 vtkInformationVector **inputVector,
1003 vtkInformationVector *outputVector)
1004 {
1005 // Let the superclass do the heavy lifting.
1006 if (!this->Superclass::RequestInformation(request, inputVector, outputVector))
1007 {
1008 return 0;
1009 }
1010
1011 // Superclass understands structured data, but we have to handle unstructured
1012 // "extents" (pieces).
1013 vtkInformation *outInfo = outputVector->GetInformationObject(0);
1014 vtkDataObject *output = vtkDataObject::GetData(outInfo);
1015 if (output)
1016 {
1017 if (output->GetExtentType() != VTK_3D_EXTENT)
1018 {
1019 outInfo->Set(
1020 CAN_HANDLE_PIECE_REQUEST(), 1);
1021 }
1022 else
1023 {
1024 outInfo->Set(CAN_PRODUCE_SUB_EXTENT(), 1);
1025 }
1026 }
1027 else
1028 {
1029 return 0;
1030 }
1031 return 1;
1032 }
1033
1034 //-----------------------------------------------------------------------------
RequestData(vtkInformation * request,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)1035 int vtkNetCDFCFReader::RequestData(vtkInformation *request,
1036 vtkInformationVector **inputVector,
1037 vtkInformationVector *outputVector)
1038 {
1039 // If the output does not directly support 3D extents, then we have to make
1040 // some from the piece information so the superclass knows what portion of
1041 // arrays to load.
1042 vtkDataObject *output = vtkDataObject::GetData(outputVector);
1043 if (output)
1044 {
1045 if (output->GetExtentType() == VTK_3D_EXTENT)
1046 {
1047 // Do nothing. 3D extents already set.
1048 }
1049 else if (output->GetExtentType() == VTK_PIECES_EXTENT)
1050 {
1051 int pieceNumber, numberOfPieces, ghostLevels;
1052 vtkInformation *outInfo = outputVector->GetInformationObject(0);
1053 pieceNumber = outInfo->Get(
1054 vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
1055 numberOfPieces = outInfo->Get(
1056 vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
1057 ghostLevels = outInfo->Get(
1058 vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
1059
1060 int extent[6];
1061 this->ExtentForDimensionsAndPiece(pieceNumber,
1062 numberOfPieces,
1063 ghostLevels,
1064 extent);
1065
1066 // Store the update extent in the output's information object to make it
1067 // easy to find whenever loading data for this object.
1068 output->GetInformation()->Set(
1069 vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent, 6);
1070 }
1071 else
1072 {
1073 vtkWarningMacro(<< "Invalid extent type encountered. Data arrays may"
1074 << " be loaded incorrectly.");
1075 }
1076 }
1077 else // output == nullptr
1078 {
1079 vtkErrorMacro(<< "No output object.");
1080 return 0;
1081 }
1082
1083 // Let the superclass do the heavy lifting.
1084 if (!this->Superclass::RequestData(request, inputVector, outputVector))
1085 {
1086 return 0;
1087 }
1088
1089 // Add spacing information defined by the COARDS conventions.
1090
1091 vtkImageData *imageOutput = vtkImageData::GetData(outputVector);
1092 if (imageOutput)
1093 {
1094 this->AddRectilinearCoordinates(imageOutput);
1095 }
1096
1097 vtkRectilinearGrid *rectilinearOutput
1098 = vtkRectilinearGrid::GetData(outputVector);
1099 if (rectilinearOutput)
1100 {
1101 switch (this->CoordinateType(this->LoadingDimensions))
1102 {
1103 case COORDS_EUCLIDEAN_PSIDED_CELLS:
1104 case COORDS_SPHERICAL_PSIDED_CELLS:
1105 // There is no sensible way to store p-sided cells in a structured grid.
1106 // Just fake some coordinates (related to ParaView bug #11543).
1107 this->FakeRectilinearCoordinates(rectilinearOutput);
1108 break;
1109 case COORDS_UNIFORM_RECTILINEAR:
1110 case COORDS_NONUNIFORM_RECTILINEAR:
1111 case COORDS_REGULAR_SPHERICAL:
1112 case COORDS_2D_EUCLIDEAN:
1113 case COORDS_2D_SPHERICAL:
1114 case COORDS_EUCLIDEAN_4SIDED_CELLS:
1115 case COORDS_SPHERICAL_4SIDED_CELLS:
1116 default:
1117 this->AddRectilinearCoordinates(rectilinearOutput);
1118 }
1119 }
1120
1121 vtkStructuredGrid *structuredOutput
1122 = vtkStructuredGrid::GetData(outputVector);
1123 if (structuredOutput)
1124 {
1125 switch (this->CoordinateType(this->LoadingDimensions))
1126 {
1127 case COORDS_UNIFORM_RECTILINEAR:
1128 case COORDS_NONUNIFORM_RECTILINEAR:
1129 this->Add1DRectilinearCoordinates(structuredOutput);
1130 break;
1131 case COORDS_REGULAR_SPHERICAL:
1132 this->Add1DSphericalCoordinates(structuredOutput);
1133 break;
1134 case COORDS_2D_EUCLIDEAN:
1135 case COORDS_EUCLIDEAN_4SIDED_CELLS:
1136 this->Add2DRectilinearCoordinates(structuredOutput);
1137 break;
1138 case COORDS_2D_SPHERICAL:
1139 case COORDS_SPHERICAL_4SIDED_CELLS:
1140 this->Add2DSphericalCoordinates(structuredOutput);
1141 break;
1142 case COORDS_EUCLIDEAN_PSIDED_CELLS:
1143 case COORDS_SPHERICAL_PSIDED_CELLS:
1144 // There is no sensible way to store p-sided cells in a structured grid.
1145 // Just fake some coordinates (ParaView bug #11543).
1146 this->FakeStructuredCoordinates(structuredOutput);
1147 break;
1148 default:
1149 vtkErrorMacro("Internal error: unknown coordinate type.");
1150 return 0;
1151 }
1152 }
1153
1154 vtkUnstructuredGrid *unstructuredOutput
1155 = vtkUnstructuredGrid::GetData(outputVector);
1156 if (unstructuredOutput)
1157 {
1158 int extent[6];
1159 this->GetUpdateExtentForOutput(unstructuredOutput, extent);
1160
1161 switch (this->CoordinateType(this->LoadingDimensions))
1162 {
1163 case COORDS_UNIFORM_RECTILINEAR:
1164 case COORDS_NONUNIFORM_RECTILINEAR:
1165 this->Add1DRectilinearCoordinates(unstructuredOutput, extent);
1166 break;
1167 case COORDS_REGULAR_SPHERICAL:
1168 this->Add1DSphericalCoordinates(unstructuredOutput, extent);
1169 break;
1170 case COORDS_2D_EUCLIDEAN:
1171 case COORDS_EUCLIDEAN_4SIDED_CELLS:
1172 this->Add2DRectilinearCoordinates(unstructuredOutput, extent);
1173 break;
1174 case COORDS_2D_SPHERICAL:
1175 case COORDS_SPHERICAL_4SIDED_CELLS:
1176 this->Add2DSphericalCoordinates(unstructuredOutput, extent);
1177 break;
1178 case COORDS_EUCLIDEAN_PSIDED_CELLS:
1179 this->AddUnstructuredRectilinearCoordinates(unstructuredOutput, extent);
1180 break;
1181 case COORDS_SPHERICAL_PSIDED_CELLS:
1182 this->AddUnstructuredSphericalCoordinates(unstructuredOutput, extent);
1183 break;
1184 default:
1185 vtkErrorMacro("Internal error: unknown coordinate type.");
1186 return 0;
1187 }
1188 }
1189
1190 return 1;
1191 }
1192
1193 //-----------------------------------------------------------------------------
ExtentForDimensionsAndPiece(int pieceNumber,int numberOfPieces,int ghostLevels,int extent[6])1194 void vtkNetCDFCFReader::ExtentForDimensionsAndPiece(int pieceNumber,
1195 int numberOfPieces,
1196 int ghostLevels,
1197 int extent[6])
1198 {
1199 VTK_CREATE(vtkExtentTranslator, extentTranslator);
1200
1201 extentTranslator->SetWholeExtent(this->WholeExtent);
1202 extentTranslator->SetPiece(pieceNumber);
1203 extentTranslator->SetNumberOfPieces(numberOfPieces);
1204 extentTranslator->SetGhostLevel(ghostLevels);
1205
1206 extentTranslator->PieceToExtent();
1207
1208 extentTranslator->GetExtent(extent);
1209 }
1210
1211 //-----------------------------------------------------------------------------
GetUpdateExtentForOutput(vtkDataSet * output,int extent[6])1212 void vtkNetCDFCFReader::GetUpdateExtentForOutput(vtkDataSet *output,
1213 int extent[6])
1214 {
1215 vtkInformation *info = output->GetInformation();
1216 if (info->Has(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT()))
1217 {
1218 info->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent);
1219 }
1220 else
1221 {
1222 return this->Superclass::GetUpdateExtentForOutput(output, extent);
1223 }
1224 }
1225
1226 //-----------------------------------------------------------------------------
AddRectilinearCoordinates(vtkImageData * imageOutput)1227 void vtkNetCDFCFReader::AddRectilinearCoordinates(vtkImageData *imageOutput)
1228 {
1229 double origin[3];
1230 origin[0] = origin[1] = origin[2] = 0.0;
1231 double spacing[3];
1232 spacing[0] = spacing[1] = spacing[2] = 1.0;
1233
1234 int numDim = this->LoadingDimensions->GetNumberOfTuples();
1235 if (numDim >= 3) numDim = 3;
1236
1237 for (int i = 0; i < numDim; i++)
1238 {
1239 // Remember that netCDF dimension ordering is backward from VTK.
1240 int dim = this->LoadingDimensions->GetValue(numDim-i-1);
1241 vtkDimensionInfo *dimInfo = this->GetDimensionInfo(dim);
1242 origin[i] = dimInfo->GetOrigin();
1243 spacing[i] = dimInfo->GetSpacing();
1244 }
1245
1246 imageOutput->SetOrigin(origin);
1247 imageOutput->SetSpacing(spacing);
1248 }
1249
1250 //-----------------------------------------------------------------------------
AddRectilinearCoordinates(vtkRectilinearGrid * rectilinearOutput)1251 void vtkNetCDFCFReader::AddRectilinearCoordinates(
1252 vtkRectilinearGrid *rectilinearOutput)
1253 {
1254 int extent[6];
1255 rectilinearOutput->GetExtent(extent);
1256
1257 int numDim = this->LoadingDimensions->GetNumberOfTuples();
1258 for (int i = 0; i < 3; i++)
1259 {
1260 vtkSmartPointer<vtkDoubleArray> coords;
1261 if (i < numDim)
1262 {
1263 // Remember that netCDF dimension ordering is backward from VTK.
1264 int dim = this->LoadingDimensions->GetValue(numDim-i-1);
1265 coords = this->GetDimensionInfo(dim)->GetCoordinates();
1266 int extLow = extent[2*i];
1267 int extHi = extent[2*i+1];
1268 if ((extLow != 0) || (extHi != coords->GetNumberOfTuples()-1))
1269 {
1270 // Getting a subset of this dimension.
1271 VTK_CREATE(vtkDoubleArray, newcoords);
1272 newcoords->SetNumberOfComponents(1);
1273 newcoords->SetNumberOfTuples(extHi-extLow+1);
1274 memcpy(newcoords->GetPointer(0), coords->GetPointer(extLow),
1275 (extHi-extLow+1)*sizeof(double));
1276 coords = newcoords;
1277 }
1278 }
1279 else
1280 {
1281 coords = vtkSmartPointer<vtkDoubleArray>::New();
1282 coords->SetNumberOfTuples(1);
1283 coords->SetComponent(0, 0, 0.0);
1284 }
1285 switch (i)
1286 {
1287 case 0: rectilinearOutput->SetXCoordinates(coords); break;
1288 case 1: rectilinearOutput->SetYCoordinates(coords); break;
1289 case 2: rectilinearOutput->SetZCoordinates(coords); break;
1290 default: vtkErrorMacro("Sanity check failed!"); break;
1291 }
1292 }
1293 }
1294
1295 //-----------------------------------------------------------------------------
FakeRectilinearCoordinates(vtkRectilinearGrid * rectilinearOutput)1296 void vtkNetCDFCFReader::FakeRectilinearCoordinates(
1297 vtkRectilinearGrid *rectilinearOutput)
1298 {
1299 int dimensionSizes[3];
1300 rectilinearOutput->GetDimensions(dimensionSizes);
1301
1302 for (int dim = 0; dim < 3; dim++)
1303 {
1304 vtkNew<vtkDoubleArray> coordinate;
1305 coordinate->SetNumberOfComponents(1);
1306 coordinate->SetNumberOfTuples(dimensionSizes[dim]);
1307
1308 for (int index = 0; index < dimensionSizes[dim]; index++)
1309 {
1310 coordinate->SetComponent(index, 0, static_cast<double>(index));
1311 }
1312 switch(dim)
1313 {
1314 case 0: rectilinearOutput->SetXCoordinates(coordinate);break;
1315 case 1: rectilinearOutput->SetYCoordinates(coordinate);break;
1316 case 2: rectilinearOutput->SetZCoordinates(coordinate);break;
1317 default: vtkErrorMacro("Sanity check failed!"); break;
1318 }
1319 }
1320 }
1321
1322 //-----------------------------------------------------------------------------
Add1DRectilinearCoordinates(vtkPoints * points,const int extent[6])1323 void vtkNetCDFCFReader::Add1DRectilinearCoordinates(vtkPoints *points,
1324 const int extent[6])
1325 {
1326 points->SetDataTypeToDouble();
1327 points->SetNumberOfPoints( (extent[1]-extent[0]+1)
1328 * (extent[3]-extent[2]+1)
1329 * (extent[5]-extent[4]+1) );
1330 vtkDataArray *pointData = points->GetData();
1331
1332 int numDimNetCDF = this->LoadingDimensions->GetNumberOfTuples();
1333 for (int dimVTK = 0; dimVTK < 3; dimVTK++)
1334 {
1335 vtkSmartPointer<vtkDoubleArray> coords;
1336 if (dimVTK < numDimNetCDF)
1337 {
1338 // Remember that netCDF dimension ordering is backward from VTK.
1339 int dimNetCDF = this->LoadingDimensions->GetValue(numDimNetCDF-dimVTK-1);
1340 coords = this->GetDimensionInfo(dimNetCDF)->GetCoordinates();
1341
1342 int ijk[3];
1343 vtkIdType pointIdx = 0;
1344 for (ijk[2] = extent[2*2]; ijk[2] <= extent[2*2+1]; ijk[2]++)
1345 {
1346 for (ijk[1] = extent[1*2]; ijk[1] <= extent[1*2+1]; ijk[1]++)
1347 {
1348 for (ijk[0] = extent[0*2]; ijk[0] <= extent[0*2+1]; ijk[0]++)
1349 {
1350 pointData->SetComponent(pointIdx, dimVTK,
1351 coords->GetValue(ijk[dimVTK]));
1352 pointIdx++;
1353 }
1354 }
1355 }
1356 }
1357 else
1358 {
1359 int ijk[3];
1360 vtkIdType pointIdx = 0;
1361 for (ijk[2] = extent[2*2]; ijk[2] <= extent[2*2+1]; ijk[2]++)
1362 {
1363 for (ijk[1] = extent[1*2]; ijk[1] <= extent[1*2+1]; ijk[1]++)
1364 {
1365 for (ijk[0] = extent[0*2]; ijk[0] <= extent[0*2+1]; ijk[0]++)
1366 {
1367 pointData->SetComponent(pointIdx, dimVTK, 0.0);
1368 pointIdx++;
1369 }
1370 }
1371 }
1372 }
1373 }
1374 }
1375
1376 //-----------------------------------------------------------------------------
Add2DRectilinearCoordinates(vtkPoints * points,const int extent[6])1377 void vtkNetCDFCFReader::Add2DRectilinearCoordinates(vtkPoints *points,
1378 const int extent[6])
1379 {
1380 points->SetDataTypeToDouble();
1381 points->Allocate( (extent[1]-extent[0]+1)
1382 * (extent[3]-extent[2]+1)
1383 * (extent[5]-extent[4]+1) );
1384
1385 vtkDependentDimensionInfo *info
1386 = this->FindDependentDimensionInfo(this->LoadingDimensions);
1387
1388 vtkDoubleArray *longitudeCoordinates = info->GetLongitudeCoordinates();
1389 vtkDoubleArray *latitudeCoordinates = info->GetLatitudeCoordinates();
1390
1391 vtkDoubleArray *verticalCoordinates = nullptr;
1392 if (this->LoadingDimensions->GetNumberOfTuples() == 3)
1393 {
1394 int vertDim = this->LoadingDimensions->GetValue(0);
1395 if (info->GetHasBounds())
1396 {
1397 verticalCoordinates = this->GetDimensionInfo(vertDim)->GetBounds();
1398 }
1399 else
1400 {
1401 verticalCoordinates = this->GetDimensionInfo(vertDim)->GetCoordinates();
1402 }
1403 }
1404
1405 for (int k = extent[4]; k <= extent[5]; k++)
1406 {
1407 double h;
1408 if (verticalCoordinates)
1409 {
1410 h = verticalCoordinates->GetValue(k);
1411 }
1412 else
1413 {
1414 h = 0.0;
1415 }
1416 for (int j = extent[2]; j <= extent[3]; j++)
1417 {
1418 for (int i = extent[0]; i <= extent[1]; i++)
1419 {
1420 double lon = longitudeCoordinates->GetComponent(j, i);
1421 double lat = latitudeCoordinates->GetComponent(j, i);
1422 points->InsertNextPoint(lon, lat, h);
1423 }
1424 }
1425 }
1426 }
1427
1428 //-----------------------------------------------------------------------------
Add1DRectilinearCoordinates(vtkStructuredGrid * structuredOutput)1429 void vtkNetCDFCFReader::Add1DRectilinearCoordinates(
1430 vtkStructuredGrid *structuredOutput)
1431 {
1432 int extent[6];
1433 structuredOutput->GetExtent(extent);
1434
1435 VTK_CREATE(vtkPoints, points);
1436 this->Add1DRectilinearCoordinates(points, extent);
1437 structuredOutput->SetPoints(points);
1438 }
1439
1440 //-----------------------------------------------------------------------------
Add2DRectilinearCoordinates(vtkStructuredGrid * structuredOutput)1441 void vtkNetCDFCFReader::Add2DRectilinearCoordinates(
1442 vtkStructuredGrid *structuredOutput)
1443 {
1444 int extent[6];
1445 structuredOutput->GetExtent(extent);
1446
1447 VTK_CREATE(vtkPoints, points);
1448 this->Add2DRectilinearCoordinates(points, extent);
1449 structuredOutput->SetPoints(points);
1450 }
1451
1452 //-----------------------------------------------------------------------------
FakeStructuredCoordinates(vtkStructuredGrid * structuredOutput)1453 void vtkNetCDFCFReader::FakeStructuredCoordinates(
1454 vtkStructuredGrid *structuredOutput)
1455 {
1456 int extent[6];
1457 structuredOutput->GetExtent(extent);
1458
1459 vtkNew<vtkPoints> points;
1460 points->SetDataTypeToDouble();
1461 points->Allocate( (extent[1]-extent[0]+1)
1462 * (extent[3]-extent[2]+1)
1463 * (extent[5]-extent[4]+1) );
1464
1465 for (int kIndex = extent[4]; kIndex <= extent[5]; kIndex++)
1466 {
1467 for (int jIndex = extent[2]; jIndex <= extent[3]; jIndex++)
1468 {
1469 for (int iIndex = extent[0]; iIndex <= extent[1]; iIndex++)
1470 {
1471 points->InsertNextPoint(iIndex, jIndex, kIndex);
1472 }
1473 }
1474 }
1475
1476 structuredOutput->SetPoints(points);
1477 }
1478
1479 //-----------------------------------------------------------------------------
Add1DRectilinearCoordinates(vtkUnstructuredGrid * unstructuredOutput,const int extent[6])1480 void vtkNetCDFCFReader::Add1DRectilinearCoordinates(
1481 vtkUnstructuredGrid *unstructuredOutput,
1482 const int extent[6])
1483 {
1484 VTK_CREATE(vtkPoints, points);
1485 this->Add1DRectilinearCoordinates(points, extent);
1486 unstructuredOutput->SetPoints(points);
1487
1488 this->AddStructuredCells(unstructuredOutput, extent);
1489 }
1490
1491 //-----------------------------------------------------------------------------
Add2DRectilinearCoordinates(vtkUnstructuredGrid * unstructuredOutput,const int extent[6])1492 void vtkNetCDFCFReader::Add2DRectilinearCoordinates(
1493 vtkUnstructuredGrid *unstructuredOutput,
1494 const int extent[6])
1495 {
1496 VTK_CREATE(vtkPoints, points);
1497 this->Add2DRectilinearCoordinates(points, extent);
1498 unstructuredOutput->SetPoints(points);
1499
1500 this->AddStructuredCells(unstructuredOutput, extent);
1501 }
1502
1503 //-----------------------------------------------------------------------------
Add1DSphericalCoordinates(vtkPoints * points,const int extent[6])1504 void vtkNetCDFCFReader::Add1DSphericalCoordinates(vtkPoints *points,
1505 const int extent[6])
1506 {
1507 points->SetDataTypeToDouble();
1508 points->Allocate( (extent[1]-extent[0]+1)
1509 * (extent[3]-extent[2]+1)
1510 * (extent[5]-extent[4]+1) );
1511
1512 vtkDoubleArray *coordArrays[3];
1513 for (vtkIdType i = 0; i < this->LoadingDimensions->GetNumberOfTuples(); i++)
1514 {
1515 int dim = this->LoadingDimensions->GetValue(i);
1516 coordArrays[i] = this->GetDimensionInfo(dim)->GetBounds();
1517 }
1518
1519 int longitudeDim, latitudeDim, verticalDim;
1520 this->IdentifySphericalCoordinates(this->LoadingDimensions,
1521 longitudeDim, latitudeDim, verticalDim);
1522
1523 if ((longitudeDim < 0) || (latitudeDim < 0))
1524 {
1525 vtkErrorMacro(<< "Internal error: treating non spherical coordinates as if"
1526 << " they were spherical.");
1527 return;
1528 }
1529
1530 // Check the height scale and bias.
1531 double vertScale = this->VerticalScale;
1532 double vertBias = this->VerticalBias;
1533 if (verticalDim >= 0)
1534 {
1535 double *verticalRange = coordArrays[verticalDim]->GetRange();
1536 if ( (verticalRange[0]*vertScale + vertBias < 0)
1537 || (verticalRange[1]*vertScale + vertBias < 0) )
1538 {
1539 vertBias = -std::min(verticalRange[0], verticalRange[1])*vertScale;
1540 }
1541 }
1542 else
1543 {
1544 if (vertScale + vertBias <= 0)
1545 {
1546 vertScale = 1.0;
1547 vertBias = 0.0;
1548 }
1549 }
1550
1551 int ijk[3];
1552 for (ijk[0] = extent[4]; ijk[0] <= extent[5]; ijk[0]++)
1553 {
1554 for (ijk[1] = extent[2]; ijk[1] <= extent[3]; ijk[1]++)
1555 {
1556 for (ijk[2] = extent[0]; ijk[2] <= extent[1]; ijk[2]++)
1557 {
1558 double lon, lat, h;
1559 if (verticalDim >= 0)
1560 {
1561 lon = coordArrays[longitudeDim]->GetValue(ijk[longitudeDim]);
1562 lat = coordArrays[latitudeDim]->GetValue(ijk[latitudeDim]);
1563 h = coordArrays[verticalDim]->GetValue(ijk[verticalDim]);
1564 }
1565 else
1566 {
1567 lon = coordArrays[longitudeDim]->GetValue(ijk[longitudeDim+1]);
1568 lat = coordArrays[latitudeDim]->GetValue(ijk[latitudeDim+1]);
1569 h = 1.0;
1570 }
1571 lon = vtkMath::RadiansFromDegrees(lon);
1572 lat = vtkMath::RadiansFromDegrees(lat);
1573 h = h*vertScale + vertBias;
1574
1575 double cartesianCoord[3];
1576 cartesianCoord[0] = h*cos(lon)*cos(lat);
1577 cartesianCoord[1] = h*sin(lon)*cos(lat);
1578 cartesianCoord[2] = h*sin(lat);
1579 points->InsertNextPoint(cartesianCoord);
1580 }
1581 }
1582 }
1583 }
1584
1585 //-----------------------------------------------------------------------------
Add2DSphericalCoordinates(vtkPoints * points,const int extent[6])1586 void vtkNetCDFCFReader::Add2DSphericalCoordinates(vtkPoints *points,
1587 const int extent[6])
1588 {
1589 points->SetDataTypeToDouble();
1590 points->Allocate( (extent[1]-extent[0]+1)
1591 * (extent[3]-extent[2]+1)
1592 * (extent[5]-extent[4]+1) );
1593
1594 vtkDependentDimensionInfo *info
1595 = this->FindDependentDimensionInfo(this->LoadingDimensions);
1596
1597 vtkDoubleArray *longitudeCoordinates = info->GetLongitudeCoordinates();
1598 vtkDoubleArray *latitudeCoordinates = info->GetLatitudeCoordinates();
1599
1600 vtkDoubleArray *verticalCoordinates = nullptr;
1601 if (this->LoadingDimensions->GetNumberOfTuples() == 3)
1602 {
1603 int vertDim = this->LoadingDimensions->GetValue(0);
1604 if (info->GetHasBounds())
1605 {
1606 verticalCoordinates = this->GetDimensionInfo(vertDim)->GetBounds();
1607 }
1608 else
1609 {
1610 verticalCoordinates = this->GetDimensionInfo(vertDim)->GetCoordinates();
1611 }
1612 }
1613
1614 // Check the height scale and bias.
1615 double vertScale = this->VerticalScale;
1616 double vertBias = this->VerticalBias;
1617 if (verticalCoordinates)
1618 {
1619 double *verticalRange = verticalCoordinates->GetRange();
1620 if ( (verticalRange[0]*vertScale + vertBias < 0)
1621 || (verticalRange[1]*vertScale + vertBias < 0) )
1622 {
1623 vertBias = -std::min(verticalRange[0], verticalRange[1])*vertScale;
1624 }
1625 }
1626 else
1627 {
1628 if (vertScale + vertBias <= 0)
1629 {
1630 vertScale = 1.0;
1631 vertBias = 0.0;
1632 }
1633 }
1634
1635 for (int k = extent[4]; k <= extent[5]; k++)
1636 {
1637 double h;
1638 if (verticalCoordinates)
1639 {
1640 h = verticalCoordinates->GetValue(k)*vertScale + vertBias;
1641 }
1642 else
1643 {
1644 h = vertScale + vertBias;
1645 }
1646 for (int j = extent[2]; j <= extent[3]; j++)
1647 {
1648 for (int i = extent[0]; i <= extent[1]; i++)
1649 {
1650 double lon = longitudeCoordinates->GetComponent(j, i);
1651 double lat = latitudeCoordinates->GetComponent(j, i);
1652 lon = vtkMath::RadiansFromDegrees(lon);
1653 lat = vtkMath::RadiansFromDegrees(lat);
1654
1655 double cartesianCoord[3];
1656 cartesianCoord[0] = h*cos(lon)*cos(lat);
1657 cartesianCoord[1] = h*sin(lon)*cos(lat);
1658 cartesianCoord[2] = h*sin(lat);
1659 points->InsertNextPoint(cartesianCoord);
1660 }
1661 }
1662 }
1663 }
1664
1665 //-----------------------------------------------------------------------------
Add1DSphericalCoordinates(vtkStructuredGrid * structuredOutput)1666 void vtkNetCDFCFReader::Add1DSphericalCoordinates(
1667 vtkStructuredGrid *structuredOutput)
1668 {
1669 int extent[6];
1670 structuredOutput->GetExtent(extent);
1671
1672 VTK_CREATE(vtkPoints, points);
1673 this->Add1DSphericalCoordinates(points, extent);
1674 structuredOutput->SetPoints(points);
1675 }
1676
1677 //-----------------------------------------------------------------------------
Add2DSphericalCoordinates(vtkStructuredGrid * structuredOutput)1678 void vtkNetCDFCFReader::Add2DSphericalCoordinates(
1679 vtkStructuredGrid *structuredOutput)
1680 {
1681 int extent[6];
1682 structuredOutput->GetExtent(extent);
1683
1684 VTK_CREATE(vtkPoints, points);
1685 this->Add2DSphericalCoordinates(points, extent);
1686 structuredOutput->SetPoints(points);
1687 }
1688
1689 //-----------------------------------------------------------------------------
Add1DSphericalCoordinates(vtkUnstructuredGrid * unstructuredOutput,const int extent[6])1690 void vtkNetCDFCFReader::Add1DSphericalCoordinates(
1691 vtkUnstructuredGrid *unstructuredOutput,
1692 const int extent[6])
1693 {
1694 VTK_CREATE(vtkPoints, points);
1695 this->Add1DSphericalCoordinates(points, extent);
1696 unstructuredOutput->SetPoints(points);
1697
1698 this->AddStructuredCells(unstructuredOutput, extent);
1699 }
1700
1701 //-----------------------------------------------------------------------------
Add2DSphericalCoordinates(vtkUnstructuredGrid * unstructuredOutput,const int extent[6])1702 void vtkNetCDFCFReader::Add2DSphericalCoordinates(
1703 vtkUnstructuredGrid *unstructuredOutput,
1704 const int extent[6])
1705 {
1706 VTK_CREATE(vtkPoints, points);
1707 this->Add2DSphericalCoordinates(points, extent);
1708 unstructuredOutput->SetPoints(points);
1709
1710 this->AddStructuredCells(unstructuredOutput, extent);
1711 }
1712
1713 //-----------------------------------------------------------------------------
AddStructuredCells(vtkUnstructuredGrid * unstructuredOutput,const int extent[6])1714 void vtkNetCDFCFReader::AddStructuredCells(
1715 vtkUnstructuredGrid *unstructuredOutput,
1716 const int extent[6])
1717 {
1718 vtkIdType numPoints[3];
1719 numPoints[0] = extent[1] - extent[0] + 1;
1720 numPoints[1] = extent[3] - extent[2] + 1;
1721 numPoints[2] = extent[5] - extent[4] + 1;
1722
1723 vtkIdType numCells[3];
1724 numCells[0] = numPoints[0] - 1;
1725 numCells[1] = numPoints[1] - 1;
1726 numCells[2] = numPoints[2] - 1;
1727
1728 vtkIdType nextPointRow = numPoints[0];
1729 vtkIdType nextPointSlab = nextPointRow*numPoints[1];
1730
1731 bool extentIs2D = (numCells[2] < 1);
1732
1733 if (extentIs2D)
1734 {
1735 vtkIdType totalNumCells = numCells[0]*numCells[1];
1736 unstructuredOutput->Allocate(totalNumCells);
1737 vtkCellArray *cells = unstructuredOutput->GetCells();
1738 cells->Allocate(cells->EstimateSize(totalNumCells, 4));
1739
1740 for (int j = 0; j < numCells[1]; j++)
1741 {
1742 vtkIdType rowStart = j*nextPointRow;
1743 for (int i = 0; i < numCells[0]; i++)
1744 {
1745 vtkIdType lowCellPoint = rowStart + i;
1746
1747 vtkIdType pointIds[4];
1748 pointIds[0] = lowCellPoint;
1749 pointIds[1] = lowCellPoint + 1;
1750 pointIds[2] = lowCellPoint + nextPointRow + 1;
1751 pointIds[3] = lowCellPoint + nextPointRow;
1752
1753 unstructuredOutput->InsertNextCell(VTK_QUAD, 4, pointIds);
1754 }
1755 }
1756 }
1757 else // !extentIs2D
1758 {
1759 vtkIdType totalNumCells = numCells[0]*numCells[1]*numCells[2];
1760 unstructuredOutput->Allocate(totalNumCells);
1761 vtkCellArray *cells = unstructuredOutput->GetCells();
1762 cells->Allocate(cells->EstimateSize(totalNumCells, 8));
1763
1764 for (int k = 0; k < numCells[2]; k++)
1765 {
1766 vtkIdType slabStart = k*nextPointSlab;
1767 for (int j = 0; j < numCells[1]; j++)
1768 {
1769 vtkIdType rowStart = slabStart + j*nextPointRow;
1770 for (int i = 0; i < numCells[0]; i++)
1771 {
1772 vtkIdType lowCellPoint = rowStart + i;
1773
1774 // This code is assuming that all axis are scaling up. If that is
1775 // not the case, this will probably make inverted hexahedra.
1776 vtkIdType pointIds[8];
1777 pointIds[0] = lowCellPoint;
1778 pointIds[1] = lowCellPoint + 1;
1779 pointIds[2] = lowCellPoint + nextPointRow + 1;
1780 pointIds[3] = lowCellPoint + nextPointRow;
1781 pointIds[4] = lowCellPoint + nextPointSlab;
1782 pointIds[5] = lowCellPoint + nextPointSlab + 1;
1783 pointIds[6] = lowCellPoint + nextPointSlab + nextPointRow + 1;
1784 pointIds[7] = lowCellPoint + nextPointSlab + nextPointRow;
1785
1786 unstructuredOutput->InsertNextCell(VTK_HEXAHEDRON, 8, pointIds);
1787 }
1788 }
1789 }
1790 }
1791 }
1792
1793 //-----------------------------------------------------------------------------
AddUnstructuredRectilinearCoordinates(vtkUnstructuredGrid * unstructuredOutput,const int extent[6])1794 void vtkNetCDFCFReader::AddUnstructuredRectilinearCoordinates(
1795 vtkUnstructuredGrid *unstructuredOutput,
1796 const int extent[6])
1797 {
1798 vtkDependentDimensionInfo *info
1799 = this->FindDependentDimensionInfo(this->LoadingDimensions);
1800
1801 vtkDoubleArray *longitudeCoordinates = info->GetLongitudeCoordinates();
1802 vtkDoubleArray *latitudeCoordinates = info->GetLatitudeCoordinates();
1803
1804 int numPointsPerCell = longitudeCoordinates->GetNumberOfComponents();
1805 vtkIdType totalNumCells = longitudeCoordinates->GetNumberOfTuples();
1806
1807 double bounds[6];
1808 GetRangeOfAllComponents(longitudeCoordinates, bounds+0);
1809 GetRangeOfAllComponents(latitudeCoordinates, bounds+2);
1810 bounds[4] = bounds[5] = 0.0;
1811
1812 VTK_CREATE(vtkPoints, points);
1813 points->SetDataTypeToDouble();
1814 points->Allocate(totalNumCells);
1815
1816 VTK_CREATE(vtkMergePoints, locator);
1817 locator->InitPointInsertion(points, bounds);
1818
1819 // Make space in output unstructured grid.
1820 unstructuredOutput->Allocate(extent[1] - extent[0]);
1821 vtkCellArray *cells = unstructuredOutput->GetCells();
1822 cells->Allocate(cells->EstimateSize(extent[1]-extent[0], numPointsPerCell));
1823
1824 std::vector<vtkIdType> cellPoints(numPointsPerCell);
1825
1826 // This is a rather lame way to break up cells amongst processes. It will be
1827 // slow and ghost cells are totally screwed up.
1828 for (int cellId = extent[0]; cellId < extent[1]; cellId++)
1829 {
1830 for (int cellPointId = 0; cellPointId < numPointsPerCell; cellPointId++)
1831 {
1832 double coord[3];
1833 coord[0] = longitudeCoordinates->GetComponent(cellId, cellPointId);
1834 coord[1] = latitudeCoordinates->GetComponent(cellId, cellPointId);
1835 coord[2] = 0.0;
1836
1837 vtkIdType pointId;
1838 locator->InsertUniquePoint(coord, pointId);
1839
1840 cellPoints[cellPointId] = pointId;
1841 }
1842 unstructuredOutput
1843 ->InsertNextCell(VTK_POLYGON, numPointsPerCell, &cellPoints.at(0));
1844 }
1845
1846 points->Squeeze();
1847 unstructuredOutput->SetPoints(points);
1848 }
1849
1850 //-----------------------------------------------------------------------------
AddUnstructuredSphericalCoordinates(vtkUnstructuredGrid * unstructuredOutput,const int extent[6])1851 void vtkNetCDFCFReader::AddUnstructuredSphericalCoordinates(
1852 vtkUnstructuredGrid *unstructuredOutput,
1853 const int extent[6])
1854 {
1855 // First load the data as rectilinear coordinates, and then convert them
1856 // to spherical coordinates. Not only does this reuse code, but it also
1857 // probably makes the locator more efficient this way.
1858 this->AddUnstructuredRectilinearCoordinates(unstructuredOutput, extent);
1859
1860 double height = 1.0*this->VerticalScale + this->VerticalBias;
1861 if (height <= 0.0)
1862 {
1863 height = 1.0;
1864 }
1865
1866 vtkPoints *points = unstructuredOutput->GetPoints();
1867 vtkIdType numPoints = points->GetNumberOfPoints();
1868 for (vtkIdType pointId = 0; pointId < numPoints; pointId++)
1869 {
1870 double lonLat[3];
1871 points->GetPoint(pointId, lonLat);
1872 double lon = vtkMath::RadiansFromDegrees(lonLat[0]);
1873 double lat = vtkMath::RadiansFromDegrees(lonLat[1]);
1874
1875 double cartesianCoord[3];
1876 cartesianCoord[0] = height*cos(lon)*cos(lat);
1877 cartesianCoord[1] = height*sin(lon)*cos(lat);
1878 cartesianCoord[2] = height*sin(lat);
1879 points->SetPoint(pointId, cartesianCoord);
1880 }
1881 }
1882
1883 //-----------------------------------------------------------------------------
ReadMetaData(int ncFD)1884 int vtkNetCDFCFReader::ReadMetaData(int ncFD)
1885 {
1886 vtkDebugMacro("ReadMetaData");
1887
1888 int numDimensions;
1889 CALL_NETCDF(nc_inq_ndims(ncFD, &numDimensions));
1890 this->DimensionInfo->v.resize(numDimensions);
1891
1892 std::set<vtkStdString> specialVariables;
1893
1894 for (int i = 0; i < numDimensions; i++)
1895 {
1896 this->DimensionInfo->v[i] = vtkDimensionInfo(ncFD, i);
1897
1898 // Record any special variables for this dimension.
1899 vtkStringArray* dimensionVariables
1900 = this->DimensionInfo->v[i].GetSpecialVariables();
1901 for (vtkIdType j = 0; j < dimensionVariables->GetNumberOfValues(); j++)
1902 {
1903 specialVariables.insert(dimensionVariables->GetValue(j));
1904 }
1905 }
1906
1907 int numVariables;
1908 CALL_NETCDF(nc_inq_nvars(ncFD, &numVariables));
1909
1910 // Check all variables for special 2D coordinates.
1911 for (int i = 0; i < numVariables; i++)
1912 {
1913 vtkDependentDimensionInfo info(ncFD, i, this);
1914 if (!info.GetValid()) continue;
1915 if (this->FindDependentDimensionInfo(info.GetGridDimensions()) != nullptr)
1916 {
1917 continue;
1918 }
1919
1920 this->DependentDimensionInfo->v.push_back(info);
1921
1922 // Record any special variables.
1923 vtkStringArray* dimensionVariables = info.GetSpecialVariables();
1924 for (vtkIdType j = 0; j < dimensionVariables->GetNumberOfValues(); j++)
1925 {
1926 specialVariables.insert(dimensionVariables->GetValue(j));
1927 }
1928 }
1929
1930 // Look at all variables and record them so that the user can select which
1931 // ones he wants. This oddness of adding and removing from
1932 // VariableArraySelection is to preserve any current settings for variables.
1933 typedef std::set<vtkStdString> stringSet;
1934 stringSet variablesToAdd;
1935 stringSet variablesToRemove;
1936
1937 // Initialize variablesToRemove with all the variables. Then remove them from
1938 // the list as we find them.
1939 for (int i = 0; i < this->VariableArraySelection->GetNumberOfArrays(); i++)
1940 {
1941 variablesToRemove.insert(this->VariableArraySelection->GetArrayName(i));
1942 }
1943
1944 for (int i = 0; i < numVariables; i++)
1945 {
1946 char name[NC_MAX_NAME+1];
1947 CALL_NETCDF(nc_inq_varname(ncFD, i, name));
1948 if (specialVariables.find(name) == specialVariables.end())
1949 {
1950 if (variablesToRemove.find(name) == variablesToRemove.end())
1951 {
1952 // Variable not already here. Insert it in the variables to add.
1953 variablesToAdd.insert(name);
1954 }
1955 else
1956 {
1957 // Variable already exists. Leave it be. Remove it from the
1958 // variablesToRemove list.
1959 variablesToRemove.erase(name);
1960 }
1961 }
1962 }
1963
1964 // Add and remove variables. This will be a no-op if the variables have not
1965 // changed.
1966 for (stringSet::iterator removeItr = variablesToRemove.begin();
1967 removeItr != variablesToRemove.end(); ++removeItr)
1968 {
1969 this->VariableArraySelection->RemoveArrayByName(removeItr->c_str());
1970 }
1971 for (stringSet::iterator addItr = variablesToAdd.begin();
1972 addItr != variablesToAdd.end(); ++addItr)
1973 {
1974 this->VariableArraySelection->AddArray(addItr->c_str());
1975 }
1976
1977 return 1;
1978 }
1979
1980 //-----------------------------------------------------------------------------
IsTimeDimension(int vtkNotUsed (ncFD),int dimId)1981 int vtkNetCDFCFReader::IsTimeDimension(int vtkNotUsed(ncFD), int dimId)
1982 {
1983 return ( this->GetDimensionInfo(dimId)->GetUnits()
1984 == vtkDimensionInfo::TIME_UNITS );
1985 }
1986
1987 //-----------------------------------------------------------------------------
GetTimeValues(int vtkNotUsed (ncFD),int dimId)1988 vtkSmartPointer<vtkDoubleArray> vtkNetCDFCFReader::GetTimeValues(
1989 int vtkNotUsed(ncFD), int dimId)
1990 {
1991 return this->GetDimensionInfo(dimId)->GetCoordinates();
1992 }
1993
1994 //-----------------------------------------------------------------------------
1995 inline vtkNetCDFCFReader::vtkDimensionInfo *
GetDimensionInfo(int dimension)1996 vtkNetCDFCFReader::GetDimensionInfo(int dimension)
1997 {
1998 return &(this->DimensionInfo->v.at(dimension));
1999 }
2000
2001 //-----------------------------------------------------------------------------
2002 vtkNetCDFCFReader::vtkDependentDimensionInfo *
FindDependentDimensionInfo(vtkIntArray * dims)2003 vtkNetCDFCFReader::FindDependentDimensionInfo(vtkIntArray *dims)
2004 {
2005 for (size_t i = 0; i < this->DependentDimensionInfo->v.size(); i++)
2006 {
2007 vtkIntArray *dependentDims
2008 = this->DependentDimensionInfo->v[i].GetGridDimensions();
2009 if (dims->GetNumberOfTuples() == dependentDims->GetNumberOfTuples())
2010 {
2011 bool same = true;
2012 for (vtkIdType j = 0; j < dims->GetNumberOfTuples(); j++)
2013 {
2014 if (dims->GetValue(j) != dependentDims->GetValue(j))
2015 {
2016 same = false;
2017 break;
2018 }
2019 }
2020 if (same) return &(this->DependentDimensionInfo->v[i]);
2021 }
2022 }
2023 return nullptr;
2024 }
2025
2026 //-----------------------------------------------------------------------------
IdentifySphericalCoordinates(vtkIntArray * dimensions,int & longitudeDim,int & latitudeDim,int & verticalDim)2027 void vtkNetCDFCFReader::IdentifySphericalCoordinates(vtkIntArray *dimensions,
2028 int &longitudeDim,
2029 int &latitudeDim,
2030 int &verticalDim)
2031 {
2032 longitudeDim = latitudeDim = verticalDim = -1;
2033 for (int i = 0; i < dimensions->GetNumberOfTuples(); i++)
2034 {
2035 switch (this->GetDimensionInfo(dimensions->GetValue(i))->GetUnits())
2036 {
2037 case vtkDimensionInfo::LONGITUDE_UNITS:
2038 longitudeDim = i;
2039 break;
2040 case vtkDimensionInfo::LATITUDE_UNITS:
2041 latitudeDim = i;
2042 break;
2043 case vtkDimensionInfo::UNDEFINED_UNITS:
2044 case vtkDimensionInfo::TIME_UNITS:
2045 case vtkDimensionInfo::VERTICAL_UNITS:
2046 default:
2047 verticalDim = i;
2048 break;
2049 }
2050 }
2051 }
2052
2053 //-----------------------------------------------------------------------------
2054 vtkNetCDFCFReader::CoordinateTypesEnum
CoordinateType(vtkIntArray * dimensions)2055 vtkNetCDFCFReader::CoordinateType(vtkIntArray *dimensions)
2056 {
2057 vtkDependentDimensionInfo *dependentDimInfo
2058 = this->FindDependentDimensionInfo(dimensions);
2059
2060 // Check to see if using p-sided cells.
2061 if (dependentDimInfo && dependentDimInfo->GetCellsUnstructured())
2062 {
2063 if (this->SphericalCoordinates)
2064 {
2065 return COORDS_SPHERICAL_PSIDED_CELLS;
2066 }
2067 else
2068 {
2069 return COORDS_EUCLIDEAN_PSIDED_CELLS;
2070 }
2071 }
2072
2073 // Check to see if using 4-sided cells.
2074 if ( dependentDimInfo
2075 && !dependentDimInfo->GetCellsUnstructured()
2076 && dependentDimInfo->GetHasBounds() )
2077 {
2078 if (this->SphericalCoordinates)
2079 {
2080 return COORDS_SPHERICAL_4SIDED_CELLS;
2081 }
2082 else
2083 {
2084 return COORDS_EUCLIDEAN_4SIDED_CELLS;
2085 }
2086 }
2087
2088 // Check to see if using 2D coordinate lookup.
2089 if ( dependentDimInfo
2090 && !dependentDimInfo->GetCellsUnstructured()
2091 && !dependentDimInfo->GetHasBounds() )
2092 {
2093 if (this->SphericalCoordinates)
2094 {
2095 return COORDS_2D_SPHERICAL;
2096 }
2097 else
2098 {
2099 return COORDS_2D_EUCLIDEAN;
2100 }
2101 }
2102
2103 // Check to see if we should (otherwise) be using spherical coordinates.
2104 if (this->SphericalCoordinates)
2105 {
2106 int longitudeDim, latitudeDim, verticalDim;
2107 this->IdentifySphericalCoordinates(dimensions,
2108 longitudeDim,
2109 latitudeDim,
2110 verticalDim);
2111 if ( (longitudeDim != -1) && (latitudeDim != -1)
2112 && ((dimensions->GetNumberOfTuples() == 2) || (verticalDim != -1)) )
2113 {
2114 return COORDS_REGULAR_SPHERICAL;
2115 }
2116 }
2117
2118 // Check to see if any dimension as irregular spacing.
2119 for (int i = 0; i < dimensions->GetNumberOfTuples(); i++)
2120 {
2121 int dimId = dimensions->GetValue(i);
2122 if (!this->GetDimensionInfo(dimId)->GetHasRegularSpacing())
2123 {
2124 return COORDS_NONUNIFORM_RECTILINEAR;
2125 }
2126 }
2127
2128 // All dimensions appear to be uniform rectilinear.
2129 return COORDS_UNIFORM_RECTILINEAR;
2130 }
2131
2132 //-----------------------------------------------------------------------------
DimensionsAreForPointData(vtkIntArray * dimensions)2133 bool vtkNetCDFCFReader::DimensionsAreForPointData(vtkIntArray *dimensions)
2134 {
2135 switch (this->CoordinateType(dimensions))
2136 {
2137 case COORDS_UNIFORM_RECTILINEAR: return true;
2138 case COORDS_NONUNIFORM_RECTILINEAR: return true;
2139 case COORDS_REGULAR_SPHERICAL: return false;
2140 case COORDS_2D_EUCLIDEAN: return true;
2141 case COORDS_2D_SPHERICAL: return true;
2142 case COORDS_EUCLIDEAN_4SIDED_CELLS: return false;
2143 case COORDS_SPHERICAL_4SIDED_CELLS: return false;
2144 case COORDS_EUCLIDEAN_PSIDED_CELLS: return false;
2145 case COORDS_SPHERICAL_PSIDED_CELLS: return false;
2146 default:
2147 vtkErrorMacro("Internal error: unknown coordinate type.");
2148 return true;
2149 }
2150 }
2151