diff --git a/Plugins/CDIReader/Reader/CMakeLists.txt b/Plugins/CDIReader/Reader/CMakeLists.txt index 89c6722ed3e4cdc1920efc9fae0766c5153ae602..717990b6a4b6a051be69569e8b51764e1684e973 100644 --- a/Plugins/CDIReader/Reader/CMakeLists.txt +++ b/Plugins/CDIReader/Reader/CMakeLists.txt @@ -4,6 +4,7 @@ set(classes set (sources cdi_tools.cxx projections.cxx + DataSource.cxx ) set (private_headers cdi_tools.h diff --git a/Plugins/CDIReader/Reader/DataSource.cxx b/Plugins/CDIReader/Reader/DataSource.cxx new file mode 100644 index 0000000000000000000000000000000000000000..4f5075bb58215ca77549eff9fc117508ee4f8513 --- /dev/null +++ b/Plugins/CDIReader/Reader/DataSource.cxx @@ -0,0 +1,56 @@ +//#include "cdi_tools.h" +#include "DataSource.h" +#include "cdi.h" + +namespace DataSource +{ + +CDIObject::CDIObject(std::string newURI) +{ + this->openURI(newURI); +} + +int CDIObject::openURI(std::string newURI) +{ + this->setVoid(); + this->URI = newURI; + + // check if we got either *.Grib or *.nc data + std::string check = this->URI.substr((URI.size() - 4), this->URI.size()); + if (check == "grib" || check == ".grb") + { + this->type = GRIB; + } + else + { + this->type = NC; + } + + this->StreamID = streamOpenRead(this->URI.c_str()); + if (this->StreamID < 0) + { + this->setVoid(); + return 0; + } + + this->VListID = streamInqVlist(this->StreamID); + return 1; +} + +void CDIObject::setVoid() +{ + if (this->StreamID > -1) + streamClose(this->StreamID); + + StreamID = -1; + VListID = -1; + type = VOID; +} + +CDIObject::~CDIObject() +{ + URI = ""; + this->setVoid(); +} + +}; diff --git a/Plugins/CDIReader/Reader/DataSource.h b/Plugins/CDIReader/Reader/DataSource.h new file mode 100644 index 0000000000000000000000000000000000000000..f3749f2f211c3444a9b4f8902dae6fe3656ef9ca --- /dev/null +++ b/Plugins/CDIReader/Reader/DataSource.h @@ -0,0 +1,39 @@ +#ifndef CDI_DATA_SOURCE +#define CDI_DATA_SOURCE + +#include + +namespace DataSource +{ +class CDIObject +{ + enum stype + { + VOID, + GRIB, + NC + }; + std::string URI; + int StreamID; + int VListID; + stype type; + +public: + void setVoid(); + CDIObject(std::string URI); + CDIObject() + { + this->StreamID = -1; + this->setVoid(); + } + ~CDIObject(); + int openURI(std::string URI); + std::string getURI() const { return URI; } + int getStreamID() const { return StreamID; } + int getVListID() const { return VListID; } + stype getType() const { return type; } + bool isVoid() const { return this->type == VOID; } +}; + +}; +#endif // CDI_DATA_SOURCE diff --git a/Plugins/CDIReader/Reader/vtkCDIReader.cxx b/Plugins/CDIReader/Reader/vtkCDIReader.cxx index 34335bd35f950721b5e89269a88f2d8f037b309b..be04cf478ea6a7ca0482702855c18fc383a864aa 100644 --- a/Plugins/CDIReader/Reader/vtkCDIReader.cxx +++ b/Plugins/CDIReader/Reader/vtkCDIReader.cxx @@ -36,6 +36,7 @@ #include "vtksys/FStream.hxx" #include "vtksys/SystemTools.hxx" +#include "DataSource.h" #include "cdi_tools.h" #include @@ -155,8 +156,6 @@ vtkCDIReader::vtkCDIReader() this->SetNumberOfInputPorts(0); this->SetNumberOfOutputPorts(1); - this->StreamID = -1; - this->VListID = -1; this->VariableDimensions = vtkSmartPointer::New(); this->AllDimensions = vtkSmartPointer::New(); this->AllVariableArrayNames = vtkSmartPointer::New(); @@ -210,13 +209,6 @@ void vtkCDIReader::DestroyData() vtkCDIReader::~vtkCDIReader() { vtkDebugMacro("Destructing vtkCDIReader..."); - this->SetFileName(nullptr); - - if (this->StreamID >= 0) - { - streamClose(this->StreamID); - this->StreamID = -1; - } this->DestroyData(); @@ -286,7 +278,11 @@ int vtkCDIReader::RequestInformation( vtkDebugMacro("In vtkCDIReader::RequestInformation setting VerticalLevelRange"); this->VerticalLevelRange[0] = 0; - this->VerticalLevelRange[1] = this->MaximumNVertLevels - 1; + if (VerticalLevelRange[1] != this->MaximumNVertLevels - 1) + { + this->VerticalLevelRange[1] = this->MaximumNVertLevels - 1; + this->Modified(); + } if (!this->BuildVarArrays()) { @@ -313,9 +309,9 @@ vtkSmartPointer vtkCDIReader::ReadTimeAxis() if ((this->FileSeriesNumber == 0) && (!this->TimeSet)) { - int taxisID = vlistInqTaxis(this->VListID); + int taxisID = vlistInqTaxis(this->DataFile.getVListID()); int calendar = taxisInqCalendar(taxisID); - streamInqTimestep(this->StreamID, 0); + streamInqTimestep(this->DataFile.getStreamID(), 0); int vdate = taxisInqVdate(taxisID); int vtime = taxisInqVtime(taxisID); @@ -350,9 +346,9 @@ vtkSmartPointer vtkCDIReader::ReadTimeAxis() int end = start + this->NumberOfTimeSteps; for (int step = start; step < end; step++) { - int taxisID = vlistInqTaxis(this->VListID); + int taxisID = vlistInqTaxis(this->DataFile.getVListID()); int calendar = taxisInqCalendar(taxisID); - streamInqTimestep(this->StreamID, counter); + streamInqTimestep(this->DataFile.getStreamID(), counter); int vdate = taxisInqVdate(taxisID); int vtime = taxisInqVtime(taxisID); double timevalue = date_to_julday(calendar, vdate); @@ -748,9 +744,8 @@ void vtkCDIReader::SetDefaults() this->DTime = 0; this->FileSeriesNumber = 0; this->NumberOfFiles = 1; - this->NeedHorizontalGridFile = false; this->NeedVerticalGridFile = false; - + this->GridID = -1; this->NumberOfProcesses = 1; this->BuildDomainArrays = false; @@ -765,50 +760,14 @@ void vtkCDIReader::SetDefaults() //---------------------------------------------------------------------------- // Get dimensions of key NetCDF variables //---------------------------------------------------------------------------- -int vtkCDIReader::OpenFile() -{ - // check if we got either *.Grib or *.nc data - std::string file = this->FileName; - std::string check = file.substr((file.size() - 4), file.size()); - if (check == "grib" || check == ".grb") - { - this->Grib = true; - } - else - { - this->Grib = false; - } - - if (this->StreamID >= 0) - { - streamClose(this->StreamID); - this->StreamID = -1; - this->VListID = -1; - } - - this->StreamID = streamOpenRead(this->FileNameGrid.c_str()); - if (this->StreamID < 0) - { - return 0; - } - - vtkDebugMacro("In vtkCDIReader::RequestInformation read file okay"); - this->VListID = streamInqVlist(this->StreamID); - - int nvars = vlistNvars(this->VListID); - char varname[CDI_MAX_NAME]; - for (int varID = 0; varID < nvars; ++varID) - { - vlistInqVarName(this->VListID, varID, varname); - } - - return 1; -} //---------------------------------------------------------------------------- void vtkCDIReader::GuessGridFile() { - std::string fallback = vtksys::SystemTools::GetParentDirectory(this->FileName) + "/grid.nc"; + std::string fallback = vtksys::SystemTools::GetParentDirectory(this->FileName); + if (fallback.empty()) + fallback = "."; + fallback += "/grid.nc"; std::string guess; if (!this->Grib) @@ -818,15 +777,22 @@ void vtkCDIReader::GuessGridFile() { if (vtksys::SystemTools::TestFileAccess(guess, vtksys::TEST_FILE_READ)) { - this->FileNameGrid = guess; - return; + this->GridFile.openURI(guess); + if (this->GridFile.isVoid()) + { + vtkWarningMacro("Cannot handle grid file " + << guess << " indicated by grid_file_uri attribute in " << this->FileName + << " Trying fallback guess " << fallback); + } + else + return; } else - vtkWarningMacro("Could not find grid file " + vtkWarningMacro("Cannot open grid file " << guess << " indicated by grid_file_uri attribute in " << this->FileName << " Trying fallback guess " << fallback); } - this->FileNameGrid = fallback; + this->GridFile.openURI(fallback); } //---------------------------------------------------------------------------- @@ -834,147 +800,132 @@ void vtkCDIReader::GuessGridFile() //---------------------------------------------------------------------------- int vtkCDIReader::GetDims() { - if (!this->FileName.empty()) + if (this->FileName.empty()) { - this->FileNameGrid = this->FileName; - if (this->VListID < 0 || this->StreamID < 0) - { - if (!this->OpenFile()) - { - return 0; - } - } - - this->ReadHorizontalGridData(); - if (this->NeedHorizontalGridFile) - { - // if there is no grid information in the data file, try opening - // an additional grid file named grid.nc in the same directory to - // read in the grid information - if (this->StreamID >= 0) - { - streamClose(this->StreamID); - this->StreamID = -1; - this->VListID = -1; - } + vtkErrorMacro("No file name provided. Cannot get dimensions"); + return 0; + } - char* directory = new char[strlen(this->FileName.c_str()) + 1]; - strcpy(directory, this->FileName.c_str()); + DataFile.openURI(FileName); + if (DataFile.isVoid()) + { + vtkErrorMacro("GetDims: Could not open " << DataFile.getURI()); + return 0; + } - this->GuessGridFile(); - if (!this->OpenFile()) - { - return 0; - } - if (!this->ReadHorizontalGridData()) - { - vtkErrorMacro("Couldn't open grid information in data nor in the grid file."); - return 0; - } + if (GridFile.isVoid()) + GridFile.openURI(FileName); + if (GridFile.isVoid()) + { + vtkErrorMacro("GetDims: Could not open horizontal grid file.\nTried " << GridFile.getURI()); + return 0; + } - this->FileNameGrid = this->FileName; - if (!this->OpenFile()) - { - return 0; - } + if (!this->ReadHorizontalGridData()) + { + this->GuessGridFile(); + if (!this->ReadHorizontalGridData()) + { + vtkErrorMacro("Could not get horizontal Grid. \nTried " << GridFile.getURI()); + return 0; } + } - this->ReadVerticalGridData(); - if (this->NeedVerticalGridFile) - { - // if there is no grid information in the data file, try opening - // an additional grid file named grid.nc in the same directory to - // read in the grid information - if (this->StreamID >= 0) - { - streamClose(this->StreamID); - this->StreamID = -1; - this->VListID = -1; - } + VGridFile.openURI(FileName); + int found = ReadVerticalGridData(); + if (!found) + { + VGridFile.openURI(GridFile.getURI()); + found = ReadVerticalGridData(); + } - char* directory = new char[strlen(this->FileName.c_str()) + 1]; - strcpy(directory, this->FileName.c_str()); - if (!this->OpenFile()) - { - return 0; - } + if (!found) + { + vtkErrorMacro("Could not get Vertical grid"); + return 0; + } - if (!this->ReadVerticalGridData()) - { - vtkDebugMacro("Couldn't neither open grid information within the data netCDF file, nor " - "in the grid.nc file."); - vtkErrorMacro("Couldn't neither open grid information within the data netCDF file, nor " - "in the grid.nc file."); - return 0; - } + this->FillGridDimensions(); - this->FileNameGrid = this->FileName; - if (!this->OpenFile()) + try + { + if (this->DimensionSelection >= 0) + { + if (DimensionSelection >= DimensionSets.size()) { + vtkErrorMacro("Trying to select inexistent dimensionset " + << DimensionSelection << " " << DimensionSets.size() << " are available."); return 0; } + for (int i = 0; i < Grids.size(); i++) + if (this->DimensionSets.at(this->DimensionSelection).GridSize == Grids.at(i).Size) + { + this->DimensionSets.at(this->DimensionSelection).GridID = Grids.at(i).GridID; + this->GridID = i; + } + this->ZAxisID = this->DimensionSets.at(this->DimensionSelection).ZAxisID; + vtkDebugMacro("NEW ZAxisID" << ZAxisID << " from " + << this->DimensionSets.at(this->DimensionSelection).ZAxisID); } + } + catch (const std::out_of_range& oor) + { + vtkErrorMacro("Out of Range error in GetDims trying to set Grid and ZAxisID: " << oor.what()); + return 0; + } - if (this->DimensionSelection > 0) - { - vlistNgrids(this->VListID); - int nzaxis = vlistNzaxis(this->VListID); - - this->GridID = vlistGrid(this->VListID, this->DimensionSelection / nzaxis); - this->ZAxisID = vlistZaxis( - this->VListID, this->DimensionSelection - (nzaxis * this->DimensionSelection / nzaxis)); - } - - if (this->GridID != -1) - { - this->NumberOfCells = static_cast(gridInqSize(this->GridID)); - - if (this->NumberOfPoints and - this->NumberOfPoints != static_cast(gridInqSize(this->GridID))) - vtkDebugMacro("GetDims: Changing number of points from " - << this->NumberOfPoints << " to " << static_cast(gridInqSize(this->GridID))); - this->NumberOfPoints = static_cast(gridInqSize(this->GridID)); - this->PointsPerCell = gridInqNvertex(this->GridID); - } - - int ntsteps = 0; - if (this->Grib) - { - while (streamInqTimestep(this->StreamID, ntsteps)) - ntsteps++; - } - else + try + { + if (GridID != -1 && Grids.at(this->GridID).GridID != -1) { - ntsteps = vlistNtsteps(this->VListID); - } - this->NumberOfTimeSteps = ntsteps; + this->NumberOfCells = static_cast(Grids.at(GridID).Size); - this->MaximumNVertLevels = 1; - if (this->ZAxisID != -1) - { - this->MaximumNVertLevels = zaxisInqSize(this->ZAxisID); + if (this->NumberOfPoints and this->NumberOfPoints != this->NumberOfCells) + vtkDebugMacro("GetDims: Changing number of points from " << this->NumberOfPoints << " to " + << this->NumberOfCells); + this->NumberOfPoints = this->NumberOfCells; + this->PointsPerCell = Grids.at(this->GridID).PointsPerCell; + vtkDebugMacro("GetDims: Found PointsPerCell to be " << this->PointsPerCell << " for grid " + << this->GridID); } + } + catch (const std::out_of_range& oor) + { + vtkErrorMacro("Out of Range error in GetDims trying to set NumberOfPoints " << oor.what()); + vtkErrorMacro("Grids.size " << Grids.size() << "\t GridID " << GridID); + return 0; + } - this->FillGridDimensions(); + int ntsteps = 0; + if (this->Grib) + { + while (streamInqTimestep(this->DataFile.getStreamID(), ntsteps)) + ntsteps++; } else { - vtkDebugMacro("No Filename yet set!"); + ntsteps = vlistNtsteps(this->DataFile.getVListID()); + } + this->NumberOfTimeSteps = ntsteps; + + this->MaximumNVertLevels = 1; + if (this->ZAxisID != -1) + { + this->MaximumNVertLevels = zaxisInqSize(this->ZAxisID); } return 1; } -//---------------------------------------------------------------------------- +//--------------------------------------------------------------------------------------------------- // Read Horizontal Grid Data -//---------------------------------------------------------------------------- +// Checks if there is at least one grid with >= 3 vertices, and if yes, sets GridID to this grid's +// ID +//--------------------------------------------------------------------------------------------------- int vtkCDIReader::ReadHorizontalGridData() { - int vlistID_l = this->VListID; - this->GridID = -1; - this->ZAxisID = -1; - this->SurfID = -1; - + Grids.resize(0); + int vlistID_l = this->GridFile.getVListID(); int ngrids = vlistNgrids(vlistID_l); for (int i = 0; i < ngrids; ++i) { @@ -983,17 +934,13 @@ int vtkCDIReader::ReadHorizontalGridData() if (nv >= 3) // ((nv == 3 || nv == 4)) // && gridInqType(gridID_l) == GRID_UNSTRUCTURED) { - this->GridID = gridID_l; - break; + Grid grid{ .GridID = gridID_l, .Size = gridInqSize(gridID_l), .PointsPerCell = nv }; + Grids.push_back(grid); } } - if (this->GridID == -1) - { - this->NeedHorizontalGridFile = true; + if (Grids.size() == 0) return 0; - } - return 1; } @@ -1003,37 +950,30 @@ int vtkCDIReader::ReadHorizontalGridData() int vtkCDIReader::ReadVerticalGridData() { this->ZAxisID = -1; - this->SurfID = -1; - int nzaxis = vlistNzaxis(this->VListID); - + int nzaxis = vlistNzaxis(this->VGridFile.getVListID()); + int found = 0; for (int i = 0; i < nzaxis; ++i) { - int zaxisID_l = vlistZaxis(this->VListID, i); + int zaxisID_l = vlistZaxis(this->VGridFile.getVListID(), i); if (zaxisInqSize(zaxisID_l) == 1 || zaxisInqType(zaxisID_l) == ZAXIS_SURFACE) { - this->SurfID = zaxisID_l; - this->ZAxisID = zaxisID_l; - break; + this->SurfIDs.insert(zaxisID_l); + + found = 1; } } for (int i = 0; i < nzaxis; ++i) { - int zaxisID_l = vlistZaxis(this->VListID, i); + int zaxisID_l = vlistZaxis(this->VGridFile.getVListID(), i); if (zaxisInqSize(zaxisID_l) > 1) { - this->ZAxisID = zaxisID_l; + found = 1; break; } } - if (this->ZAxisID == -1) - { - this->NeedVerticalGridFile = true; - return 0; - } - - return 1; + return found; } //---------------------------------------------------------------------------- @@ -1044,37 +984,42 @@ int vtkCDIReader::GetVars() int cellVarIndex = -1; int pointVarIndex = -1; int domainVarIndex = -1; + int numVars = vlistNvars(this->DataFile.getVListID()); + + vtkDebugMacro("Found " << numVars << " as Variables for VListID " << this->DataFile.getVListID()); - int numVars = vlistNvars(this->VListID); for (int i = 0; i < numVars; i++) { int varID = i; cdi_tools::CDIVar aVar; - aVar.StreamID = this->StreamID; + aVar.StreamID = this->DataFile.getStreamID(); aVar.VarID = varID; - aVar.GridID = vlistInqVarGrid(this->VListID, varID); - aVar.ZAxisID = vlistInqVarZaxis(this->VListID, varID); + aVar.GridID = vlistInqVarGrid(this->DataFile.getVListID(), varID); + aVar.ZAxisID = vlistInqVarZaxis(this->DataFile.getVListID(), varID); aVar.GridSize = static_cast(gridInqSize(aVar.GridID)); aVar.NLevel = zaxisInqSize(aVar.ZAxisID); aVar.Type = 0; aVar.ConstTime = 0; + vlistInqVarName(this->DataFile.getVListID(), varID, aVar.Name); + vtkDebugMacro("Processing variable " << i << '\t' << aVar.Name); // to do multiple grids: // - Check how many grids are available // - Check if all grids can be reconstructed, or if bnds are all zero // - Reform gui to load either Cell, Point or Edge data - if (vlistInqVarTsteptype(this->VListID, varID) == TIME_CONSTANT) + if (vlistInqVarTsteptype(this->DataFile.getVListID(), varID) == TIME_CONSTANT) { aVar.ConstTime = 1; } - if (aVar.ZAxisID != this->ZAxisID && aVar.ZAxisID != this->SurfID) + if (aVar.ZAxisID != this->ZAxisID && SurfIDs.count(aVar.ZAxisID) == 0) + // We are handling a different 3D Axis. { + vtkDebugMacro("Skipping " << aVar.Name << " as it has the wrong ZAxis " << aVar.ZAxisID); continue; } - vlistInqVarName(this->VListID, varID, aVar.Name); aVar.Type = 2; if (aVar.NLevel > 1) { @@ -1091,19 +1036,23 @@ int vtkCDIReader::GetVars() } else if ((aVar.GridSize < this->NumberOfCells) && (this->PointsPerCell == 3)) { + vtkDebugMacro("Skipping " << aVar.Name << " as it has the wrong GridSize " << aVar.GridSize + << " != " << this->NumberOfCells); if (this->NumberOfPoints and this->NumberOfPoints != aVar.GridSize) { vtkWarningMacro("Not adding " << aVar.Name << " as point var, as it's size " << aVar.GridSize << " does not correspond to our understanding of the correct size for 'the' point grid: " << this->NumberOfPoints); - break; + continue; } isPointData = true; this->NumberOfPoints = aVar.GridSize; } else { + vtkDebugMacro("Skipping " << aVar.Name << " as it has the wrong GridSize " << aVar.GridSize + << " != " << this->NumberOfCells); continue; } @@ -1219,7 +1168,7 @@ int vtkCDIReader::BuildVarArrays() << " NumberOfPointVars: " << this->NumberOfPointVars); if (this->NumberOfCellVars == 0) { - vtkErrorMacro("No cell variables found!"); + vtkDebugMacro("No cell variables found!"); } for (int var = 0; var < this->NumberOfPointVars; var++) @@ -1422,39 +1371,63 @@ int vtkCDIReader::ConstructGridGeometry() CHECK_NEW(this->DepthVar); vtkDebugMacro("Start reading Vertices"); - gridInqXboundsPart( - this->GridID, (this->BeginCell * this->PointsPerCell), size, cLonVertices.data()); - gridInqYboundsPart( - this->GridID, (this->BeginCell * this->PointsPerCell), size, cLatVertices.data()); + try + { + gridInqXboundsPart(Grids.at(this->GridID).GridID, (this->BeginCell * this->PointsPerCell), size, + cLonVertices.data()); + gridInqYboundsPart(Grids.at(this->GridID).GridID, (this->BeginCell * this->PointsPerCell), size, + cLatVertices.data()); + } + catch (const std::out_of_range& oor) + { + vtkErrorMacro( + "Out of Range error trying to get the grid id for reading vertices: " << oor.what()); + return 0; + } vtkDebugMacro("Done reading Vertices"); + vtkDebugMacro("Getting vertical axis" << this->ZAxisID << " expecting up to " + << this->MaximumNVertLevels << " levels."); zaxisInqLevels(this->ZAxisID, this->DepthVar); + vtkDebugMacro("Got vertical axis" << this->ZAxisID); char units[CDI_MAX_NAME]; this->OrigConnections.resize(size); int new_cells[2]; - - if (this->ProjectionMode != projection::CATALYST) + try { - gridInqXunits(this->GridID, units); - if (strncmp(units, "degree", 6) == 0) + if (this->ProjectionMode != projection::CATALYST) { - for (int i = 0; i < size; i++) + gridInqXunits(Grids.at(this->GridID).GridID, units); + if (strncmp(units, "degree", 6) == 0) { - cLonVertices[i] = vtkMath::RadiansFromDegrees(cLonVertices[i]); + for (int i = 0; i < size; i++) + { + cLonVertices[i] = vtkMath::RadiansFromDegrees(cLonVertices[i]); + } } - } - gridInqYunits(this->GridID, units); - if (strncmp(units, "degree", 6) == 0) - { - for (int i = 0; i < size; i++) + gridInqYunits(Grids.at(this->GridID).GridID, units); + if (strncmp(units, "degree", 6) == 0) { - cLatVertices[i] = vtkMath::RadiansFromDegrees(cLatVertices[i]); + for (int i = 0; i < size; i++) + { + cLatVertices[i] = vtkMath::RadiansFromDegrees(cLatVertices[i]); + } } } } + catch (const std::out_of_range& oor) + { + vtkErrorMacro("Out of Range error trying to get the grid id for getting the coordinate units " + "for projecting: " + << oor.what()); + return 0; + } // check for duplicates in the Point list and update the triangle list + vtkDebugMacro("Removing duplicates for clon/clat, size = " << size); + this->RemoveDuplicates( cLonVertices.data(), cLatVertices.data(), size, &this->OrigConnections[0], new_cells); + vtkDebugMacro("Removed duplicates for clon/clat"); this->NumberLocalCells = new_cells[0] / this->PointsPerCell; this->NumberLocalPoints = new_cells[1]; if (this->NumberOfPoints and this->NumberOfPoints != new_cells[1]) @@ -1503,31 +1476,43 @@ int vtkCDIReader::ConstructGridGeometry() if (this->Piece == 0) { int new_cells2[2]; - double clon_vert2[size2]; - double clat_vert2[size2]; - - gridInqXboundsPart(this->GridID, 0, size2, clon_vert2); - gridInqYboundsPart(this->GridID, 0, size2, clat_vert2); - - gridInqXunits(this->GridID, units); - if (strncmp(units, "degree", 6) == 0) + std::vector clon_vert2(size2); + std::vector clat_vert2(size2); + try { - for (int i = 0; i < size2; i++) + gridInqXboundsPart(Grids.at(this->GridID).GridID, 0, size2, clon_vert2.data()); + gridInqYboundsPart(Grids.at(this->GridID).GridID, 0, size2, clat_vert2.data()); + + gridInqXunits(Grids.at(this->GridID).GridID, units); + if (strncmp(units, "degree", 6) == 0) { - clon_vert2[i] = vtkMath::RadiansFromDegrees(clon_vert2[i]); + for (int i = 0; i < size2; i++) + { + clon_vert2[i] = vtkMath::RadiansFromDegrees(clon_vert2[i]); + } } - } - gridInqYunits(this->GridID, units); - if (strncmp(units, "degree", 6) == 0) - { - for (int i = 0; i < size2; i++) + gridInqYunits(Grids.at(this->GridID).GridID, units); + if (strncmp(units, "degree", 6) == 0) { - clat_vert2[i] = vtkMath::RadiansFromDegrees(clat_vert2[i]); + for (int i = 0; i < size2; i++) + { + clat_vert2[i] = vtkMath::RadiansFromDegrees(clat_vert2[i]); + } } } + catch (const std::out_of_range& oor) + { + vtkErrorMacro( + "Out of Range error trying to get the grid id for converting lat/lon in parallel: " + << oor.what()); + return 0; + } - this->RemoveDuplicates(clon_vert2, clat_vert2, size2, vertex_ids2.data(), new_cells2); + vtkDebugMacro("Removing duplicates for clon/clat2"); + + this->RemoveDuplicates( + clon_vert2.data(), clat_vert2.data(), size2, vertex_ids2.data(), new_cells2); for (int i = 1; i < this->NumPieces; i++) { this->Controller->Send(vertex_ids2.data(), size2, i, 101); @@ -1655,26 +1640,39 @@ int vtkCDIReader::LoadClonClatVars() std::vector cLon_l(this->NumberLocalCells); std::vector cLat_l(this->NumberLocalCells); - gridInqXvalsPart(this->GridID, this->BeginCell, this->NumberLocalCells, cLon_l.data()); - gridInqYvalsPart(this->GridID, this->BeginCell, this->NumberLocalCells, cLat_l.data()); + gridInqXvalsPart( + Grids.at(this->GridID).GridID, this->BeginCell, this->NumberLocalCells, cLon_l.data()); + gridInqYvalsPart( + Grids.at(this->GridID).GridID, this->BeginCell, this->NumberLocalCells, cLat_l.data()); char units[CDI_MAX_NAME]; - gridInqXunits(this->GridID, units); - if (strncmp(units, "degree", 6) == 0) + + try { - for (int i = 0; i < this->NumberLocalCells; i++) + gridInqXunits(Grids.at(this->GridID).GridID, units); + if (strncmp(units, "degree", 6) == 0) { - cLon_l[i] = vtkMath::RadiansFromDegrees(cLon_l[i]); + for (int i = 0; i < this->NumberLocalCells; i++) + { + cLon_l[i] = vtkMath::RadiansFromDegrees(cLon_l[i]); + } } - } - gridInqYunits(this->GridID, units); - if (strncmp(units, "degree", 6) == 0) - { - for (int i = 0; i < this->NumberLocalCells; i++) + gridInqYunits(Grids.at(this->GridID).GridID, units); + if (strncmp(units, "degree", 6) == 0) { - cLat_l[i] = vtkMath::RadiansFromDegrees(cLat_l[i]); + for (int i = 0; i < this->NumberLocalCells; i++) + { + cLat_l[i] = vtkMath::RadiansFromDegrees(cLat_l[i]); + } } } + catch (const std::out_of_range& oor) + { + vtkErrorMacro( + "Out of Range error trying to get the grid id for converting lat/lon in LoadClonClatVars: " + << oor.what()); + return 0; + } if (this->ShowMultilayerView) { @@ -1802,7 +1800,7 @@ int vtkCDIReader::CheckForMaskData() { const double maskVal = this->UseCustomMaskValue ? this->CustomMaskValue - : vlistInqVarMissval(this->VListID, this->Internals->CellVars[mask_pos].VarID); + : vlistInqVarMissval(this->DataFile.getVListID(), this->Internals->CellVars[mask_pos].VarID); cdi_tools::CDIVar* cdiVar = &(this->Internals->CellVars[mask_pos]); if (this->ShowMultilayerView) @@ -1902,7 +1900,7 @@ bool vtkCDIReader::BuildDomainCellVars() CHECK_NEW(this->DomainCellVar); double val = 0; int mask_pos = 0; - int numVars = vlistNvars(this->VListID); + int numVars = vlistNvars(this->DataFile.getVListID()); for (int i = 0; i < numVars; i++) { @@ -2759,7 +2757,7 @@ int vtkCDIReader::LoadCellVarDataTemplate( //------------------------------------------------------------------------------ int vtkCDIReader::ReplaceFillWithNan(const int varID, vtkDataArray* dataArray) { - double miss = vlistInqVarMissval(this->VListID, varID); + double miss = vlistInqVarMissval(this->DataFile.getVListID(), varID); // NaN only available with float and double. if (dataArray->GetDataType() == VTK_FLOAT) @@ -3072,11 +3070,11 @@ int vtkCDIReader::LoadDomainVarData(int variableIndex) //----------------------------------------------------------------------------- int vtkCDIReader::FillGridDimensions() { - int ngrids = vlistNgrids(this->VListID); - int nzaxis = vlistNzaxis(this->VListID); - int nvars = vlistNvars(this->VListID); - this->AllDimensions->SetNumberOfValues(0); - this->VariableDimensions->SetNumberOfValues(ngrids * nzaxis); + this->DimensionSets.resize(0); + + int ngrids = vlistNgrids(this->DataFile.getVListID()); + int nzaxis = vlistNzaxis(this->DataFile.getVListID()); + int nvars = vlistNvars(this->DataFile.getVListID()); char nameGridX[CDI_MAX_NAME]; char nameGridY[CDI_MAX_NAME]; char nameLev[CDI_MAX_NAME]; @@ -3085,40 +3083,59 @@ int vtkCDIReader::FillGridDimensions() for (int k = 0; k < nvars; k++) { - int i = vlistInqVarGrid(this->VListID, k); - int j = vlistInqVarZaxis(this->VListID, k); + int i = vlistInqVarGrid(this->DataFile.getVListID(), k); + int j = vlistInqVarZaxis(this->DataFile.getVListID(), k); hits.insert(std::to_string(i) + "x" + std::to_string(j)); // IDs are not 0 to n-1 but can be 30-ish for a file with 3 grids. // they map to the gridID_l and zaxisID_l values below. // Thus we need to a map to catch rather unpredictable values. } - + size_t counter = 0; for (int i = 0; i < ngrids; ++i) { for (int j = 0; j < nzaxis; ++j) { std::string dimEncoding("("); - int gridID_l = vlistGrid(this->VListID, i); + int gridID_l = vlistGrid(this->DataFile.getVListID(), i); gridInqXname(gridID_l, nameGridX); gridInqYname(gridID_l, nameGridY); dimEncoding += nameGridX; dimEncoding += ", "; dimEncoding += nameGridY; dimEncoding += ", "; - int zaxisID_l = vlistZaxis(this->VListID, j); + int zaxisID_l = vlistZaxis(this->DataFile.getVListID(), j); zaxisInqName(zaxisID_l, nameLev); dimEncoding += nameLev; dimEncoding += ")"; if (hits.count(std::to_string(gridID_l) + "x" + std::to_string(zaxisID_l)) == 0) { + vtkDebugMacro("vtkCDIReader::FillGridDimensions: i, j, dimEncoding: " + << i << '\t' << j << "\t" << gridID_l << '\t' << zaxisID_l << "\t" << dimEncoding + << " - has no hits.\n"); continue; // skip empty grid combinations } + vtkDebugMacro("vtkCDIReader::FillGridDimensions: i, j, GridID, ZAxisID, dimEncoding: " + << i << '\t' << j << "\t" << gridID_l << '\t' << zaxisID_l << "\t" << dimEncoding + << " - has hits.\n"); - this->AllDimensions->InsertNextValue(dimEncoding); - this->VariableDimensions->SetValue(i * nzaxis + j, dimEncoding.c_str()); + dimset ds{ .DimsetID = counter, + .GridID = -1, + .ZAxisID = zaxisID_l, + .GridSize = gridInqSize(gridID_l), + .NLevel = zaxisInqSize(zaxisID_l), + .label = dimEncoding }; + DimensionSets.push_back(ds); + counter++; } } + this->AllDimensions->SetNumberOfValues(0); + this->VariableDimensions->SetNumberOfValues(counter); + for (int i = 0; i < counter; i++) + { + this->AllDimensions->InsertNextValue(DimensionSets[i].label); + this->VariableDimensions->SetValue(i, DimensionSets[i].label.c_str()); + } return 1; } @@ -3298,18 +3315,13 @@ const char* vtkCDIReader::GetDomainArrayName(int index) } //---------------------------------------------------------------------------- -// Set to lat/lon (equidistant cylindrical) projection. +// Load a new file. //---------------------------------------------------------------------------- void vtkCDIReader::SetFileName(const char* val) { if (this->FileName.empty() || val == nullptr || strcmp(this->FileName.c_str(), val) != 0) { - if (this->StreamID >= 0) - { - streamClose(this->StreamID); - this->StreamID = -1; - this->VListID = -1; - } + this->DataFile.setVoid(); this->Modified(); if (val == nullptr) { @@ -3330,13 +3342,21 @@ void vtkCDIReader::SetVerticalLevel(int level) { if (this->VerticalLevelSelected != level) { - if (level < 0 || level > this->MaximumNVertLevels - 1) + if (level < 0) { vtkErrorMacro("Requested inexistent vertical level: " << level << ".\nThe level must be the in range [ 0 ; " << this->MaximumNVertLevels - 1 << " ]."); return; } + + if (level > this->MaximumNVertLevels - 1) + { + vtkWarningMacro("Requested inexistent vertical level: " + << level << ".\nThe level must be the in range [ 0 ; " << this->MaximumNVertLevels - 1 + << " ]. \nTriying with 0."); + level = 0; + } this->VerticalLevelSelected = level; this->Modified(); vtkDebugMacro("Set VerticalLevelSelected to: " << level); diff --git a/Plugins/CDIReader/Reader/vtkCDIReader.h b/Plugins/CDIReader/Reader/vtkCDIReader.h index 82bfac79455523bc4f13609d8b33f6f84b12fdad..7e2b27fd26607c5f2e88e849a89d89b8b937851e 100644 --- a/Plugins/CDIReader/Reader/vtkCDIReader.h +++ b/Plugins/CDIReader/Reader/vtkCDIReader.h @@ -50,7 +50,9 @@ #include "projections.h" // for projection enum +#include "DataSource.h" #include // for unique_ptr +#include #include // for std::vector class vtkCallbackCommand; @@ -58,6 +60,23 @@ class vtkDoubleArray; class vtkFieldData; class vtkMultiProcessController; +struct dimset +{ + size_t DimsetID; + int GridID; + int ZAxisID; + size_t GridSize; + int NLevel; + std::string label; +}; + +struct Grid +{ + int GridID; + size_t Size; + int PointsPerCell; +}; + class VTKCDIREADER_EXPORT vtkCDIReader : public vtkUnstructuredGridAlgorithm { public: @@ -305,6 +324,8 @@ protected: double Layer0OffsetRange[2]; int DimensionSelection; + std::vector DimensionSets; + std::vector Grids; bool InvertZAxis; bool AddCoordinateVars; projection::Projection ProjectionMode; @@ -349,11 +370,10 @@ protected: int NumberOfDomainVars; bool GridReconstructed; - int StreamID; - int VListID; + DataSource::CDIObject DataFile, GridFile, VGridFile; int GridID; int ZAxisID; - int SurfID; + std::unordered_set SurfIDs; std::string TimeUnits; std::string Calendar;