diff --git a/README.md b/README.md index 6729f10..69848c2 100644 --- a/README.md +++ b/README.md @@ -27,7 +27,7 @@ To install from conda use the command line: # Installation via CMake (Recommended) -TempestExtremes can be built and installed on various systems using CMake. Our new script, `./quick_make_unix.sh`, automatically detects your platform(UNIX-based systems only) and loads any required modules before building. +TempestExtremes can be built and installed on various systems using CMake. Our new script, `./quick_make_unix.sh` in the TempestExtremes root directory, automatically detects your operating system (Unix-based systems only) and loads any required modules before building. To install all TempestExtremes executables to the ./bin directory, just run the quick make script. ## General CMake Configuration diff --git a/src/base/DataOp.cpp b/src/base/DataOp.cpp index 4a0123b..cdc340f 100644 --- a/src/base/DataOp.cpp +++ b/src/base/DataOp.cpp @@ -436,9 +436,14 @@ bool DataOp_VECMAG::Apply( const DataArray1D & dataRight = *(vecArgData[1]); for (int i = 0; i < dataout.GetRows(); i++) { - dataout[i] = - sqrt(dataLeft[i] * dataLeft[i] - + dataRight[i] * dataRight[i]); + bool fIsFillValue = dataLeft.IsFillValueAtIx(i); + if (!fIsFillValue){ + dataout[i] = + sqrt(dataLeft[i] * dataLeft[i] + + dataRight[i] * dataRight[i]); + } else { + dataout[i]=std::numeric_limits::quiet_NaN(); + } } return true; @@ -699,7 +704,12 @@ bool DataOp_DIFF::Apply( const DataArray1D & dataRight = *(vecArgData[1]); for (int i = 0; i < dataout.GetRows(); i++) { - dataout[i] = dataLeft[i] - dataRight[i]; + bool fIsFillValue = (dataLeft.IsFillValueAtIx(i) || dataRight.IsFillValueAtIx(i)); + if (!fIsFillValue){ + dataout[i] = dataLeft[i] - dataRight[i]; + } else { + dataout[i]=std::numeric_limits::quiet_NaN(); + } } } @@ -1702,7 +1712,6 @@ void BuildCurlOperator( ) { opCurlE.Clear(); opCurlN.Clear(); - int iRef = 0; // Scaling factor used in Curl calculation @@ -1955,7 +1964,6 @@ bool DataOp_CURL::Apply( m_fInitialized = true; } - m_opCurlE.Apply(*(vecArgData[0]), dataout, true); m_opCurlN.Apply(*(vecArgData[1]), dataout, false); diff --git a/src/base/SparseMatrix.h b/src/base/SparseMatrix.h index e8bbc31..9165260 100644 --- a/src/base/SparseMatrix.h +++ b/src/base/SparseMatrix.h @@ -18,7 +18,6 @@ #define _SPARSEMATRIX_H_ #include "DataArray1D.h" - #include /////////////////////////////////////////////////////////////////////////////// @@ -190,8 +189,13 @@ class SparseMatrix { SparseMapConstIterator iter = m_mapEntries.begin(); for (; iter != m_mapEntries.end(); iter++) { - dataVectorOut[iter->first.first] += + if (dataVectorIn[iter->first.second]>1e14){ + // dataVectorOut[iter->first.first] += std::numeric_limits::quiet_NaN(); + continue; + } else { + dataVectorOut[iter->first.first] += iter->second * dataVectorIn[iter->first.second]; + } } } diff --git a/src/nodes/NodeFileFilter.cpp b/src/nodes/NodeFileFilter.cpp index b015a99..7b5cf64 100644 --- a/src/nodes/NodeFileFilter.cpp +++ b/src/nodes/NodeFileFilter.cpp @@ -69,7 +69,7 @@ class NearbyBlobsOp { m_varix(InvalidVariableIndex), m_dDistanceDeg(0.0), m_eOp(GreaterThan), - m_dValue(0.0), + m_strValue(""), m_dMaxDistDeg(180.0) { } @@ -133,7 +133,7 @@ class NearbyBlobsOp { // Read in value } else if (eReadMode == ReadMode_Value) { - m_dValue = atof(strSubStr.c_str()); + m_strValue = strSubStr; iLast = i + 1; eReadMode = ReadMode_MaxDist; @@ -191,12 +191,7 @@ class NearbyBlobsOp { strDescription += " is not equal to "; } - if (fabs(m_dValue) < 1.0e-4) { - snprintf(szBuffer, 128, "%e", m_dValue); - } else { - snprintf(szBuffer, 128, "%f", m_dValue); - } - strDescription += szBuffer; + strDescription += m_strValue; snprintf(szBuffer, 128, "%f", m_dMaxDistDeg); strDescription += std::string(" (max dist ") + szBuffer + " degrees)"; @@ -209,35 +204,36 @@ class NearbyBlobsOp { /// Determine if a particular value satisfies this threshold. /// bool SatisfiedBy( - double dValue + double dValue, + double dOpValue ) const { if (m_eOp == GreaterThan) { - if (dValue > m_dValue) { + if (dValue > dOpValue) { return true; } } else if (m_eOp == LessThan) { - if (dValue < m_dValue) { + if (dValue < dOpValue) { return true; } } else if (m_eOp == GreaterThanEqualTo) { - if (dValue >= m_dValue) { + if (dValue >= dOpValue) { return true; } } else if (m_eOp == LessThanEqualTo) { - if (dValue <= m_dValue) { + if (dValue <= dOpValue) { return true; } } else if (m_eOp == EqualTo) { - if (dValue == m_dValue) { + if (dValue == dOpValue) { return true; } } else if (m_eOp == NotEqualTo) { - if (dValue != m_dValue) { + if (dValue != dOpValue) { return true; } @@ -267,7 +263,7 @@ class NearbyBlobsOp { /// /// Threshold value. /// - double m_dValue; + std::string m_strValue; /// /// Maximum distance in degrees. @@ -441,6 +437,7 @@ void BuildMask_ByContour( template void BuildMask_NearbyBlobs( const SimpleGrid & grid, + const NodeFile & nodefile, const DataArray1D & dataState, std::vector & vecLonRad, std::vector & vecLatRad, @@ -450,6 +447,16 @@ void BuildMask_NearbyBlobs( // Get the variable const double dDistDeg = nearbyblobsop.m_dDistanceDeg; const double dMaxDistDeg = nearbyblobsop.m_dMaxDistDeg; + std::string strValue = nearbyblobsop.m_strValue; + double dOpValue = 0.0; + + std::vector VecValue; + if (STLStringHelper::IsFloat(strValue)) { + dOpValue = std::stod(strValue); + } else { + nodefile.InterpolatedColumnDouble( + strValue, VecValue); + } _ASSERT(vecLonRad.size() == vecLatRad.size()); _ASSERT(dataState.GetRows() == grid.GetSize()); @@ -472,7 +479,7 @@ void BuildMask_NearbyBlobs( std::queue queueNodes; queueNodes.push(ix0); - // Set of nodes that have already beenvisited + // Set of nodes that have already been visited std::set setNodesVisited; // Loop through all elements @@ -504,8 +511,12 @@ void BuildMask_NearbyBlobs( queueNodes.push(grid.m_vecConnectivity[ix][n]); } + if (!VecValue.empty()){ + dOpValue=VecValue[j]; + } + // Check if this point satisfies the nearbyblobs criteria - if (!nearbyblobsop.SatisfiedBy(static_cast(dataState[ix]))) { + if (!nearbyblobsop.SatisfiedBy(static_cast(dataState[ix]),dOpValue)) { continue; } @@ -543,7 +554,7 @@ void BuildMask_NearbyBlobs( } // Verify this point satisfies the condition - if (!nearbyblobsop.SatisfiedBy(static_cast(dataState[ixblob]))) { + if (!nearbyblobsop.SatisfiedBy(static_cast(dataState[ixblob]),dOpValue)) { // Isn't part of the blob, but add it to the list of // nodes to visit. @@ -562,9 +573,9 @@ void BuildMask_NearbyBlobs( queueThresholdedNodes.push(grid.m_vecConnectivity[ixblob][n]); } - dataMask[ixblob] = 1.0; + dataMask[ixblob] += 1.0; } - } + } } } @@ -976,23 +987,36 @@ void NodeFileFilter( AnnounceEndBlock("Done"); } + int tc = 0; + if (vecNearbyBlobsOp.size() != 0) { AnnounceStartBlock("Building mask (nearbyblobs)"); - Variable & varOp = varreg.Get(vecNearbyBlobsOp[0].m_varix); - vecInFiles.SetConstantTimeIx(t); - varOp.LoadGridData(varreg, vecInFiles, grid); - const DataArray1D & dataState = varOp.GetData(); - _ASSERT(dataState.GetRows() == grid.GetSize()); - - BuildMask_NearbyBlobs( - grid, - dataState, - vecLonRad, - vecLatRad, - vecNearbyBlobsOp[0], - dataMask); - - AnnounceEndBlock("Done"); + for (tc = 0; tc < vecNearbyBlobsOp.size(); tc++) { + Variable & varOp = varreg.Get(vecNearbyBlobsOp[tc].m_varix); + vecInFiles.SetConstantTimeIx(t); + varOp.LoadGridData(varreg, vecInFiles, grid); + const DataArray1D & dataState = varOp.GetData(); + _ASSERT(dataState.GetRows() == grid.GetSize()); + int vecsize = vecNearbyBlobsOp.size(); + BuildMask_NearbyBlobs( + grid, + nodefile, + dataState, + vecLonRad, + vecLatRad, + vecNearbyBlobsOp[tc], + dataMask); + + AnnounceEndBlock("Done"); + } + } + + if (tc > 1){ + for (int i = 0; i < dataMask.GetRows(); i++) { + if(dataMask[i] > 0){ + dataMask[i] = dataMask[i] - (tc-1); + } + } } if (fInvert) { @@ -1381,9 +1405,23 @@ try { if (strNearbyBlobs != "") { AnnounceStartBlock("Parsing --nearbyblobs"); - NearbyBlobsOp op; - op.Parse(varreg, strNearbyBlobs); - vecNearbyBlobsOp.push_back(op); + int iLast = 0; + for (int i = 0; i <= strNearbyBlobs.length(); i++) { + + if ((i == strNearbyBlobs.length()) || + (strNearbyBlobs[i] == ';') || + (strNearbyBlobs[i] == ':') + ) { + std::string strSubStr = + strNearbyBlobs.substr(iLast, i - iLast); + + int iNextOp = (int)(vecNearbyBlobsOp.size()); + vecNearbyBlobsOp.resize(iNextOp + 1); + vecNearbyBlobsOp[iNextOp].Parse(varreg, strSubStr); + + iLast = i + 1; + } + } AnnounceEndBlock(NULL); }