diff --git a/.travis.yml b/.travis.yml index 45d829e9e..7afe8f8a7 100644 --- a/.travis.yml +++ b/.travis.yml @@ -19,12 +19,13 @@ env: - BUILD_TYPE=Debug - TECA_DIR=/travis_teca_dir - TECA_PYTHON_VERSION=3 - - TECA_DATA_REVISION=163 + - TECA_DATA_REVISION=164 + - TECA_TELITE_REVISION=188dc11d jobs: - DOCKER_IMAGE=ubuntu IMAGE_VERSION=22.04 IMAGE_NAME=ubuntu_22_04 REQUIRE_NETCDF_MPI=TRUE - - DOCKER_IMAGE=ubuntu IMAGE_VERSION=22.04 IMAGE_NAME=ubuntu_22_04 REQUIRE_NETCDF_MPI=FALSE + #- DOCKER_IMAGE=ubuntu IMAGE_VERSION=22.04 IMAGE_NAME=ubuntu_22_04 REQUIRE_NETCDF_MPI=FALSE - DOCKER_IMAGE=fedora IMAGE_VERSION=37 IMAGE_NAME=fedora_37 REQUIRE_NETCDF_MPI=TRUE - - DOCKER_IMAGE=fedora IMAGE_VERSION=37 IMAGE_NAME=fedora_37 REQUIRE_NETCDF_MPI=FALSE + #- DOCKER_IMAGE=fedora IMAGE_VERSION=37 IMAGE_NAME=fedora_37 REQUIRE_NETCDF_MPI=FALSE - NO_DOCKER=TRUE jobs: @@ -66,6 +67,7 @@ install: "export TECA_PYTHON_VERSION=${TECA_PYTHON_VERSION} && export TECA_DATA_REVISION=${TECA_DATA_REVISION} && export REQUIRE_NETCDF_MPI=${REQUIRE_NETCDF_MPI} && + export TECA_TELITE_REVISION=${TECA_TELITE_REVISION} && ${TECA_DIR}/test/travis_ci/install_${IMAGE_NAME}.sh"; fi @@ -83,5 +85,6 @@ script: export DOCKER_IMAGE=${DOCKER_IMAGE} && export IMAGE_VERSION=${IMAGE_VERSION} && export REQUIRE_NETCDF_MPI=${REQUIRE_NETCDF_MPI} && + export TECA_TELITE_REVISION=${TECA_TELITE_REVISION} && ${TECA_DIR}/test/travis_ci/ctest_linux.sh"; fi diff --git a/CMake/FindNetCDF.cmake b/CMake/FindNetCDF.cmake index 7f7a9957b..43757f05d 100644 --- a/CMake/FindNetCDF.cmake +++ b/CMake/FindNetCDF.cmake @@ -95,6 +95,7 @@ if (NOT NC_TMP_FOUND OR NOT NC_TMP_LINK_LIBRARIES OR NOT NC_TMP_LIBRARY_DIRS OR endif() # look for header file that indicates MPI support +set(NETCDF_INCLUDE_DIR ${NETCDF_INCLUDE_DIR} /usr/lib/x86_64-linux-gnu/netcdf/mpi/include) set(NETCDF_IS_PARALLEL FALSE) find_file(NETCDF_PAR_INCLUDE_DIR netcdf_par.h PATHS ${NETCDF_INCLUDE_DIR} NO_DEFAULT_PATH) diff --git a/CMakeLists.txt b/CMakeLists.txt index bb9d99809..ed35974ca 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -281,6 +281,19 @@ if (NOT TECA_HAS_PARAVIEW) endif() set(TECA_HAS_VTK ${tmp} CACHE BOOL "VTK features") +# configure for TELite +set(tmp OFF) +find_package(TELite QUIET) +if (TELite_FOUND AND ((DEFINED TECA_HAS_TELITE AND TECA_HAS_TELITE) OR (NOT DEFINED TECA_HAS_TELITE))) + message(STATUS "TELite features -- enabled") + set(tmp ON) +elseif (REQUIRE_TELITE) + message(FATAL_ERROR "TELite features -- required but not found. set TELite_DIR to enable.") +else() + message(STATUS "TELite features -- not found. set TELite_DIR to enable.") +endif() +set(TECA_HAS_TELITE ${tmp} CACHE BOOL "TELite features") + #configure for Boost set(tmp OFF) find_package(Boost QUIET COMPONENTS program_options) diff --git a/alg/CMakeLists.txt b/alg/CMakeLists.txt index 2b2932d97..8fe7657fb 100644 --- a/alg/CMakeLists.txt +++ b/alg/CMakeLists.txt @@ -18,6 +18,7 @@ set(teca_alg_cxx_srcs teca_cartesian_mesh_source.cxx teca_cartesian_mesh_subset.cxx teca_cartesian_mesh_regrid.cxx + teca_mesh_join.cxx teca_connected_components.cxx teca_2d_component_area.cxx teca_component_area_filter.cxx @@ -65,8 +66,11 @@ set(teca_alg_cxx_srcs teca_vorticity.cxx teca_dataset_diff.cxx teca_temporal_reduction.cxx + teca_detect_nodes.cxx + teca_stitch_nodes.cxx ) + set(teca_alg_cuda_srcs) if (TECA_HAS_CUDA) set(teca_alg_cuda_srcs @@ -107,6 +111,10 @@ if (TECA_HAS_BOOST) list(APPEND teca_alg_link ${Boost_LIBRARIES}) endif() +if (TECA_HAS_TELITE) + list(APPEND teca_alg_link TELite) +endif() + if (TECA_HAS_CUDA) set_source_files_properties(${teca_alg_cxx_srcs} PROPERTIES LANGUAGE CUDA) endif() diff --git a/alg/cuCompactor.cuh b/alg/cuCompactor.cuh new file mode 100644 index 000000000..2d2ea0cef --- /dev/null +++ b/alg/cuCompactor.cuh @@ -0,0 +1,165 @@ +/* + * cuCompactor.h + * + * Created on: 21/mag/2015 + * Author: knotman + */ + +#ifndef CUCOMPACTOR_H_ +#define CUCOMPACTOR_H_ + +#include +#include + +#include +#include + +// Define this to turn on error checking +#define CUDA_ERROR_CHECK + +#define CUDASAFECALL( err ) __cudaSafeCall( err, __FILE__, __LINE__ ) + +inline void __cudaSafeCall( cudaError err, const char *file, const int line ) +{ +#ifdef CUDA_ERROR_CHECK + if ( cudaSuccess != err ) + { + fprintf( stderr, "cudaSafeCall() failed at %s:%i : %s\n", + file, line, cudaGetErrorString( err ) ); + + fprintf( stdout, "cudaSafeCall() failed at %s:%i : %s\n", + file, line, cudaGetErrorString( err ) ); + exit( -1 ); + } +#endif + + return; +} + +namespace cuCompactor { + +#define warpSize (32) +#define FULL_MASK 0xffffffff + +__host__ __device__ int divup(int x, int y) { return x / y + (x % y ? 1 : 0); } + +__device__ __inline__ int pow2i (int e){ + return 1< +__global__ void computeBlockCounts(T* d_input,int length,int*d_BlockCounts,Predicate predicate){ + int idx = threadIdx.x + blockIdx.x*blockDim.x; + if(idx < length){ + int pred = predicate(d_input[idx]); + int BC=__syncthreads_count(pred); + + if(threadIdx.x==0){ + d_BlockCounts[blockIdx.x]=BC; // BC will contain the number of valid elements in all threads of this thread block + } + } +} + + + +template +__global__ void compactK(T* d_input,int length, T* d_output,int* d_BlocksOffset,Predicate predicate ){ + int idx = threadIdx.x + blockIdx.x*blockDim.x; + extern __shared__ int warpTotals[]; + if(idx < length){ + int pred = predicate(d_input[idx]); + int w_i = threadIdx.x/warpSize; //warp index + int w_l = idx % warpSize;//thread index within a warp + + // compute exclusive prefix sum based on predicate validity to get output offset for thread in warp + int t_m = FULL_MASK >> (warpSize-w_l); //thread mask + #if (CUDART_VERSION < 9000) + int b = __ballot(pred) & t_m; //ballot result = number whose ith bit is one if the ith's thread pred is true masked up to the current index in warp + #else + int b = __ballot_sync(FULL_MASK,pred) & t_m; + #endif + int t_u = __popc(b); // popc count the number of bit one. simply count the number predicated true BEFORE MY INDEX + + // last thread in warp computes total valid counts for the warp + if(w_l==warpSize-1){ + warpTotals[w_i]=t_u+pred; + } + + // need all warps in thread block to fill in warpTotals before proceeding + __syncthreads(); + + // first numWarps threads in first warp compute exclusive prefix sum to get output offset for each warp in thread block + int numWarps = blockDim.x/warpSize; + unsigned int numWarpsMask = FULL_MASK >> (warpSize-numWarps); + if(w_i==0 && w_l +__global__ void printArray_GPU(T* hd_data, int size,int newline){ + int w=0; + for(int i=0;i ",w); + w++; + } + printf("%i ",hd_data[i]); + } + printf("\n"); +} + +template +int compact(T* d_input,T* d_output,int length, Predicate predicate, int blockSize){ + int numBlocks = divup(length,blockSize); + int* d_BlocksCount; + int* d_BlocksOffset; + CUDASAFECALL (cudaMalloc(&d_BlocksCount,sizeof(int)*numBlocks)); + CUDASAFECALL (cudaMalloc(&d_BlocksOffset,sizeof(int)*numBlocks)); + thrust::device_ptr thrustPrt_bCount(d_BlocksCount); + thrust::device_ptr thrustPrt_bOffset(d_BlocksOffset); + + //phase 1: count number of valid elements in each thread block + computeBlockCounts<<>>(d_input,length,d_BlocksCount,predicate); + + //phase 2: compute exclusive prefix sum of valid block counts to get output offset for each thread block in grid + thrust::exclusive_scan(thrustPrt_bCount, thrustPrt_bCount + numBlocks, thrustPrt_bOffset); + + //phase 3: compute output offset for each thread in warp and each warp in thread block, then output valid elements + compactK<<>>(d_input,length,d_output,d_BlocksOffset,predicate); + + // determine number of elements in the compacted list + int compact_length = thrustPrt_bOffset[numBlocks-1] + thrustPrt_bCount[numBlocks-1]; + + cudaFree(d_BlocksCount); + cudaFree(d_BlocksOffset); + + return compact_length; +} + + + +} /* namespace cuCompactor */ +#endif /* CUCOMPACTOR_H_ */ diff --git a/alg/teca_detect_nodes.cxx b/alg/teca_detect_nodes.cxx new file mode 100644 index 000000000..994fa4b87 --- /dev/null +++ b/alg/teca_detect_nodes.cxx @@ -0,0 +1,1937 @@ +#include "teca_detect_nodes.h" + +#include "teca_variant_array.h" +#include "teca_variant_array_impl.h" + +#include +#include +#include +#include +#include +#include +#include +#include "STLStringHelper.h" + +#include +#include +#include +#include +#include + +#if defined(TECA_HAS_BOOST) +#include +#endif + +#if defined(TECA_HAS_CUDA) +#include "teca_cuda_util.h" +#include "cuCompactor.cuh" +#endif + +#define TECA_DEBUG 1 + +using std::string; +using std::vector; +using std::cerr; +using std::endl; +using std::cos; + +// PIMPL idiom hides internals +// defines the API for reduction operators +class teca_detect_nodes::internals_t +{ +public: + internals_t() : vecConnectivity(nullptr) {} + ~internals_t() {} + + template + static + int string_parsing(std::set &req_arrays, + VariableRegistry &varreg, + std::string &cmd, + std::vector &vec); + +public: + VariableRegistry varreg; + + std::mutex m_mutex; + SimpleGrid grid; + p_teca_variant_array vecConnectivity; + + std::vector vec_closed_contour_op; + std::vector vec_no_closed_contour_op; + std::vector vec_threshold_op; + std::vector vec_output_op; + + std::string str_search_by; + bool f_search_by_minima; + + std::set dependent_variables; +}; + +// -------------------------------------------------------------------------- +/** + * Determine if the given field has a closed contour about this point. + */ +template +bool has_closed_contour( + const SimpleGrid &grid, + const DataArray1D &dataState, + const int ix0, + double dDeltaAmt, + double dDeltaDist, + double dMinMaxDist) +{ + // Verify arguments + if (dDeltaAmt == 0.0) + { + _EXCEPTIONT("Closed contour amount must be non-zero"); + } + if (dDeltaDist <= 0.0) + { + _EXCEPTIONT("Closed contour distance must be positive"); + } + + // Find min/max near point + int ixOrigin; + + if (dMinMaxDist == 0.0) + { + ixOrigin = ix0; + } + // Find a local minimum / maximum + else + { + real dValue; + float dR; + + FindLocalMinMax( + grid, + (dDeltaAmt > 0.0), + dataState, + ix0, + dMinMaxDist, + ixOrigin, + dValue, + dR); + } + + // Set of visited nodes + std::set setNodesVisited; + + // Set of nodes to visit + std::queue queueToVisit; + queueToVisit.push(ixOrigin); + + // Reference value + real dRefValue = dataState[ixOrigin]; + + const double dLat0 = grid.m_dLat[ixOrigin]; + const double dLon0 = grid.m_dLon[ixOrigin]; + + Announce(2, "Checking (%lu) : (%1.5f %1.5f)", + ixOrigin, dLat0, dLon0); + + // Build up nodes + while (queueToVisit.size() != 0) + { + int ix = queueToVisit.front(); + queueToVisit.pop(); + + if (setNodesVisited.find(ix) != setNodesVisited.end()) + { + continue; + } + + setNodesVisited.insert(ix); + + // Great circle distance to this element + double dLatThis = grid.m_dLat[ix]; + double dLonThis = grid.m_dLon[ix]; + + double dR = + sin(dLat0) * sin(dLatThis) + + cos(dLat0) * cos(dLatThis) * cos(dLonThis - dLon0); + + if (dR >= 1.0) + { + dR = 0.0; + } + else if (dR <= -1.0) + { + dR = 180.0; + } + else + { + dR = 180.0 / M_PI * acos(dR); + } + if (dR != dR) + { + _EXCEPTIONT("NaN value detected"); + } + + Announce(2, "-- (%lu) : (%1.5f %1.5f) : dx %1.5f", + ix, dLatThis, dLonThis, dR); + + // Check great circle distance + if (dR > dDeltaDist) + { + Announce(2, "Failed criteria; returning"); + AnnounceEndBlock(2, NULL); + return false; + } + // Verify sufficient increase in value + if (dDeltaAmt > 0.0) + { + if (dataState[ix] - dRefValue >= dDeltaAmt) + { + continue; + } + } + // Verify sufficient decrease in value + else + { + if (dRefValue - dataState[ix] >= -dDeltaAmt) + { + continue; + } + } + + // Add all neighbors of this point + for (long unsigned int n = 0; n < grid.m_vecConnectivity[ix].size(); ++n) + { + queueToVisit.push(grid.m_vecConnectivity[ix][n]); + } + } + + // Report success with criteria + Announce(2, "Passed criteria; returning"); + AnnounceEndBlock(2, NULL); + return true; +} + +// -------------------------------------------------------------------------- +/** + * Determine if the given field satisfies the threshold. + */ +template +bool satisfies_threshold( + const SimpleGrid &grid, + const DataArray1D &dataState, + const int ix0, + const ThresholdOp::Operation op, + const double dTargetValue, + const double dMaxDist) +{ + // Verify that dMaxDist is less than 180.0 + if (dMaxDist > 180.0) + { + _EXCEPTIONT("MaxDist must be less than 180.0"); + } + + // Queue of nodes that remain to be visited + std::queue queueNodes; + queueNodes.push(ix0); + + // Set of nodes that have already been visited + std::set setNodesVisited; + + // Latitude and longitude at the origin + double dLat0 = grid.m_dLat[ix0]; + double dLon0 = grid.m_dLon[ix0]; + + // Loop through all latlon elements + while (queueNodes.size() != 0) + { + int ix = queueNodes.front(); + queueNodes.pop(); + + if (setNodesVisited.find(ix) != setNodesVisited.end()) + { + continue; + } + + setNodesVisited.insert(ix); + + // Great circle distance to this element + double dLatThis = grid.m_dLat[ix]; + double dLonThis = grid.m_dLon[ix]; + + double dR = + sin(dLat0) * sin(dLatThis) + + cos(dLat0) * cos(dLatThis) * cos(dLonThis - dLon0); + + if (dR >= 1.0) + { + dR = 0.0; + } + else if (dR <= -1.0) + { + dR = 180.0; + } + else + { + dR = 180.0 / M_PI * acos(dR); + } + if (dR != dR) + { + _EXCEPTIONT("NaN value detected"); + } + if ((ix != ix0) && (dR > dMaxDist)) + { + continue; + } + + // Value at this location + double dValue = dataState[ix]; + + // Apply operator + if (op == ThresholdOp::GreaterThan) + { + if (dValue > dTargetValue) + { + return true; + } + } + else if (op == ThresholdOp::LessThan) + { + if (dValue < dTargetValue) + { + return true; + } + } + else if (op == ThresholdOp::GreaterThanEqualTo) + { + if (dValue >= dTargetValue) + { + return true; + } + } + else if (op == ThresholdOp::LessThanEqualTo) + { + if (dValue <= dTargetValue) + { + return true; + } + } + else if (op == ThresholdOp::EqualTo) + { + if (dValue == dTargetValue) + { + return true; + } + } else if (op == ThresholdOp::NotEqualTo) + { + if (dValue != dTargetValue) + { + return true; + } + } + else + { + _EXCEPTIONT("Invalid operation"); + } + + // Special case: zero distance + if (dMaxDist == 0.0) + { + return false; + } + + // Add all neighbors of this point + for (long unsigned int n = 0; n < grid.m_vecConnectivity[ix].size(); ++n) + { + queueNodes.push(grid.m_vecConnectivity[ix][n]); + } + } + + return false; +} + +#if defined(TECA_HAS_CUDA) +struct int_predicate +{ + __host__ __device__ + bool operator()(const int x) + { + return x>=0; + } +}; +namespace cuda_gpu +{ +// -------------------------------------------------------------------------- +template +__global__ +void generate_rectilinear_connectivity( + T *p_vecConnectivity, + unsigned long nLat, + unsigned long nLon, + bool fRegional, + bool fDiagonalConnectivity) +{ + unsigned long q = teca_cuda_util::thread_id_to_array_index(); + + if (q >= nLat*nLon) + return; + + unsigned long i = q % nLon; //lon index + unsigned long j = q / nLon; //lat index + + int neighbors = 0; + + // Connectivity in eight directions + if (fDiagonalConnectivity) + { + if (fRegional) + { + for (int ix = -1; ix <= 1; ix++) + { + for (int jx = -1; jx <= 1; jx++) + { + if ((ix == 0) && (jx == 0)) + { + continue; + } + + int inew = i + ix; + int jnew = j + jx; + + if ((inew < 0) || (inew >= nLon)) + { + neighbors++; + continue; + } + if ((jnew < 0) || (jnew >= nLat)) + { + neighbors++; + continue; + } + p_vecConnectivity[neighbors * nLat * nLon + q] = jnew * nLon + inew; + neighbors++; + } + } + } else { + for (int ix = -1; ix <= 1; ix++) + { + for (int jx = -1; jx <= 1; jx++) + { + if ((ix == 0) && (jx == 0)) + { + continue; + } + + int inew = i + ix; + int jnew = j + jx; + + if ((jnew < 0) || (jnew >= nLat)) + { + neighbors++; + continue; + } + if (inew < 0) + { + inew += nLon; + } + if (inew >= nLon) + { + inew -= nLon; + } + + p_vecConnectivity[neighbors * nLat * nLon + q] = jnew * nLon + inew; + neighbors++; + } + } + } + + // Connectivity in the four primary directions + } else { + if (j != 0) + { + p_vecConnectivity[neighbors * nLat * nLon + q] = (j-1) * nLon + i; + neighbors++; + } + + if (j != nLat-1) + { + p_vecConnectivity[neighbors * nLat * nLon + q] = (j+1) * nLon + i; + neighbors++; + } + + if ((!fRegional) || ((i != 0) && (i != nLon-1))) + { + p_vecConnectivity[neighbors * nLat * nLon + q] = j * nLon + ((i + 1) % nLon); + neighbors++; + p_vecConnectivity[neighbors * nLat * nLon + q] = j * nLon + ((i + nLon - 1) % nLon); + } + //if (fRegional) + //{ + // if (i != 0) + // { + // p_vecConnectivity[neighbors * nLat * nLon + q] = q - 1; + // } + // neighbors++; + + // if (i != nLon-1) + // { + // p_vecConnectivity[neighbors * nLat * nLon + q] = q + 1; + // } + //} else { + // p_vecConnectivity[neighbors * nLat * nLon + q] = j * nLon + ((i + nLon - 1) % nLon); + // neighbors++; + // p_vecConnectivity[neighbors * nLat * nLon + q] = j * nLon + ((i + 1) % nLon); + //} + } +} +// -------------------------------------------------------------------------- +template +__global__ +void find_all_local_minima( + T *p_vecConnectivity, + const T_SEARCH *p_data, + T_SEARCH *p_setMinima, + int sNeighbors, + unsigned long nLat, + unsigned long nLon) +{ + unsigned long q = teca_cuda_util::thread_id_to_array_index(); + + if (q >= nLat*nLon) + return; + + bool fMinimum = true; + + for (int n = 0; n < sNeighbors; n++) + { + if (p_vecConnectivity[n * nLat * nLon + q] >= 0) + { + if (p_data[int(p_vecConnectivity[n * nLat * nLon + q])] < p_data[q]) + { + fMinimum = false; + break; + } + } + } + + if (fMinimum) + { + p_setMinima[q] = q; + } +} +// -------------------------------------------------------------------------- +template +__global__ +void find_all_local_maxima( + T *p_vecConnectivity, + const T_SEARCH *p_data, + T_SEARCH *p_setMaxima, + int sNeighbors, + unsigned long nLat, + unsigned long nLon) +{ + unsigned long q = teca_cuda_util::thread_id_to_array_index(); + + if (q >= nLat*nLon) + return; + + bool fMaximum = true; + + for (int n = 0; n < sNeighbors; n++) + { + if (p_vecConnectivity[n * nLat * nLon + q] >= 0) + { + if (p_data[int(p_vecConnectivity[n * nLat * nLon + q])] > p_data[q]) + { + fMaximum = false; + break; + } + } + } + + if (fMaximum) + { + p_setMaxima[q] = q; + } +} +// -------------------------------------------------------------------------- +template +__global__ +void find_all_local_minmax_with_threshold( + T *p_vecConnectivity, + const T_SEARCH *p_data, + T_SEARCH *p_setMinima, + int sNeighbors, + unsigned long nLat, + unsigned long nLon, + bool fMinima, + ThresholdOp::Operation opThreshold, + double dThresholdValue) +{ + unsigned long q = teca_cuda_util::thread_id_to_array_index(); + + if (q >= nLat*nLon) + return; + + T_SEARCH dSign = (fMinima)?(-1.0):(1.0); + + T_SEARCH dValue = p_data[q]; + + bool fThreshold = true; + if (opThreshold != ThresholdOp::NoThreshold) + { + if (opThreshold == ThresholdOp::GreaterThan) + { + if (dValue <= static_cast(dThresholdValue)) + { + fThreshold = false; + } + + } else if (opThreshold == ThresholdOp::LessThan) + { + if (dValue >= static_cast(dThresholdValue)) + { + fThreshold = false; + } + } else if (opThreshold == ThresholdOp::GreaterThanEqualTo) + { + if (dValue < static_cast(dThresholdValue)) + { + fThreshold = false; + } + } else if (opThreshold == ThresholdOp::LessThanEqualTo) + { + if (dValue > static_cast(dThresholdValue)) + { + fThreshold = false; + } + } else if (opThreshold == ThresholdOp::EqualTo) + { + if (dValue != static_cast(dThresholdValue)) + { + fThreshold = false; + } + } else if (opThreshold == ThresholdOp::NotEqualTo) + { + if (dValue == static_cast(dThresholdValue)) + { + fThreshold = false; + } + } else { + fThreshold = false; + } + } + + if (fThreshold) + { + bool fExtrema = true; + for (int n = 0; n < sNeighbors; n++) + { + if (p_vecConnectivity[n * nLat * nLon + q] >= 0) + { + if (dSign * p_data[int(p_vecConnectivity[n * nLat * nLon + q])] > dSign * dValue) + { + fExtrema = false; + break; + } + } + } + + if (fExtrema) + { + p_setMinima[q] = q; + } + } +} +// -------------------------------------------------------------------------- +template +int generate_rectilinear_connectivity( + int device_id, + T *p_vecConnectivity, + unsigned long nLat, + unsigned long nLon, + bool fRegional, + bool fDiagonalConnectivity) +{ + // determine kernel launch parameters + int n_blocks = 0; + dim3 block_grid; + dim3 thread_grid; + + if (teca_cuda_util::partition_thread_blocks(device_id, + nLat*nLon, 8, block_grid, n_blocks, thread_grid)) + { + TECA_ERROR("Failed to partition thread blocks") + return -1; + } + + cudaError_t ierr = cudaSuccess; + generate_rectilinear_connectivity<<>>(p_vecConnectivity, + nLat, nLon, fRegional, fDiagonalConnectivity); + if ((ierr = cudaGetLastError()) != cudaSuccess) + { + TECA_ERROR("Failed to launch the generate_rectilinear_connectivity CUDA kernel: " + << cudaGetErrorString(ierr)) + return -1; + } + + return 0; +} +// -------------------------------------------------------------------------- +template +int find_all_local_minima( + int device_id, + T *p_vecConnectivity, + const T_SEARCH *p_data, + T_SEARCH *p_setMinima, + T_SEARCH *p_setMinimaCompact, + int sNeighbors, + unsigned long nLat, + unsigned long nLon) +{ + // determine kernel launch parameters + int n_blocks = 0; + dim3 block_grid; + dim3 thread_grid; + + if (teca_cuda_util::partition_thread_blocks(device_id, + nLat*nLon, 8, block_grid, n_blocks, thread_grid)) + { + TECA_ERROR("Failed to partition thread blocks") + return -1; + } + + cudaError_t ierr = cudaSuccess; + find_all_local_minima<<>>(p_vecConnectivity, p_data, + p_setMinima, sNeighbors, nLat, nLon); + if ((ierr = cudaGetLastError()) != cudaSuccess) + { + TECA_ERROR("Failed to launch the find_all_local_minima CUDA kernel: " + << cudaGetErrorString(ierr)) + return -1; + } + + int compact_length = cuCompactor::compact(p_setMinima, p_setMinimaCompact, + nLat*nLon, int_predicate(), thread_grid.x); + + return compact_length; +} +// -------------------------------------------------------------------------- +template +int find_all_local_maxima( + int device_id, + T *p_vecConnectivity, + const T_SEARCH *p_data, + T_SEARCH *p_setMaxima, + T_SEARCH *p_setMaximaCompact, + int sNeighbors, + unsigned long nLat, + unsigned long nLon) +{ + // determine kernel launch parameters + int n_blocks = 0; + dim3 block_grid; + dim3 thread_grid; + + if (teca_cuda_util::partition_thread_blocks(device_id, + nLat*nLon, 8, block_grid, n_blocks, thread_grid)) + { + TECA_ERROR("Failed to partition thread blocks") + return -1; + } + + cudaError_t ierr = cudaSuccess; + find_all_local_maxima<<>>(p_vecConnectivity, p_data, + p_setMaxima, sNeighbors, nLat, nLon); + if ((ierr = cudaGetLastError()) != cudaSuccess) + { + TECA_ERROR("Failed to launch the find_all_local_maxima CUDA kernel: " + << cudaGetErrorString(ierr)) + return -1; + } + + int compact_length = cuCompactor::compact(p_setMaxima, p_setMaximaCompact, + nLat*nLon, int_predicate(), thread_grid.x); + + return compact_length; +} +// -------------------------------------------------------------------------- +template +int find_all_local_minmax_with_threshold( + int device_id, + T *p_vecConnectivity, + const T_SEARCH *p_data, + T_SEARCH *p_setMinima, + T_SEARCH *p_setMinimaCompact, + int sNeighbors, + unsigned long nLat, + unsigned long nLon, + bool fMinima, + std::string search_by_threshold) +{ + if (search_by_threshold.length() != 0) + { + if (search_by_threshold.length() < 2) + { + TECA_ERROR("Threshold operator \"" + << search_by_threshold.c_str() + << "\" must be of form \"\""); + return -1; + } + + ThresholdOp::Operation opThreshold = ThresholdOp::NoThreshold; + double dThresholdValue = 0.0; + + std::string strThreshold1 = search_by_threshold.substr(0,1); + std::string strThreshold2 = search_by_threshold.substr(0,2); + std::string strThresholdValue; + + if (strThreshold2 == ">=") + { + opThreshold = ThresholdOp::GreaterThanEqualTo; + strThresholdValue = search_by_threshold.substr(2); + } else if (strThreshold2 == "<=") + { + opThreshold = ThresholdOp::LessThanEqualTo; + strThresholdValue = search_by_threshold.substr(2); + } else if (strThreshold2 == "!=") + { + opThreshold = ThresholdOp::NotEqualTo; + strThresholdValue = search_by_threshold.substr(2); + } else if (strThreshold1 == ">") + { + opThreshold = ThresholdOp::GreaterThan; + strThresholdValue = search_by_threshold.substr(1); + } else if (strThreshold1 == "<") + { + opThreshold = ThresholdOp::LessThan; + strThresholdValue = search_by_threshold.substr(1); + } else if (strThreshold1 == "=") + { + opThreshold = ThresholdOp::EqualTo; + strThresholdValue = search_by_threshold.substr(1); + } else { + TECA_ERROR("Threshold operator \"" + << search_by_threshold.c_str() + << "\" must be of form \"\"" + << "where is one of >=,<=,!=,>,<,="); + return -1; + } + + // Extract value + if (!STLStringHelper::IsFloat(strThresholdValue)) + { + TECA_ERROR("Threshold operator \"" + << search_by_threshold.c_str() + << "\" must be of form \"\"" + << "where is a floating point number"); + return -1; + } + dThresholdValue = static_cast(std::stod(strThresholdValue)); + + // determine kernel launch parameters + int n_blocks = 0; + dim3 block_grid; + dim3 thread_grid; + + if (teca_cuda_util::partition_thread_blocks(device_id, + nLat*nLon, 8, block_grid, n_blocks, thread_grid)) + { + TECA_ERROR("Failed to partition thread blocks") + return -1; + } + + cudaError_t ierr = cudaSuccess; + find_all_local_minmax_with_threshold<<>>(p_vecConnectivity, + p_data, p_setMinima, sNeighbors, nLat, nLon, + fMinima, opThreshold, dThresholdValue); + if ((ierr = cudaGetLastError()) != cudaSuccess) + { + TECA_ERROR("Failed to launch the " + << "find_all_local_minmax_with_threshold CUDA kernel: " + << cudaGetErrorString(ierr)) + return -1; + } + + int compact_length = cuCompactor::compact(p_setMinima, p_setMinimaCompact, + nLat*nLon, int_predicate(), thread_grid.x); + + return compact_length; + } + + return -1; +} +} //end of namespace cuda +#endif + +// -------------------------------------------------------------------------- +int teca_detect_nodes::detect_cyclones_unstructured( + int device_id, + const_p_teca_cartesian_mesh mesh, + SimpleGrid &grid, + std::set &set_candidates) +{ + // get coordinate arrays + const_p_teca_variant_array y = mesh->get_y_coordinates(); + const_p_teca_variant_array x = mesh->get_x_coordinates(); + + if (!y || !x) + { + TECA_FATAL_ERROR("mesh coordinates are missing.") + return -1; + } + + // get search_by array + const_p_teca_variant_array search_by = + mesh->get_point_arrays()->get(this->internals->str_search_by); + + if (!search_by) + { + TECA_FATAL_ERROR("Dataset missing search_by variable \"" + << this->internals->str_search_by << "\"") + return -1; + } + + VARIANT_ARRAY_DISPATCH_FP(y.get(), + { + std::lock_guard lock(this->internals->m_mutex); + // initialize grid on the first pass. + if (!grid.IsInitialized()) + { + auto [sp_y, p_y] = get_host_accessible(y); + DataArray1D vec_lat(y->size(), false); + vec_lat.AttachToData((void*)p_y); + for (long unsigned int j = 0; j < y->size(); ++j) + { + vec_lat[j] *= M_PI / 180.0; + } + + assert_type(x); + auto [sp_x, p_x] = get_host_accessible(x); + DataArray1D vec_lon(x->size(), false); + vec_lon.AttachToData((void*)p_x); + for (long unsigned int i = 0; i < x->size(); ++i) + { + vec_lon[i] *= M_PI / 180.0; + } + + // No connectivity file; check for latitude/longitude dimension + if (this->in_connect == "") + { + AnnounceStartBlock("Generating RLL grid data"); + grid.GenerateLatitudeLongitude( + vec_lat, vec_lon, this->regional, this->diag_connect, false); + AnnounceEndBlock("Done"); + } + // Check for connectivity file + else + { + TECA_ERROR("Loading grid data from connectivity file" + " is not supported") + return -1; + } + } + } +#if defined(TECA_HAS_CUDA) + if (device_id >= 0) + { + // GPU + // No connectivity file; check for latitude/longitude dimension + if (this->in_connect == "") + { + int neighbors; + if (this->diag_connect) + { + neighbors = 8; + } else { + neighbors = 4; + } + + { + std::lock_guard lock(this->internals->m_mutex); + // initialize grid on the first pass. + if (this->internals->vecConnectivity == nullptr) + { + this->internals->vecConnectivity = teca_variant_array_impl::New( + y->size()*x->size()*neighbors, -1., allocator::cuda_async); + auto [p_vecConnectivity] = data(this->internals->vecConnectivity); + + // Generate connectivity + cuda_gpu::generate_rectilinear_connectivity(device_id, p_vecConnectivity, + y->size(), x->size(), this->regional, this->diag_connect); + } + } + auto [p_vecConnectivity] = data(this->internals->vecConnectivity); + + NESTED_VARIANT_ARRAY_DISPATCH_FP(search_by.get(), _SEARCH, + + auto [sp_data_search, p_data_search] = get_cuda_accessible(search_by); + + p_teca_variant_array candidates = teca_variant_array_impl::New( + y->size()*x->size(), -1., allocator::cuda_async); + auto [p_candidates] = data(candidates); + + int compact_length; + p_teca_variant_array candidates_compact = teca_variant_array_impl::New( + y->size()*x->size(), -1., allocator::cuda_async); + auto [p_candidates_compact] = data(candidates_compact); + + // Tag all minima + //AnnounceStartBlock("FindAllLocalMinMax"); + if (this->search_by_threshold == "") + { + if (this->internals->f_search_by_minima) + compact_length = cuda_gpu::find_all_local_minima( + device_id, p_vecConnectivity, p_data_search, + p_candidates, p_candidates_compact, + neighbors, y->size(), x->size()); + else + compact_length = cuda_gpu::find_all_local_maxima( + device_id, p_vecConnectivity, p_data_search, + p_candidates, p_candidates_compact, + neighbors, y->size(), x->size()); + } else { + compact_length = cuda_gpu::find_all_local_minmax_with_threshold( + device_id, p_vecConnectivity, p_data_search, + p_candidates, p_candidates_compact, + neighbors, y->size(), x->size(), + this->internals->f_search_by_minima, + this->search_by_threshold); + } + + if (compact_length < 0) + { + TECA_ERROR("Failed to find all local min/max") + return -1; + } + + auto [sp_cand, p_cand] = get_host_accessible(candidates_compact); + sync_host_access_any(candidates_compact); + set_candidates.insert(p_cand, p_cand+(compact_length)); + ) + } + } + else + { +#endif + NESTED_VARIANT_ARRAY_DISPATCH_FP(search_by.get(), _SEARCH, + + // CPU + DataArray1D data_search(search_by->size(), false); + auto [sp_search_by, p_search_by] = get_host_accessible(search_by); + data_search.AttachToData((void*)p_search_by); + + // Tag all minima + //AnnounceStartBlock("FindAllLocalMinMax"); + if (this->search_by_threshold == "") + { + if (this->internals->f_search_by_minima) + FindAllLocalMinima(grid, data_search, set_candidates); + else + FindAllLocalMaxima(grid, data_search, set_candidates); + } else { + FindAllLocalMinMaxWithThreshold( + grid, + data_search, + this->internals->f_search_by_minima, + this->search_by_threshold, + set_candidates); + } + //AnnounceEndBlock("Done"); + ) +#if defined(TECA_HAS_CUDA) + } +#endif + ) + + int n_rejected_merge = 0; + VARIANT_ARRAY_DISPATCH_FP(search_by.get(), + + DataArray1D data_search(search_by->size(), false); + auto [sp_search_by, p_search_by] = get_host_accessible(search_by); + data_search.AttachToData((void*)p_search_by); + + // Eliminate based on merge distance + //AnnounceStartBlock("Eliminate based on merge distance"); + if (this->merge_dist != 0.0) + { + std::set set_new_candidates; + + // Calculate chord distance + double d_sph_dist = 2.0 * sin(0.5 * this->merge_dist / 180.0 * M_PI); + + // Create a new KD Tree containing all nodes + kdtree * kdMerge = kd_create(3); + if (kdMerge == NULL) + { + _EXCEPTIONT("kd_create(3) failed"); + } + + std::set::const_iterator iter_candidate = set_candidates.begin(); + for (; iter_candidate != set_candidates.end(); ++iter_candidate) + { + double d_lat = grid.m_dLat[*iter_candidate]; + double d_lon = grid.m_dLon[*iter_candidate]; + + double dx = cos(d_lon) * cos(d_lat); + double dy = sin(d_lon) * cos(d_lat); + double dz = sin(d_lat); + + kd_insert3(kdMerge, dx, dy, dz, (void*)(&(*iter_candidate))); + } + + // Loop through all candidates find set of nearest neighbors + iter_candidate = set_candidates.begin(); + for (; iter_candidate != set_candidates.end(); ++iter_candidate) + { + double d_lat = grid.m_dLat[*iter_candidate]; + double d_lon = grid.m_dLon[*iter_candidate]; + + double dx = cos(d_lon) * cos(d_lat); + double dy = sin(d_lon) * cos(d_lat); + double dz = sin(d_lat); + + // Find all neighbors within d_sph_dist + kdres * kdresMerge = + kd_nearest_range3(kdMerge, dx, dy, dz, d_sph_dist); + + // Number of neighbors + int n_neighbors = kd_res_size(kdresMerge); + if (n_neighbors == 0) + { + set_new_candidates.insert(*iter_candidate); + } + else + { + double d_value = + static_cast(data_search[*iter_candidate]); + + bool f_extrema = true; + for (;;) + { + int * ppr = (int *)(kd_res_item_data(kdresMerge)); + + if (this->internals->f_search_by_minima) + { + if (static_cast(data_search[*ppr]) < d_value) + { + f_extrema = false; + break; + } + + } + else + { + if (static_cast(data_search[*ppr]) > d_value) + { + f_extrema = false; + break; + } + } + + int i_has_more = kd_res_next(kdresMerge); + if (!i_has_more) + { + break; + } + } + + if (f_extrema) + { + set_new_candidates.insert(*iter_candidate); + } + else + { + n_rejected_merge++; + } + } + + kd_res_free(kdresMerge); + } + + // Destroy the KD Tree + kd_free(kdMerge); + + // Update set of pressure minima + set_candidates = set_new_candidates; + } + //AnnounceEndBlock("Done"); + ) + + // Eliminate based on interval + //AnnounceStartBlock("Eliminate based on interval"); + int n_rejected_location = 0; + if ((this->min_lat != this->max_lat) || + (this->min_lon != this->max_lon) || + (this->min_abs_lat != 0.0)) + { + std::set set_new_candidates; + + std::set::const_iterator iter_candidate = set_candidates.begin(); + for (; iter_candidate != set_candidates.end(); ++iter_candidate) + { + double d_lat = grid.m_dLat[*iter_candidate]; + double d_lon = grid.m_dLon[*iter_candidate]; + + if (this->min_lat != this->max_lat) + { + if (d_lat < this->min_lat) + { + n_rejected_location++; + continue; + } + if (d_lat > this->max_lat) + { + n_rejected_location++; + continue; + } + } + if (this->min_lon != this->max_lon) + { + if (d_lon < 0.0) + { + int i_lon_shift = static_cast(d_lon / (2.0 * M_PI)); + d_lon += static_cast(i_lon_shift + 1) * 2.0 * M_PI; + } + if (d_lon >= 2.0 * M_PI) + { + int i_lon_shift = static_cast(d_lon / (2.0 * M_PI)); + d_lon -= static_cast(i_lon_shift - 1) * 2.0 * M_PI; + } + if (this->min_lon < this->max_lon) + { + if (d_lon < this->min_lon) + { + n_rejected_location++; + continue; + } + if (d_lon > this->max_lon) + { + n_rejected_location++; + continue; + } + } + else + { + if ((d_lon > this->max_lon) && + (d_lon < this->min_lon)) + { + n_rejected_location++; + continue; + } + } + } + if (this->min_abs_lat != 0.0) + { + if (fabs(d_lat) < this->min_abs_lat) + { + n_rejected_location++; + continue; + } + } + set_new_candidates.insert(*iter_candidate); + } + set_candidates = set_new_candidates; + } + //AnnounceEndBlock("Done"); + + // Eliminate based on thresholds + //AnnounceStartBlock("Eliminate based on thresholds"); + DataArray1D vec_rejected_threshold( + this->internals->vec_threshold_op.size()); + for (long unsigned int tc = 0; tc < this->internals->vec_threshold_op.size(); ++tc) + { + std::set set_new_candidates; + + // Load the search variable data + Variable &var = + this->internals->varreg.Get( + this->internals->vec_threshold_op[tc].m_varix); + const_p_teca_variant_array threshold_var = + mesh->get_point_arrays()->get(var.GetName()); + + if (!threshold_var) + { + TECA_FATAL_ERROR("Dataset missing variable \"" + << var.GetName() << "\"") + return -1; + } + + VARIANT_ARRAY_DISPATCH_FP(threshold_var.get(), + + DataArray1D data_state(threshold_var->size(), false); + auto [sp_threshold_var, p_threshold_var] = get_host_accessible(threshold_var); + data_state.AttachToData((void*)p_threshold_var); + + // Loop through all pressure minima + std::set::const_iterator iter_candidate = set_candidates.begin(); + for (; iter_candidate != set_candidates.end(); ++iter_candidate) + { + // Determine if the threshold is satisfied + bool f_satisfies_threshold = satisfies_threshold( + grid, data_state, *iter_candidate, + this->internals->vec_threshold_op[tc].m_eOp, + this->internals->vec_threshold_op[tc].m_dValue, + this->internals->vec_threshold_op[tc].m_dDistance); + + // If not rejected, add to new pressure minima array + if (f_satisfies_threshold) + { + set_new_candidates.insert(*iter_candidate); + } + else + { + vec_rejected_threshold[tc]++; + } + } + + set_candidates = set_new_candidates; + ) + } + //AnnounceEndBlock("Done"); + + // Eliminate based on closed contours + //AnnounceStartBlock("Eliminate based on closed contours"); + DataArray1D vec_rejected_closed_contour( + this->internals->vec_closed_contour_op.size()); + for (long unsigned int ccc = 0; ccc < this->internals->vec_closed_contour_op.size(); ++ccc) + { + std::set set_new_candidates; + + // Load the search variable data + Variable &var = + this->internals->varreg.Get( + this->internals->vec_closed_contour_op[ccc].m_varix); + const_p_teca_variant_array closed_contour_var = + mesh->get_point_arrays()->get(var.GetName()); + + if (!closed_contour_var) + { + TECA_FATAL_ERROR("Dataset missing variable \"" + << var.GetName() << "\"") + return -1; + } + + VARIANT_ARRAY_DISPATCH_FP(closed_contour_var.get(), + + DataArray1D data_state(closed_contour_var->size(), false); + auto [sp_closed_contour_var, p_closed_contour_var] = get_host_accessible(closed_contour_var); + data_state.AttachToData((void*)p_closed_contour_var); + + // Loop through all pressure minima + std::set::const_iterator iter_candidate = set_candidates.begin(); + for (; iter_candidate != set_candidates.end(); ++iter_candidate) + { + // Determine if a closed contour is present + bool f_has_closed_contour = has_closed_contour( + grid, data_state, *iter_candidate, + this->internals->vec_closed_contour_op[ccc].m_dDeltaAmount, + this->internals->vec_closed_contour_op[ccc].m_dDistance, + this->internals->vec_closed_contour_op[ccc].m_dMinMaxDist); + + // If not rejected, add to new pressure minima array + if (f_has_closed_contour) + { + set_new_candidates.insert(*iter_candidate); + } + else + { + vec_rejected_closed_contour[ccc]++; + } + } + + set_candidates = set_new_candidates; + ) + } + //AnnounceEndBlock("Done"); + + // Eliminate based on no closed contours + //AnnounceStartBlock("Eliminate based on no closed contours"); + DataArray1D vec_rejected_no_closed_contour( + this->internals->vec_no_closed_contour_op.size()); + for (long unsigned int ccc = 0; ccc < this->internals->vec_no_closed_contour_op.size(); ++ccc) + { + std::set set_new_candidates; + + // Load the search variable data + Variable &var = + this->internals->varreg.Get( + this->internals->vec_no_closed_contour_op[ccc].m_varix); + const_p_teca_variant_array no_closed_contour_var = + mesh->get_point_arrays()->get(var.GetName()); + + if (!no_closed_contour_var) + { + TECA_FATAL_ERROR("Dataset missing variable \"" + << var.GetName() << "\"") + return -1; + } + + VARIANT_ARRAY_DISPATCH_FP(no_closed_contour_var.get(), + + DataArray1D data_state(no_closed_contour_var->size(), false); + auto [sp_no_closed_contour_var, p_no_closed_contour_var] = get_host_accessible(no_closed_contour_var); + data_state.AttachToData((void*)p_no_closed_contour_var); + + // Loop through all pressure minima + std::set::const_iterator iter_candidate = set_candidates.begin(); + for (; iter_candidate != set_candidates.end(); ++iter_candidate) + { + // Determine if a closed contour is present + bool f_has_closed_contour = has_closed_contour( + grid, + data_state, + *iter_candidate, + this->internals->vec_no_closed_contour_op[ccc].m_dDeltaAmount, + this->internals->vec_no_closed_contour_op[ccc].m_dDistance, + this->internals->vec_no_closed_contour_op[ccc].m_dMinMaxDist); + + // If a closed contour is present, reject this candidate + if (f_has_closed_contour) + { + vec_rejected_no_closed_contour[ccc]++; + } + else + { + set_new_candidates.insert(*iter_candidate); + } + } + + set_candidates = set_new_candidates; + ) + } + //AnnounceEndBlock("Done"); + + Announce("Total candidates: %i", set_candidates.size()); + Announce("Rejected ( location): %i", n_rejected_location); + Announce("Rejected ( merged): %i", n_rejected_merge); + + for (long unsigned int tc = 0; tc < vec_rejected_threshold.GetRows(); ++tc) + { + Announce("Rejected (thresh. %s): %i", + this->internals->varreg.GetVariableString( + this->internals->vec_threshold_op[tc].m_varix).c_str(), + vec_rejected_threshold[tc]); + } + + for (long unsigned int ccc = 0; ccc < vec_rejected_closed_contour.GetRows(); ++ccc) + { + Announce("Rejected (contour %s): %i", + this->internals->varreg.GetVariableString( + this->internals->vec_closed_contour_op[ccc].m_varix).c_str(), + vec_rejected_closed_contour[ccc]); + } + + for (long unsigned int ccc = 0; ccc < vec_rejected_no_closed_contour.GetRows(); ++ccc) + { + Announce("Rejected (nocontour %s): %i", + this->internals->varreg.GetVariableString( + this->internals->vec_no_closed_contour_op[ccc].m_varix).c_str(), + vec_rejected_no_closed_contour[ccc]); + } + + return 0; +} + +// -------------------------------------------------------------------------- +teca_detect_nodes::teca_detect_nodes() : + in_connect(""), + search_by_min(""), + search_by_max(""), + closed_contour_cmd(""), + no_closed_contour_cmd(""), + candidate_threshold_cmd(""), + output_cmd(""), + search_by_threshold(""), + min_lon(0.0), + max_lon(0.0), + min_lat(0.0), + max_lat(0.0), + min_abs_lat(0.0), + merge_dist(6.0), + diag_connect(false), + regional(true), + out_header(true) +{ + this->set_number_of_input_connections(1); + this->set_number_of_output_ports(1); + + this->internals = new teca_detect_nodes::internals_t; +} + +// -------------------------------------------------------------------------- +teca_detect_nodes::~teca_detect_nodes() +{ + delete this->internals; +} + +#if defined(TECA_HAS_BOOST) +// -------------------------------------------------------------------------- +void teca_detect_nodes::get_properties_description( + const std::string &prefix, options_description &opts) +{ + options_description ard_opts("Options for " + + (prefix.empty()?"teca_tc_candidates":prefix)); + + ard_opts.add_options() + TECA_POPTS_GET(std::string, prefix, in_connect, + "Connectivity file that describes the unstructured grid") + TECA_POPTS_GET(std::string, prefix, search_by_min, + "Variable to search for the minimum") + TECA_POPTS_GET(std::string, prefix, search_by_max, + "Variable to search for the maximum") + TECA_POPTS_GET(std::string, prefix, closed_contour_cmd, + "Closed contour commands [var,delta,dist,minmaxdist;...]") + TECA_POPTS_GET(std::string, prefix, no_closed_contour_cmd, + "No closed contour commands [var,delta,dist,minmaxdist;...]") + TECA_POPTS_GET(std::string, prefix, candidate_threshold_cmd, + "Threshold commands for candidates [var,op,value,dist;...]") + TECA_POPTS_GET(std::string, prefix, output_cmd, + "Candidates output commands [var,op,dist;...]") + TECA_POPTS_GET(std::string, prefix, search_by_threshold, + "Threshold for search operation in the form ") + TECA_POPTS_GET(double, prefix, min_lon, + "Minimum longitude in degrees for detection") + TECA_POPTS_GET(double, prefix, max_lon, + "Maximum longitude in degrees for detection") + TECA_POPTS_GET(double, prefix, min_lat, + "Minimum latitude in degrees for detection") + TECA_POPTS_GET(double, prefix, max_lat, + "Maximum latitude in degrees for detection") + TECA_POPTS_GET(double, prefix, min_abs_lat, + "Minimum absolute value of latitude in degrees for detection") + TECA_POPTS_GET(double, prefix, merge_dist, + "Minimum allowable distance between two candidates in degrees") + TECA_POPTS_GET(bool, prefix, diag_connect, + "Diagonal connectivity for RLL grids") + TECA_POPTS_GET(bool, prefix, regional, + "Regional (do not wrap longitudinal boundaries)") + TECA_POPTS_GET(bool, prefix, out_header, + "Output header at the beginning of the output file") + ; + + this->teca_algorithm::get_properties_description(prefix, opts); + + opts.add(ard_opts); +} + +// -------------------------------------------------------------------------- +void teca_detect_nodes::set_properties( + const std::string &prefix, variables_map &opts) +{ + this->teca_algorithm::set_properties(prefix, opts); + + TECA_POPTS_SET(opts, std::string, prefix, in_connect) + TECA_POPTS_SET(opts, std::string, prefix, search_by_min) + TECA_POPTS_SET(opts, std::string, prefix, search_by_max) + TECA_POPTS_SET(opts, std::string, prefix, closed_contour_cmd) + TECA_POPTS_SET(opts, std::string, prefix, no_closed_contour_cmd) + TECA_POPTS_SET(opts, std::string, prefix, candidate_threshold_cmd) + TECA_POPTS_SET(opts, std::string, prefix, output_cmd) + TECA_POPTS_SET(opts, std::string, prefix, search_by_threshold) + TECA_POPTS_SET(opts, double, prefix, min_lon) + TECA_POPTS_SET(opts, double, prefix, max_lon) + TECA_POPTS_SET(opts, double, prefix, min_lat) + TECA_POPTS_SET(opts, double, prefix, max_lat) + TECA_POPTS_SET(opts, double, prefix, min_abs_lat) + TECA_POPTS_SET(opts, double, prefix, merge_dist) + TECA_POPTS_SET(opts, bool, prefix, diag_connect) + TECA_POPTS_SET(opts, bool, prefix, regional) + TECA_POPTS_SET(opts, bool, prefix, out_header) +} +#endif + +// -------------------------------------------------------------------------- +template +int teca_detect_nodes::internals_t::string_parsing( + std::set &dep_vars, + VariableRegistry &varreg, + std::string &cmd, + std::vector &vec) +{ + int i_last = 0; + for (long unsigned int i = 0; i <= cmd.length(); ++i) + { + if ((i == cmd.length()) || + (cmd[i] == ';') || + (cmd[i] == ':')) + { + std::string strSubStr = cmd.substr(i_last, i - i_last); + + int i_next_op = (int)(vec.size()); + vec.resize(i_next_op + 1); + vec[i_next_op].Parse(varreg, strSubStr); + + // Load the search variable data + Variable &var = varreg.Get(vec[i_next_op].m_varix); + // Get the data directly from a variable + if (!var.IsOp()) + { + dep_vars.insert(var.GetName()); + } + // Evaluate a data operator to get the contents of this variable + else + { + TECA_ERROR("Data operator is not supported") + return -1; + } + + i_last = i + 1; + } + } + return 0; +} + +// -------------------------------------------------------------------------- +int teca_detect_nodes::initialize() +{ + std::set dep_vars; + + // Only one of search by min or search by max should be specified + if ((this->search_by_min == "") && (this->search_by_max == "")) + { + this->internals->str_search_by = "PSL"; + } + if ((this->search_by_min != "") && (this->search_by_max != "")) + { + TECA_ERROR("Only one of --searchbymin or --searchbymax can" + " be specified"); + return -1; + } + + this->internals->f_search_by_minima = true; + if (this->search_by_min != "") + { + this->internals->str_search_by = this->search_by_min; + this->internals->f_search_by_minima = true; + } + if (this->search_by_max != "") + { + this->internals->str_search_by = this->search_by_max; + this->internals->f_search_by_minima = false; + } + dep_vars.insert(this->internals->str_search_by); + + // Parse the closed contour command string + if (this->closed_contour_cmd != "") + { + if (internals_t::string_parsing(dep_vars, this->internals->varreg, + this->closed_contour_cmd, this->internals->vec_closed_contour_op)) + return -1; + } + + // Parse the no closed contour command string + if (this->no_closed_contour_cmd != "") + { + if (internals_t::string_parsing(dep_vars, this->internals->varreg, + this->no_closed_contour_cmd, this->internals->vec_no_closed_contour_op)) + return -1; + } + + // Parse the threshold operator command string + if (this->candidate_threshold_cmd != "") + { + if (internals_t::string_parsing(dep_vars, this->internals->varreg, + this->candidate_threshold_cmd, this->internals->vec_threshold_op)) + return -1; + } + + // Parse the output operator command string + if (this->output_cmd != "") + { + if (internals_t::string_parsing(dep_vars, this->internals->varreg, + this->output_cmd, this->internals->vec_output_op)) + return -1; + } + + this->internals->dependent_variables = std::move(dep_vars); + + this->min_lon *= M_PI / 180.0; + this->max_lon *= M_PI / 180.0; + this->min_lat *= M_PI / 180.0; + this->max_lat *= M_PI / 180.0; + this->min_abs_lat *= M_PI / 180.0; + + return 0; +} + +// -------------------------------------------------------------------------- +std::vector teca_detect_nodes::get_upstream_request( + unsigned int port, + const std::vector &md_in, + const teca_metadata &req_in) +{ + if (this->get_verbose() > 1) + { + std::cerr << teca_parallel_id() + << "teca_detect_nodes::get_upstream_request" << std::endl; + } + + (void)port; + + std::vector up_reqs; + teca_metadata md = md_in[0]; + + // locate the extents of the user supplied region of + // interest + teca_metadata coords; + if (md.get("coordinates", coords)) + { + TECA_FATAL_ERROR("metadata is missing \"coordinates\"") + return up_reqs; + } + + p_teca_variant_array in_x, in_y, in_z; + if (!(in_x = coords.get("x")) || + !(in_y = coords.get("y")) || + !(in_z = coords.get("z"))) + { + TECA_FATAL_ERROR("metadata missing coordinate arrays") + return up_reqs; + } + + unsigned long extent[6] = {0}; + double req_bounds[6] = {0.0}; + req_bounds[0] = this->min_lon * 180.0 / M_PI; + req_bounds[1] = this->max_lon * 180.0 / M_PI; + req_bounds[2] = this->min_lat * 180.0 / M_PI; + req_bounds[3] = this->max_lat * 180.0 / M_PI; + + if (teca_coordinate_util::bounds_to_extent(req_bounds, + in_x, in_y, in_z, extent) || + teca_coordinate_util::validate_extent(in_x->size(), + in_y->size(), in_z->size(), extent, true)) + { + TECA_FATAL_ERROR("failed to determine the active extent") + return up_reqs; + } + + if (this->get_verbose() > 1) + { + cerr << teca_parallel_id() << "active_bound = " + << this->min_lon<< ", " << this->max_lon + << ", " << this->min_lat << ", " << this->max_lat + << endl; + cerr << teca_parallel_id() << "active_extent = " + << extent[0] << ", " << extent[1] << ", " << extent[2] << ", " + << extent[3] << ", " << extent[4] << ", " << extent[5] << endl; + } + + // get the requested arrays + std::set req_arrays; + req_in.get("arrays", req_arrays); + + req_arrays.insert(this->internals->dependent_variables.begin(), + this->internals->dependent_variables.end()); + + teca_metadata req(req_in); + req.set("arrays", req_arrays); + req.set("extent", extent); + up_reqs.push_back(req); + + return up_reqs; +} + +// -------------------------------------------------------------------------- +const_p_teca_dataset teca_detect_nodes::execute( + unsigned int port, + const std::vector &input_data, + const teca_metadata &req_in) +{ + (void)port; + (void)req_in; + + if (!input_data.size()) + { + TECA_FATAL_ERROR("empty input") + return nullptr; + } + + // get the assigned GPU or CPU + int device_id = -1; + req_in.get("device_id", device_id); + +#if defined(TECA_HAS_CUDA) + if (device_id >= 0) + { + if (teca_cuda_util::set_device(device_id)) + return nullptr; + } +#endif + + if (this->get_verbose() > 1) + { + std::cerr << teca_parallel_id() + << "teca_detect_nodes::execute " + << " device " << device_id << std::endl; + } + + // get the input dataset + const_p_teca_cartesian_mesh mesh + = std::dynamic_pointer_cast(input_data[0]); + if (!mesh) + { + TECA_FATAL_ERROR("teca_cartesian_mesh is required") + return nullptr; + } + + // get time step + unsigned long time_step; + mesh->get_time_step(time_step); + + // get temporal offset of the current timestep + double time_offset = 0.0; + mesh->get_time(time_offset); + + // get offset units + std::string time_units; + mesh->get_time_units(time_units); + + // get offset calendar + std::string calendar; + mesh->get_calendar(calendar); + + std::set set_candidates; + + if (this->detect_cyclones_unstructured(device_id, mesh, this->internals->grid, set_candidates)) + { + TECA_FATAL_ERROR("TC detector encountered an error") + return nullptr; + } + + // Write candidate information + int i_candidate_ix = 0; + + // Apply output operators + std::vector> vec_output_value; + vec_output_value.resize(set_candidates.size()); + for (long unsigned int i = 0; i < set_candidates.size(); ++i) + { + vec_output_value[i].resize(this->internals->vec_output_op.size()); + } + + for (long unsigned int outc = 0; outc < this->internals->vec_output_op.size(); ++outc) + { + Variable &var = + this->internals->varreg.Get( + this->internals->vec_output_op[outc].m_varix); + const_p_teca_variant_array output_var = + mesh->get_point_arrays()->get(var.GetName()); + + if (!output_var) + { + TECA_FATAL_ERROR("Dataset missing variable \"" << + var.GetName() << "\"") + return nullptr; + } + + VARIANT_ARRAY_DISPATCH_FP(output_var.get(), + + DataArray1D data_state(output_var->size(), false); + auto [sp_output_var, p_output_var] = get_host_accessible(output_var); + data_state.AttachToData((void*)p_output_var); + + // Loop through all pressure minima + i_candidate_ix = 0; + std::set::const_iterator iter_candidate = set_candidates.begin(); + for (; iter_candidate != set_candidates.end(); ++iter_candidate) + { + ApplyNodeOutputOp( + this->internals->vec_output_op[outc], + this->internals->grid, + data_state, + *iter_candidate, + vec_output_value[i_candidate_ix][outc]); + + i_candidate_ix++; + } + ) + } + + // Output all candidates + i_candidate_ix = 0; + + double * lat_array = new double[set_candidates.size()]; + double * lon_array = new double[set_candidates.size()]; + int * i_array = new int[set_candidates.size()]; + int * j_array = new int[set_candidates.size()]; + + std::set::const_iterator iter_candidate = set_candidates.begin(); + for (; iter_candidate != set_candidates.end(); ++iter_candidate) + { + if (this->internals->grid.m_nGridDim.size() == 1) + { + i_array[i_candidate_ix] = *iter_candidate; + } + else if (this->internals->grid.m_nGridDim.size() == 2) + { + i_array[i_candidate_ix] = + (*iter_candidate) % static_cast(this->internals->grid.m_nGridDim[1]); + j_array[i_candidate_ix] = + (*iter_candidate) / static_cast(this->internals->grid.m_nGridDim[1]); + } + lat_array[i_candidate_ix] = this->internals->grid.m_dLat[*iter_candidate] * 180.0 / M_PI; + lon_array[i_candidate_ix] = this->internals->grid.m_dLon[*iter_candidate] * 180.0 / M_PI; + + i_candidate_ix++; + } + + // build the output + p_teca_table out_table = teca_table::New(); + out_table->set_calendar(calendar); + out_table->set_time_units(time_units); + + // add time stamp + out_table->declare_columns("step", long(), "time", double(), "ncandidates", long()); + for (unsigned long i = 0; i < set_candidates.size(); ++i) + out_table << time_step << time_offset << set_candidates.size(); + + // put the arrays into the table + p_teca_variant_array_impl latitude = + teca_variant_array_impl::New(set_candidates.size(), lat_array); + p_teca_variant_array_impl longitude = + teca_variant_array_impl::New(set_candidates.size(), lon_array); + p_teca_variant_array_impl i = + teca_variant_array_impl::New(set_candidates.size(), i_array); + + out_table->append_column("i", i); + if (this->internals->grid.m_nGridDim.size() == 2) + { + p_teca_variant_array_impl j = + teca_variant_array_impl::New(set_candidates.size(), j_array); + out_table->append_column("j", j); + } + out_table->append_column("lat", latitude); + out_table->append_column("lon", longitude); + + for (long unsigned int outc = 0; outc < this->internals->vec_output_op.size(); ++outc) + { + i_candidate_ix = 0; + double * output_array = new double[set_candidates.size()]; + std::set::const_iterator iter_candidate = set_candidates.begin(); + for (; iter_candidate != set_candidates.end(); ++iter_candidate) + { + output_array[i_candidate_ix] = std::stod(vec_output_value[i_candidate_ix][outc]); + i_candidate_ix++; + } + + p_teca_variant_array_impl output = + teca_variant_array_impl::New(set_candidates.size(), + output_array); + + Variable &var = this->internals->varreg.Get( + this->internals->vec_output_op[outc].m_varix); + out_table->append_column(var.ToString(this->internals->varreg).c_str(), + output); + delete [] output_array; + } + +#if TECA_DEBUG > 1 + out_table->to_stream(cerr); + cerr << std::endl; +#endif + + delete [] lat_array; + delete [] lon_array; + delete [] i_array; + delete [] j_array; + + return out_table; +} diff --git a/alg/teca_detect_nodes.h b/alg/teca_detect_nodes.h new file mode 100644 index 000000000..9de576aca --- /dev/null +++ b/alg/teca_detect_nodes.h @@ -0,0 +1,213 @@ +#ifndef teca_detect_nodes_h +#define teca_detect_nodes_h + +#include "teca_shared_object.h" +#include "teca_algorithm.h" +#include "teca_metadata.h" +#include "teca_table.h" +#include "teca_coordinate_util.h" + +#include + +#include +#include + +TECA_SHARED_OBJECT_FORWARD_DECL(teca_detect_nodes) + +/** + * Porting the DetectNodes function from TempestExtremes to TECA. + * detect_nodes is used for the detection of nodal features. + * Candidate points are selected based on information at a single time slice. + * Typically detect_nodes is followed by stitch_nodes + * to connect candidate points in time. + */ +class TECA_EXPORT teca_detect_nodes : public teca_algorithm +{ +public: + TECA_ALGORITHM_STATIC_NEW(teca_detect_nodes) + TECA_ALGORITHM_DELETE_COPY_ASSIGN(teca_detect_nodes) + TECA_ALGORITHM_CLASS_NAME(teca_detect_nodes) + ~teca_detect_nodes(); + + // report/initialize to/from Boost program options + // objects. + TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION() + TECA_SET_ALGORITHM_PROPERTIES() + + int initialize(); + + /** @name in_connet + * Set the connectivity file + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, in_connect) + ///@} + + /** @name search_by_min + * Set variable to search for the minimum + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, search_by_min) + ///@} + + /** @name search_by_max + * Set variable to search for the maximum + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, search_by_max) + ///@} + + /** @name closed_contour_cmd + * Set the closed contour commands + * [var,delta,dist,minmaxdist;...] + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, closed_contour_cmd) + ///@} + + /** @name no_closed_contour_cmd + * Set the no closed contour commands + * [var,delta,dist,minmaxdist;...] + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, no_closed_contour_cmd) + ///@} + + /** @name candidate_threshold_cmd + * Set the threshold commands + * [var,op,value,dist;...] + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, candidate_threshold_cmd) + ///@} + + /** @name output_cmd + * Set the output commands + * [var,op,dist;...] + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, output_cmd) + ///@} + + /** @name search_by_threshold + * Set threshold for search operation + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, search_by_threshold) + ///@} + + /** @name min_lon + * Set minimum longitude in degrees for detection + * default 0.0 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, min_lon) + ///@} + + /** @name max_lon + * Set maximum longitude in degrees for detection + * default 0.0 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, max_lon) + ///@} + + /** @name min_lat + * Set minimum latitude in degrees for detection + * default 0.0 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, min_lat) + ///@} + + /** @name max_lat + * Set maximum latitude in degrees for detection + * default 0.0 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, max_lat) + ///@} + + /** @name min_abs_lat + * Set minimum absolute value of latitude in degrees for detection + * default 0.0 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, min_abs_lat) + ///@} + + /** @name merge_dist + * Set minimum allowable distance between two candidates in degrees + * default 6.0 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, merge_dist) + ///@} + + /** @name diag_connect + * Set diagonal connectivity for RLL grids + * default false + */ + ///@{ + TECA_ALGORITHM_PROPERTY(bool, diag_connect) + ///@} + + /** @name regional + * Set regional (do not wrap longitudinal boundaries) + * default true + */ + ///@{ + TECA_ALGORITHM_PROPERTY(bool, regional) + ///@} + + /** @name out_header + * Set output header + * default true + */ + ///@{ + TECA_ALGORITHM_PROPERTY(bool, out_header) + ///@} + +protected: + teca_detect_nodes(); + +private: + int detect_cyclones_unstructured( + int device_id, + const_p_teca_cartesian_mesh mesh, + SimpleGrid &grid, + std::set &setCandidates); + + std::vector get_upstream_request( + unsigned int port, + const std::vector &md_in, + const teca_metadata &req_in) override; + + const_p_teca_dataset execute( + unsigned int port, + const std::vector &input_data, + const teca_metadata &req_in) override; + +private: + std::string in_connect; + std::string search_by_min; + std::string search_by_max; + std::string closed_contour_cmd; + std::string no_closed_contour_cmd; + std::string candidate_threshold_cmd; + std::string output_cmd; + std::string search_by_threshold; + double min_lon; + double max_lon; + double min_lat; + double max_lat; + double min_abs_lat; + double merge_dist; + bool diag_connect; + bool regional; + bool out_header; + + class internals_t; + internals_t *internals; +}; +#endif diff --git a/alg/teca_mesh_join.cxx b/alg/teca_mesh_join.cxx new file mode 100644 index 000000000..1f4961727 --- /dev/null +++ b/alg/teca_mesh_join.cxx @@ -0,0 +1,360 @@ +#include "teca_mesh_join.h" + +#include "teca_cartesian_mesh.h" +#include "teca_array_collection.h" +#include "teca_variant_array.h" +#include "teca_variant_array_impl.h" +#include "teca_variant_array_util.h" +#include "teca_metadata.h" +#include "teca_array_attributes.h" +#include "teca_coordinate_util.h" + +#include +#include +#include +#include +#include + +#if defined(TECA_HAS_BOOST) +#include +#endif + +using namespace teca_variant_array_util; +using allocator = teca_variant_array::allocator; + +//#define TECA_DEBUG + +// -------------------------------------------------------------------------- +teca_mesh_join::teca_mesh_join() +{ + this->set_number_of_input_connections(2); + this->set_number_of_output_ports(1); +} + +// -------------------------------------------------------------------------- +teca_mesh_join::~teca_mesh_join() +{} + +#if defined(TECA_HAS_BOOST) +// -------------------------------------------------------------------------- +void teca_mesh_join::get_properties_description( + const std::string &prefix, options_description &global_opts) +{ + options_description opts("Options for " + + (prefix.empty()?"teca_mesh_join":prefix)); + + this->teca_algorithm::get_properties_description(prefix, opts); + + global_opts.add(opts); +} + +// -------------------------------------------------------------------------- +void teca_mesh_join::set_properties( + const std::string &prefix, variables_map &opts) +{ + this->teca_algorithm::set_properties(prefix, opts); +} +#endif + + +// -------------------------------------------------------------------------- +teca_metadata teca_mesh_join::get_output_metadata( + unsigned int port, const std::vector &input_md) +{ +#ifdef TECA_DEBUG + std::cerr << teca_parallel_id() + << "teca_mesh_join::get_output_metadata" << std::endl; +#endif + (void)port; + + // start with a copy of metadata from the target + unsigned int md_target = 0; + teca_metadata output_md(input_md[md_target]); + + // get target metadata + std::set target_vars; + input_md[md_target].get("variables", target_vars); + + teca_metadata target_atts; + input_md[md_target].get("attributes", target_atts); + + teca_metadata target_coords; + input_md[md_target].get("coordinates", target_coords); + + // work with each source + unsigned int n_in = this->get_number_of_input_connections(); + for (unsigned int i = 1; i < n_in; ++i) + { + unsigned int md_src = i; + + // get source metadata + std::vector source_vars; + input_md[md_src].get("variables", source_vars); + + teca_metadata source_atts; + input_md[md_src].get("attributes", source_atts); + + teca_metadata source_coords; + input_md[md_src].get("coordinates", source_coords); + + // merge metadata from source and target variables should be unique + // lists. attributes are indexed by variable names in the case of + // collisions, the target variable is kept, the source variable is + // ignored + size_t n_source_vars = source_vars.size(); + for (size_t i = 0; i < n_source_vars; ++i) + { + const std::string &src_var = source_vars[i]; + + auto [it, ins] = target_vars.insert(src_var); + + if (ins) + { + teca_metadata atts; + source_atts.get(src_var, atts); + target_atts.set(src_var, atts); + + teca_metadata coords; + source_coords.get(src_var, coords); + target_coords.set(src_var, coords); + } + } + } + + // update with merged lists + output_md.set("variables", target_vars); + output_md.set("attributes", target_atts); + output_md.set("coordinates", target_coords); + + return output_md; +} + +// -------------------------------------------------------------------------- +std::vector teca_mesh_join::get_upstream_request( + unsigned int port, const std::vector &input_md, + const teca_metadata &request) +{ + (void)port; + // route requests for arrays to either target or the input. + // if the array exists in both then it is take from the target + + // start by duplicating the request for each input + unsigned int n_in = this->get_number_of_input_connections(); + std::vector up_reqs(n_in, request); + + // get input metadata + std::vector> input_vars(n_in); + for (unsigned int i = 0; i < n_in; ++i) + { + input_md[i].get("variables", input_vars[i]); + } + + // get the requested arrays + std::vector req_arrays; + request.get("arrays", req_arrays); + + // route the request for each array to the most appropriate input. in the + // case of inputs providing the same array the request is sent to the lower + // input. + std::vector> up_req_arrays(n_in); + + auto it = req_arrays.begin(); + auto end = req_arrays.end(); + for (; it != end; ++it) + { + // work with each input + bool array_found = false; + for (unsigned int i = 0; i < n_in; ++i) + { + // check if the i'th input has the array + if (input_vars[i].count(*it)) + { + // request from the i'th input + up_req_arrays[i].insert(*it); + array_found = true; + break; + } + } + + // require that at least one input can provide the requested array + if (!array_found) + { + TECA_FATAL_ERROR("\"" << *it << "\" was not found on any input") + return {}; + } + } + + // update the requests + for (unsigned int i = 0; i < n_in; ++i) + up_reqs[i].set("arrays", up_req_arrays[i]); + +#ifdef TECA_DEBUG + for (unsigned int i = 0; i < n_in; ++i) + { + std::cerr << "request[" << i << "] = "; + up_reqs[i].to_stream(std::cerr); + std::cerr << std::endl; + } +#endif + + return up_reqs; +} + +// -------------------------------------------------------------------------- +const_p_teca_dataset teca_mesh_join::execute( + unsigned int port, const std::vector &input_data, + const teca_metadata &request) +{ +#ifdef TECA_DEBUG + std::cerr << teca_parallel_id() << "teca_mesh_join::execute" << std::endl; +#endif + (void)port; + + unsigned int n_in = this->get_number_of_input_connections(); + + p_teca_cartesian_mesh in_target + = std::dynamic_pointer_cast( + std::const_pointer_cast(input_data[0])); + + if (!in_target) + { + TECA_FATAL_ERROR("invalid input. target invalid") + return nullptr; + } + + // create the output + p_teca_cartesian_mesh target = teca_cartesian_mesh::New(); + target->shallow_copy(in_target); + + // get the coordinate and data arrays + const_p_teca_variant_array target_xc = target->get_x_coordinates(); + const_p_teca_variant_array target_yc = target->get_y_coordinates(); + const_p_teca_variant_array target_zc = target->get_z_coordinates(); + p_teca_array_collection target_ac = target->get_point_arrays(); + + // get the attributes + teca_metadata target_atts; + target->get_attributes(target_atts); + + unsigned long target_nx = target_xc->size(); + unsigned long target_ny = target_yc->size(); + unsigned long target_nz = target_zc->size(); + + unsigned long target_ihi = target_nx - 1; + unsigned long target_jhi = target_ny - 1; + + double tx0, tx1; + target_xc->get(0, tx0); + target_xc->get(target_ihi, tx1); + + double ty0, ty1; + target_yc->get(0, ty0); + target_yc->get(target_jhi, ty1); + + // get the list of arrays to move + std::vector req_arrays; + request.get("arrays", req_arrays); + + auto it = req_arrays.begin(); + auto end = req_arrays.end(); + for (; it != end; ++it) + { + if (!target->get_point_arrays()->has(*it)) + { + bool array_found = false; + for (unsigned int i = 1; i < n_in; ++i) + { + const_p_teca_cartesian_mesh source + = std::dynamic_pointer_cast + (input_data[i]); + + if (!source) + { + TECA_FATAL_ERROR("invalid input. source invalid") + return nullptr; + } + + if (source->get_point_arrays()->has(*it)) + { + const_p_teca_variant_array source_xc = source->get_x_coordinates(); + const_p_teca_variant_array source_yc = source->get_y_coordinates(); + const_p_teca_variant_array source_zc = source->get_z_coordinates(); + const_p_teca_array_collection source_ac = source->get_point_arrays(); + + unsigned long source_shape[4] = {0}; + source->get_array_shape(*it, source_shape); + + unsigned long source_nx = source_shape[0]; + unsigned long source_ny = source_shape[1]; + unsigned long source_nz = source_shape[2]; + unsigned long source_ihi = source_nx - 1; + unsigned long source_jhi = source_ny - 1; + + double sx0, sx1; + source_xc->get(0, sx0); + source_xc->get(source_ihi, sx1); + + double sy0, sy1; + source_yc->get(0, sy0); + source_yc->get(source_jhi, sy1); + + // copy when the input and output meshes are the same. the meshes are the + // same when they span the same world space and have the same resolution + if ((source_nx == target_nx) && (source_ny == target_ny) && (source_nz == target_nz) && + teca_coordinate_util::equal(sx0, tx0) && teca_coordinate_util::equal(sx1, tx1) && + teca_coordinate_util::equal(sy0, ty0) && teca_coordinate_util::equal(sy1, ty1)) + { + // move the array attributes + teca_metadata source_atts; + source->get_attributes(source_atts); + + teca_metadata array_atts; + if (!source_atts.get(*it, array_atts)) + target_atts.set(*it, array_atts); + + if (this->verbose) + { + TECA_STATUS("\""<< *it << "\". Identical dimensions[x,y,z]: " + << source_nx << "/" << target_nx << ", " + << source_ny << "/" << target_ny << ", " + << source_nz << "/" << target_nz + << " and identical bounds[x,y]: " + << sx0 << "/" << tx0 << ", " << sx1 << "/" << tx1 << ", " + << sy0 << "/" << ty0 << ", " << sy1 << "/" << ty1 + << ". Copying data.") + } + + const_p_teca_variant_array source_a = source_ac->get(*it); + p_teca_variant_array target_a = source_a->new_copy(); + target_ac->set(*it, target_a); + + array_found = true; + break; + } + else + { + TECA_FATAL_ERROR("\""<< *it << "\". Not identical dimensions[x,y,z]: " + << source_nx << "/" << target_nx << ", " + << source_ny << "/" << target_ny << ", " + << source_nz << "/" << target_nz + << " or not identical bounds[x,y]: " + << sx0 << "/" << tx0 << ", " << sx1 << "/" << tx1 << ", " + << sy0 << "/" << ty0 << ", " << sy1 << "/" << ty1) + return nullptr; + } + } + } + + if (!array_found) + { + TECA_FATAL_ERROR("Array \"" << *it + << "\" is not present on any of the inputs") + return nullptr; + } + } + } + + target->set_attributes(target_atts); + + return target; +} diff --git a/alg/teca_mesh_join.h b/alg/teca_mesh_join.h new file mode 100644 index 000000000..2d0afd406 --- /dev/null +++ b/alg/teca_mesh_join.h @@ -0,0 +1,60 @@ +#ifndef teca_mesh_join_h +#define teca_mesh_join_h + +#include "teca_shared_object.h" +#include "teca_algorithm.h" +#include "teca_metadata.h" + +#include +#include + +TECA_SHARED_OBJECT_FORWARD_DECL(teca_mesh_join) + + +/** + * An algorithm that joins two or more mesh data. + * It requires that the mesh geometry be identical in all files. + */ +class TECA_EXPORT teca_mesh_join : public teca_algorithm +{ +public: + TECA_ALGORITHM_STATIC_NEW(teca_mesh_join) + TECA_ALGORITHM_DELETE_COPY_ASSIGN(teca_mesh_join) + TECA_ALGORITHM_CLASS_NAME(teca_mesh_join) + ~teca_mesh_join(); + + // report/initialize to/from Boost program options + // objects. + TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION() + TECA_SET_ALGORITHM_PROPERTIES() + + /** Set the number input connections. The default is 2, this must be called + * when more than 2 meshes are to be merged. + */ + void set_number_of_input_connections(unsigned int n) + { this->teca_algorithm::set_number_of_input_connections(n); } + +protected: + teca_mesh_join(); + +private: + using teca_algorithm::get_output_metadata; + + teca_metadata get_output_metadata( + unsigned int port, + const std::vector &input_md) override; + + std::vector get_upstream_request( + unsigned int port, + const std::vector &input_md, + const teca_metadata &request) override; + + const_p_teca_dataset execute( + unsigned int port, + const std::vector &input_data, + const teca_metadata &request) override; + +private: +}; + +#endif diff --git a/alg/teca_rename_variables.cxx b/alg/teca_rename_variables.cxx index 5bb8849f3..8a45b8690 100644 --- a/alg/teca_rename_variables.cxx +++ b/alg/teca_rename_variables.cxx @@ -137,6 +137,7 @@ teca_metadata teca_rename_variables::get_output_metadata( attributes.set(this->new_variable_names[i], atts); } + out_md.set("variables", out_vars); out_md.set("attributes", attributes); return out_md; @@ -205,7 +206,6 @@ const_p_teca_dataset teca_rename_variables::execute( p_teca_mesh out_mesh = std::static_pointer_cast (std::const_pointer_cast(in_mesh)->new_shallow_copy()); - // rename the arrays if they are found p_teca_array_collection arrays = out_mesh->get_point_arrays(); @@ -220,6 +220,12 @@ const_p_teca_dataset teca_rename_variables::execute( arrays->remove(var_name); arrays->set(this->new_variable_names[i], array); } + else + { + TECA_STATUS("\"" << var_name + << "\" not found to rename to " + << this->new_variable_names[i]) + } } return out_mesh; diff --git a/alg/teca_stitch_nodes.cxx b/alg/teca_stitch_nodes.cxx new file mode 100644 index 000000000..b9bd12b06 --- /dev/null +++ b/alg/teca_stitch_nodes.cxx @@ -0,0 +1,1625 @@ +#include "teca_stitch_nodes.h" + +#include "teca_variant_array_impl.h" +#include "teca_table.h" + +#include +#include "Announce.h" +#include "TimeObj.h" +#include "CoordTransforms.h" +#include "STLStringHelper.h" + +#include +#include + +#if defined(TECA_HAS_BOOST) +#include +#endif + +#if defined(TECA_HAS_CUDA) +#include "teca_cuda_util.h" +#endif + +#define TECA_DEBUG 2 + +using std::cerr; +using std::endl; + +// -------------------------------------------------------------------------- +typedef std::vector> TimesliceCandidateInformation; +typedef std::pair TimeToCandidateInfoPair; +typedef std::map TimeToCandidateInfoMap; +typedef TimeToCandidateInfoMap::iterator TimeToCandidateInfoMapIterator; + +// -------------------------------------------------------------------------- +/** + * Read a DetectNodes output file into mapCandidateInfo + */ +int parse_detect_nodes_file( + const_p_teca_table &candidates, + const std::vector &vec_format_strings, + TimeToCandidateInfoMap &mapCandidateInfo, + Time::CalendarType caltype, + bool allow_repeated_times) +{ + // Current read state + enum ReadState + { + ReadState_Time, + ReadState_Candidate + } eReadState = ReadState_Time; + + // Current candidate at this time + int iCandidate = 0; + // Total candidates at this time + int nCandidates = 0; + + // Current candidate information + TimeToCandidateInfoMapIterator iterCurrentTime = mapCandidateInfo.end(); + + auto year = candidates->get_column("year"); + auto month = candidates->get_column("month"); + auto day = candidates->get_column("day"); + auto hour = candidates->get_column("hour"); + auto minute = candidates->get_column("minute"); + auto n_candidates = candidates->get_column("ncandidates"); + auto lon = candidates->get_column("lon"); + auto lat = candidates->get_column("lat"); + auto ii = candidates->get_column("i"); + auto jj = candidates->get_column("j"); + + unsigned long n_rows = candidates->get_number_of_rows(); + + VARIANT_ARRAY_DISPATCH(lat.get(), + + auto [sp_year, p_year] = get_host_accessible(year); + auto [sp_month, p_month] = get_host_accessible(month); + auto [sp_day, p_day] = get_host_accessible(day); + auto [sp_hour, p_hour] = get_host_accessible(hour); + auto [sp_minute, p_minute] = get_host_accessible(minute); + auto [sp_n_candidates, p_n_candidates] = get_host_accessible(n_candidates); + auto [sp_lat, p_lat] = get_host_accessible(lat); + auto [sp_lon, p_lon] = get_host_accessible(lon); + auto [sp_ii, p_ii] = get_host_accessible(ii); + auto [sp_jj, p_jj] = get_host_accessible(jj); + + for (unsigned long row = 0; row < n_rows; ++row) + { + // Parse the time + if (eReadState == ReadState_Time) + { + iCandidate = 0; + nCandidates = p_n_candidates[row]; + + Time time(caltype); + time.SetYear(p_year[row]); + time.SetMonth(p_month[row]); + time.SetDay(p_day[row]); + time.SetSecond(p_hour[row] * 3600 + p_minute[row] * 60); + + auto it = mapCandidateInfo.find(time); + //if (it == mapCandidateInfo.end()) + if (it != mapCandidateInfo.end()) + { + if (allow_repeated_times) + { + iterCurrentTime = it; + iCandidate = iterCurrentTime->second.size(); + nCandidates += iCandidate; + iterCurrentTime->second.resize(iCandidate + nCandidates); + + } + else + { + TECA_ERROR("Repated time found in candidate files on line " << row) + } + } + else + { + auto ins = mapCandidateInfo.insert(TimeToCandidateInfoPair(time, + TimesliceCandidateInformation(nCandidates))); + + if (!ins.second) + { + TECA_ERROR("Insertion of time " + << " into candidate info map failed on line " << row) + return -1; + } + + iterCurrentTime = ins.first; + } + + // Prepare to parse candidate data + if (nCandidates != 0) + { + eReadState = ReadState_Candidate; + } + } + + if (eReadState == ReadState_Candidate) + { + TimesliceCandidateInformation &tscinfo = iterCurrentTime->second; + + for (size_t i = 0; i < vec_format_strings.size(); ++i) + { + if (vec_format_strings[i] == "lat") + { + tscinfo[iCandidate].push_back(std::to_string(p_lat[row])); + } + else if (vec_format_strings[i] == "lon") + { + tscinfo[iCandidate].push_back(std::to_string(p_lon[row])); + } + else if (vec_format_strings[i] == "i") + { + tscinfo[iCandidate].push_back(std::to_string(p_ii[row])); + } + else if (vec_format_strings[i] == "j") + { + tscinfo[iCandidate].push_back(std::to_string(p_jj[row])); + } + else + { + auto var = candidates->get_column(vec_format_strings[i]); + auto [sp_var, p_var] = get_host_accessible(var); + tscinfo[iCandidate].push_back(std::to_string(p_var[row])); + } + } + + // Advance candidate number + iCandidate++; + if (iCandidate == nCandidates) + { + eReadState = ReadState_Time; + iterCurrentTime = mapCandidateInfo.end(); + } + } + } + ) + + return 0; +} + +// -------------------------------------------------------------------------- +struct Node +{ + double x; + double y; + double z; + + double lonRad; + double latRad; +}; + +// -------------------------------------------------------------------------- +/** + * A (time,candidate index) pair + */ +struct TimeCandidatePair : public std::pair +{ + TimeCandidatePair(int _timeix, int _candidate) : + std::pair(_timeix, _candidate) {} + + inline const int & timeix() const + { + return first; + } + + inline const int & candidate() const + { + return second; + } +}; + +// -------------------------------------------------------------------------- +/** + * A segment in a simple path connecting candidate indices at two times + */ +class SimplePathSegment : public std::pair +{ +public: + // Constructor. + SimplePathSegment( + int iTime0, + int iCandidate0, + int iTime1, + int iCandidate1) : std::pair( + TimeCandidatePair(iTime0, iCandidate0), + TimeCandidatePair(iTime1, iCandidate1)) {} +public: + // Comparator. + bool operator< (const SimplePathSegment & seg) const + { + return (first < seg.first); + } +}; + +typedef std::set SimplePathSegmentSet; + +// -------------------------------------------------------------------------- +/** + * A vector of times and candidate indices at those times defining a path + */ +class SimplePath +{ +public: + // Array of times. + std::vector m_iTimes; + + // Array of candidates. + std::vector m_iCandidates; +}; + +// A vector of SimplePaths. +typedef std::vector SimplePathVector; + +// -------------------------------------------------------------------------- +/** + * An operator which is used to enforce some criteria on paths + */ +class PathThresholdOp +{ +public: + // Index denoting all points along path. + static const int Count_All = (-1); + + // Index denoting only the first point along the path. + static const int Count_First = (-2); + + // Index denoting only the last point along the path. + static const int Count_Last = (-3); + + // Possible operations. + enum Operation { + GreaterThan, + LessThan, + GreaterThanEqualTo, + LessThanEqualTo, + EqualTo, + NotEqualTo, + AbsGreaterThanEqualTo, + AbsLessThanEqualTo + }; + +public: + // Parse a threshold operator string. + void Parse(const std::string &strOp, + const std::vector &vecFormatStrings) + { + // Read mode + enum { + ReadMode_Column, + ReadMode_Op, + ReadMode_Value, + ReadMode_MinCount, + ReadMode_Invalid + } eReadMode = ReadMode_Column; + + // Loop through string + int iLast = 0; + for (long unsigned int i = 0; i <= strOp.length(); i++) + { + // Comma-delineated + if ((i == strOp.length()) || (strOp[i] == ',')) + { + std::string strSubStr = strOp.substr(iLast, i - iLast); + + // Read in column name + if (eReadMode == ReadMode_Column) + { + long unsigned int j = 0; + for (; j < vecFormatStrings.size(); j++) + { + if (strSubStr == vecFormatStrings[j]) + { + m_iColumn = j; + break; + } + } + if (j == vecFormatStrings.size()) + { + TECA_ERROR("Threshold column name \"" << strSubStr << "\" " + << "not found in --format"); + } + + iLast = i + 1; + eReadMode = ReadMode_Op; + + // Read in operation + } + else if (eReadMode == ReadMode_Op) + { + if (strSubStr == ">") + { + m_eOp = GreaterThan; + } + else if (strSubStr == "<") + { + m_eOp = LessThan; + } + else if (strSubStr == ">=") + { + m_eOp = GreaterThanEqualTo; + } + else if (strSubStr == "<=") + { + m_eOp = LessThanEqualTo; + } + else if (strSubStr == "=") + { + m_eOp = EqualTo; + } + else if (strSubStr == "!=") + { + m_eOp = NotEqualTo; + } + else if (strSubStr == "|>=") + { + m_eOp = AbsGreaterThanEqualTo; + } + else if (strSubStr == "|<=") + { + m_eOp = AbsLessThanEqualTo; + } + else + { + TECA_ERROR("Threshold invalid operation \"" << strSubStr << "\" ") + } + + iLast = i + 1; + eReadMode = ReadMode_Value; + + // Read in value + } + else if (eReadMode == ReadMode_Value) + { + m_dValue = atof(strSubStr.c_str()); + iLast = i + 1; + eReadMode = ReadMode_MinCount; + + // Read in minimum count + } + else if (eReadMode == ReadMode_MinCount) + { + if (strSubStr == "all") + { + m_nMinimumCount = Count_All; + } + else if (strSubStr == "first") + { + m_nMinimumCount = Count_First; + } + else if (strSubStr == "last") + { + m_nMinimumCount = Count_Last; + } + else + { + m_nMinimumCount = atoi(strSubStr.c_str()); + } + + if (m_nMinimumCount < Count_Last) + { + TECA_ERROR("Invalid minimum count \"" + << m_nMinimumCount << "\""); + } + + iLast = i + 1; + eReadMode = ReadMode_Invalid; + + // Invalid + } + else if (eReadMode == ReadMode_Invalid) + { + TECA_ERROR("Too many entries in threshold string \"" + << strOp << "\""); + } + } + } + + if (eReadMode != ReadMode_Invalid) + { + TECA_ERROR("Insufficient entries in threshold string \"" + << strOp << "\""); + } + + // Output announcement + std::string strDescription; + strDescription += vecFormatStrings[m_iColumn]; + if (m_eOp == GreaterThan) + { + strDescription += " greater than "; + } + else if (m_eOp == LessThan) + { + strDescription += " less than "; + } + else if (m_eOp == GreaterThanEqualTo) + { + strDescription += " greater than or equal to "; + } + else if (m_eOp == LessThanEqualTo) + { + strDescription += " less than or equal to "; + } + else if (m_eOp == EqualTo) + { + strDescription += " equal to "; + } + else if (m_eOp == NotEqualTo) + { + strDescription += " not equal to "; + } + else if (m_eOp == AbsGreaterThanEqualTo) + { + strDescription += " magnitude greater than or equal to "; + } + else if (m_eOp == AbsLessThanEqualTo) + { + strDescription += " magnitude less than or equal to "; + } + + char szValue[128]; + snprintf(szValue, 128, "%f", m_dValue); + strDescription += szValue; + + char szMinCount[160]; + if (m_nMinimumCount == Count_All) + { + strDescription += " at all times"; + } + else if (m_nMinimumCount == Count_First) + { + strDescription += " at the first time"; + } + else if (m_nMinimumCount == Count_Last) + { + strDescription += " at the last time"; + } + else + { + snprintf(szMinCount, 160, " at least %i time(s)", m_nMinimumCount); + strDescription += szMinCount; + } + + Announce("%s", strDescription.c_str()); + } + + // Verify that the specified path satisfies the threshold op. + bool Apply(const SimplePath &path, + const std::vector &vecCandidates) + { + int nCount = 0; + for (long unsigned int s = 0; s < path.m_iTimes.size(); s++) + { + int t = path.m_iTimes[s]; + int i = path.m_iCandidates[s]; + + if ((m_nMinimumCount == Count_First) && (s > 0)) { + continue; + } + if ((m_nMinimumCount == Count_Last) && (s < path.m_iTimes.size()-1)) { + continue; + } + + _ASSERT((t >= 0) && (t < vecCandidates.size())); + + double dCandidateValue = + std::stod(vecCandidates[t]->second[i][m_iColumn]); + + if ((m_eOp == GreaterThan) && + (dCandidateValue > m_dValue)) + { + nCount++; + } + else if ((m_eOp == LessThan) && + (dCandidateValue < m_dValue)) + { + nCount++; + } + else if ((m_eOp == GreaterThanEqualTo) && + (dCandidateValue >= m_dValue)) + { + nCount++; + } + else if ((m_eOp == LessThanEqualTo) && + (dCandidateValue <= m_dValue)) + { + nCount++; + } + else if ((m_eOp == EqualTo) && + (dCandidateValue == m_dValue)) + { + nCount++; + } + else if ((m_eOp == NotEqualTo) && + (dCandidateValue != m_dValue)) + { + nCount++; + } + else if ((m_eOp == AbsGreaterThanEqualTo) && + (fabs(dCandidateValue) >= m_dValue)) + { + nCount++; + } + else if ((m_eOp == AbsLessThanEqualTo) && + (fabs(dCandidateValue) <= m_dValue)) + { + nCount++; + } + } + + // Check that the criteria is satisfied for all segments + if (m_nMinimumCount == Count_All) + { + if (nCount == (int)(path.m_iTimes.size())) + { + return true; + } + else + { + return false; + } + } + + // Check that the criteria are satisfied for the first or last segment + if ((m_nMinimumCount == Count_First) || (m_nMinimumCount == Count_Last)) + { + if (nCount == 1) + { + return true; + } + else + { + return false; + } + } + + // Check total count against min count + if (nCount >= m_nMinimumCount) + { + return true; + } + else + { + return false; + } + } + +protected: + // Active column. + int m_iColumn; + + // Operation. + Operation m_eOp; + + // Threshold value. + double m_dValue; + + // Minimum number of segments that need to satisfy the op. + int m_nMinimumCount; +}; + +// -------------------------------------------------------------------------- +void generate_path_segments_set_basic( + const std::vector &vecIterCandidates, + const std::vector> &vecNodes, + const std::vector &vecKDTrees, + int nMaxGapSteps, + double dMaxGapSeconds, + double dRangeDeg, + std::vector & vecPathSegmentsSet) +{ + vecPathSegmentsSet.resize(vecIterCandidates.size()-1); + + // Null pointer + int * noptr = NULL; + + // Search nodes at current time level + for (size_t t = 0; t < vecIterCandidates.size()-1; t++) + { + const Time & timeCurrent = vecIterCandidates[t]->first; + const TimesliceCandidateInformation & tscinfo = vecIterCandidates[t]->second; + + // Loop through all candidates at current time level + for (long unsigned int i = 0; i < tscinfo.size(); i++) + { + // Check future timesteps for next candidate in path + for (int g = 1; ; g++) + { + if ((nMaxGapSteps != (-1)) && (g > nMaxGapSteps+1)) + { + break; + } + if (t+g >= vecIterCandidates.size()) + { + break; + } + + if (dMaxGapSeconds >= 0.0) + { + const Time & timeNext = vecIterCandidates[t+g]->first; + double dDeltaSeconds = timeCurrent.DeltaSeconds(timeNext); + if (dDeltaSeconds > dMaxGapSeconds) + { + break; + } + } + + if (vecKDTrees[t+g] == NULL) + { + continue; + } + + kdres * set = kd_nearest3(vecKDTrees[t+g], vecNodes[t][i].x, vecNodes[t][i].y, vecNodes[t][i].z); + + if (kd_res_size(set) == 0) + { + kd_res_free(set); + break; + } + + int iRes = reinterpret_cast(kd_res_item_data(set)) + - reinterpret_cast(noptr); + + kd_res_free(set); + + // Great circle distance between points + double dRDeg = GreatCircleDistance_Deg(vecNodes[t][i].lonRad, + vecNodes[t][i].latRad, + vecNodes[t+g][iRes].lonRad, + vecNodes[t+g][iRes].latRad); + + // Verify great circle distance satisfies range requirement + if (dRDeg <= dRangeDeg) + { + // Insert new path segment into set of path segments + vecPathSegmentsSet[t].insert( + SimplePathSegment(t, i, t+g, iRes)); + + break; + } + } + } + } +} + +// -------------------------------------------------------------------------- +/** + * Priority method for generating path segments + */ +void generate_path_segments_with_priority( + const std::vector &vecIterCandidates, + const std::vector> &vecNodes, + const std::vector &vecKDTrees, + const int ixPriorityCol, + int nMaxGapSteps, + double dMaxGapSeconds, + double dRangeDeg, + std::vector &vecPathSegmentsSet) +{ + // For target candidates, priority includes both timestep delta and distance + typedef std::pair PriorityPair; + + // Null pointer + int * noptr = NULL; + + // Chord distance + const double dChordLength = ChordLengthFromGreatCircleDistance_Deg(dRangeDeg); + + // Initialize the path segment set + vecPathSegmentsSet.resize(vecIterCandidates.size()-1); + + // Get priority of all candidates + if (ixPriorityCol < 0) + { + TECA_ERROR("Invalid priority column index"); + } + std::multimap mapPriority; + for (size_t t = 0; t < vecIterCandidates.size()-1; t++) + { + const TimesliceCandidateInformation & tscinfo = vecIterCandidates[t]->second; + + for (long unsigned int i = 0; i < tscinfo.size(); i++) + { + if (tscinfo[i].size() <= ixPriorityCol) + { + TECA_ERROR("Priority column index out of range (" << tscinfo[i].size() << "<=" << ixPriorityCol << ")"); + } + + double dPriority = std::stod(tscinfo[i][ixPriorityCol]); + mapPriority.insert(std::pair(dPriority, TimeCandidatePair(t,i))); + } + } + + Announce("%lu total candidates found", mapPriority.size()); + + // Set of candidates that are already targets of another node + std::set setUsedCandidates; + + // Loop through all candidates in increasing priority + for (auto & it : mapPriority) + { + const long unsigned int t = it.second.timeix(); + const int i = it.second.candidate(); + + const Time & timeCurrent = vecIterCandidates[t]->first; + + // Find all candidates within prescribed distance + std::multimap, TimeCandidatePair> mmapTargets; + for (int g = 1; ; g++) + { + if ((nMaxGapSteps != (-1)) && (g > nMaxGapSteps+1)) + { + break; + } + if (t+g >= vecKDTrees.size()) + { + break; + } + if (dMaxGapSeconds >= 0.0) + { + const Time & timeNext = vecIterCandidates[t+g]->first; + double dDeltaSeconds = timeCurrent.DeltaSeconds(timeNext); + if (dDeltaSeconds > dMaxGapSeconds) + { + break; + } + } + if (vecKDTrees[t+g] == NULL) + { + continue; + } + + kdres * set = kd_nearest_range3(vecKDTrees[t+g], + vecNodes[t][i].x, + vecNodes[t][i].y, + vecNodes[t][i].z, + dChordLength); + if (set == NULL) + { + TECA_ERROR("Fatal exception in kd_nearest_range3"); + } + if (kd_res_size(set) == 0) + { + kd_res_free(set); + continue; + } + + do + { + int iRes = reinterpret_cast(kd_res_item_data(set)) + - reinterpret_cast(noptr); + + //printf("%i %i %i %i - %lu %lu\n", t, i, t+g, iRes); + + double dRDeg = GreatCircleDistance_Deg(vecNodes[t][i].lonRad, + vecNodes[t][i].latRad, + vecNodes[t+g][iRes].lonRad, + vecNodes[t+g][iRes].latRad); + mmapTargets.insert(std::pair(PriorityPair(g,dRDeg), + TimeCandidatePair(t+g,iRes))); + + } while (kd_res_next(set)); + + kd_res_free(set); + } + + // Find the shortest target node and insert a segment + for (auto & itTarget : mmapTargets) + { + auto itUsed = setUsedCandidates.find(itTarget.second); + if (itUsed == setUsedCandidates.end()) + { + setUsedCandidates.insert(itTarget.second); + + vecPathSegmentsSet[t].insert( + SimplePathSegment(t, i, itTarget.second.timeix(), itTarget.second.candidate())); + + break; + } + } + } +} + +// -------------------------------------------------------------------------- +/** + * First-come-first-serve method for generating paths (default). + */ +void generate_paths_basic( + const std::vector & vecIterCandidates, + std::vector & vecPathSegmentsSet, + std::vector & vecPaths) +{ + // Loop through all times + for (size_t t = 0; t < vecIterCandidates.size()-1; t++) + { + // Loop through all remaining segments at this time + while (vecPathSegmentsSet[t].size() > 0) + { + // Create a new path starting with this segment + SimplePath path; + + SimplePathSegmentSet::iterator iterSeg + = vecPathSegmentsSet[t].begin(); + + path.m_iTimes.push_back(iterSeg->first.timeix()); + path.m_iCandidates.push_back(iterSeg->first.candidate()); + + int tx = t; + + for (;;) + { + path.m_iTimes.push_back(iterSeg->second.timeix()); + path.m_iCandidates.push_back(iterSeg->second.candidate()); + + long unsigned int txnext = iterSeg->second.timeix(); + + if (txnext >= vecIterCandidates.size()-1) + { + vecPathSegmentsSet[tx].erase(iterSeg); + break; + } + + SimplePathSegment segFind( + iterSeg->second.timeix(), iterSeg->second.candidate(), 0, 0); + + vecPathSegmentsSet[tx].erase(iterSeg); + + iterSeg = vecPathSegmentsSet[txnext].find(segFind); + + if (iterSeg == vecPathSegmentsSet[txnext].end()) + { + break; + } + + tx = txnext; + } + + // Add path to array of paths + vecPaths.push_back(path); + } + } +} + +// -------------------------------------------------------------------------- +// PIMPL idiom hides internals +class teca_stitch_nodes::internals_t +{ +public: + internals_t() {} + ~internals_t() {} + + static + int parse_variable_list(std::string &cmd, + std::vector &vec); + template + static + int parse_threshold_list(std::vector &vec_format_strings, + std::string &cmd, + std::vector &vec); +public: + std::vector vec_format_strings; + std::vector vec_threshold_op; + + double min_time_seconds; + double max_gap_seconds; + int max_gap_steps; + long unsigned int min_time_steps; + Time::CalendarType caltype; +}; + +// -------------------------------------------------------------------------- +int teca_stitch_nodes::internals_t::parse_variable_list( + std::string &cmd, + std::vector &vec) +{ + long unsigned int iVarBegin = 0; + long unsigned int iVarCurrent = 0; + + // Parse variable name + for (;;) + { + if ((iVarCurrent >= cmd.length()) || + (cmd[iVarCurrent] == ',') || + (cmd[iVarCurrent] == ' ') || + (cmd[iVarCurrent] == '\t') || + (cmd[iVarCurrent] == '\n') || + (cmd[iVarCurrent] == '\r')) + { + if (iVarCurrent == iVarBegin) + { + if (iVarCurrent >= cmd.length()) + { + break; + } + + iVarCurrent++; + iVarBegin++; + continue; + } + vec.push_back(cmd.substr(iVarBegin, iVarCurrent - iVarBegin)); + iVarBegin = iVarCurrent + 1; + } + + iVarCurrent++; + } + return 0; +} + +// -------------------------------------------------------------------------- +template +int teca_stitch_nodes::internals_t::parse_threshold_list( + std::vector &vec_format_strings, + std::string &cmd, + std::vector &vec) +{ + int i_last = 0; + for (long unsigned int i = 0; i <= cmd.length(); ++i) + { + if ((i == cmd.length()) || + (cmd[i] == ';') || + (cmd[i] == ':')) + { + std::string strSubStr = cmd.substr(i_last, i - i_last); + + int i_next_op = (int)(vec.size()); + vec.resize(i_next_op + 1); + vec[i_next_op].Parse(strSubStr, vec_format_strings); + + i_last = i + 1; + } + } + return 0; +} + +// -------------------------------------------------------------------------- +teca_stitch_nodes::teca_stitch_nodes() : + in_connect(""), + in_fmt(""), + min_time("10"), + cal_type("standard"), + max_gap("3"), + track_threshold_cmd(""), + prioritize(""), + min_path_length(0), + range(8.0), + min_endpoint_distance(0.0), + min_path_distance(0.0), + allow_repeated_times(false) +{ + this->set_number_of_input_connections(1); + this->set_number_of_output_ports(1); + + this->internals = new teca_stitch_nodes::internals_t; +} + +// -------------------------------------------------------------------------- +teca_stitch_nodes::~teca_stitch_nodes() +{ + delete this->internals; +} + +#if defined(TECA_HAS_BOOST) +// -------------------------------------------------------------------------- +void teca_stitch_nodes::get_properties_description( + const std::string &prefix, options_description &global_opts) +{ + options_description opts("Options for " + + (prefix.empty()?"teca_stitch_nodes":prefix)); + + opts.add_options() + TECA_POPTS_GET(std::string, prefix, in_connect, + "Connectivity file that describes the unstructured grid") + TECA_POPTS_GET(std::string, prefix, in_fmt, + "A comma-separated list of names of the auxiliary columns" + " within the candidates output file.") + TECA_POPTS_GET(std::string, prefix, min_time, + "Minimum duration of path") + TECA_POPTS_GET(std::string, prefix, cal_type, + "Calendar type") + TECA_POPTS_GET(std::string, prefix, max_gap, + "Maximum time gap (in time steps or duration)") + TECA_POPTS_GET(std::string, prefix, track_threshold_cmd, + "Threshold commands for path [col,op,value,count;...]") + TECA_POPTS_GET(std::string, prefix, prioritize, + "Variable to use when prioritizing paths") + TECA_POPTS_GET(int, prefix, min_path_length, + "Minimum path length") + TECA_POPTS_GET(double, prefix, range, + "The maximum distance between candidates" + " along a path (in great-circle degrees") + TECA_POPTS_GET(double, prefix, min_endpoint_distance, + "Minimum distance between endpoints of path") + TECA_POPTS_GET(double, prefix, min_path_distance, + "Minimum accumulated great-circle distance" + " between nodes in a path (in degrees)") + TECA_POPTS_GET(bool, prefix, allow_repeated_times, + "Allow repeated times") + ; + + this->teca_algorithm::get_properties_description(prefix, opts); + + global_opts.add(opts); +} + +// -------------------------------------------------------------------------- +void teca_stitch_nodes::set_properties( + const std::string &prefix, variables_map &opts) +{ + this->teca_algorithm::set_properties(prefix, opts); + + TECA_POPTS_SET(opts, std::string, prefix, in_connect) + TECA_POPTS_SET(opts, std::string, prefix, in_fmt) + TECA_POPTS_SET(opts, std::string, prefix, min_time) + TECA_POPTS_SET(opts, std::string, prefix, cal_type) + TECA_POPTS_SET(opts, std::string, prefix, max_gap) + TECA_POPTS_SET(opts, std::string, prefix, track_threshold_cmd) + TECA_POPTS_SET(opts, std::string, prefix, prioritize) + TECA_POPTS_SET(opts, int, prefix, min_path_length) + TECA_POPTS_SET(opts, double, prefix, range) + TECA_POPTS_SET(opts, double, prefix, min_endpoint_distance) + TECA_POPTS_SET(opts, double, prefix, min_path_distance) + TECA_POPTS_SET(opts, bool, prefix, allow_repeated_times) +} +#endif + +// -------------------------------------------------------------------------- +int teca_stitch_nodes::initialize() +{ + // Parse minimum time + this->internals->min_time_steps = this->min_path_length; + this->internals->min_time_seconds = 0.0; + if (STLStringHelper::IsIntegerIndex(this->min_time)) + { + this->internals->min_time_steps = std::stoi(this->min_time); + if (this->internals->min_time_steps < 1) + { + TECA_ERROR("Invalid value of --min_time;" + " expected positive integer"); + return -1; + } + + } + else + { + Time timeMinTime; + timeMinTime.FromFormattedString(this->min_time); + this->internals->min_time_seconds = timeMinTime.AsSeconds(); + if (this->internals->min_time_seconds <= 0.0) + { + TECA_ERROR("Invalid value of --min_time;" + " expected positive integer"); + return -1; + } + } + + // Parse maximum gap + this->internals->max_gap_steps = (-1); + this->internals->max_gap_seconds = -1.0; + if (STLStringHelper::IsIntegerIndex(this->max_gap)) + { + this->internals->max_gap_steps = std::stoi(this->max_gap); + if (this->internals->max_gap_steps < 0) + { + TECA_ERROR("Invalid value of --max_gap;" + " expected nonnegative integer"); + return -1; + } + } + else + { + Time timeMaxGap; + timeMaxGap.FromFormattedString(this->max_gap); + this->internals->max_gap_seconds = timeMaxGap.AsSeconds(); + if (this->internals->max_gap_seconds < 0.0) + { + TECA_ERROR("Invalid value of --max_gap;" + " expected nonnegative integer"); + return -1; + } + } + + // Parse calendar type + this->internals->caltype = Time::CalendarTypeFromString(this->cal_type); + + if (this->internals->caltype == Time::CalendarNone) + { + if (this->internals->min_time_seconds > 0.0) + { + TECA_ERROR("A calendar type (--caltype) must be specified"); + return -1; + } + if (this->internals->max_gap_seconds >= 0.0) + { + TECA_ERROR("A calendar type (--caltype) must be specified"); + return -1; + } + } + + // Check for connectivity file + if (this->in_connect != "") + { + TECA_ERROR("Loading grid data from connectivity file" + "is not supported") + return -1; + } + + // Parse format string + if (internals_t::parse_variable_list(this->in_fmt, + this->internals->vec_format_strings)) + return -1; + + if (this->internals->vec_format_strings.size() == 0) + { + TECA_ERROR("No format specified"); + return -1; + } + + // Parse the threshold string + if (this->track_threshold_cmd != "") + { + if (internals_t::parse_threshold_list(this->internals->vec_format_strings, + this->track_threshold_cmd, + this->internals->vec_threshold_op)) + return -1; + } + return 0; +} + +// -------------------------------------------------------------------------- +std::vector teca_stitch_nodes::get_upstream_request( + unsigned int port, const std::vector &input_md, + const teca_metadata &request) +{ +#ifdef TECA_DEBUG + cerr << teca_parallel_id() + << "teca_stitch_nodes::get_upstream_request" << endl; +#endif + (void)port; + (void)input_md; + + std::vector up_reqs; + + teca_metadata req(request); + up_reqs.push_back(req); + + return up_reqs; +} + +// -------------------------------------------------------------------------- +const_p_teca_dataset teca_stitch_nodes::execute( + unsigned int port, const std::vector &input_data, + const teca_metadata &request) +{ +#ifdef TECA_DEBUG + cerr << teca_parallel_id() << "teca_stitch_nodes:execute" << endl; +#endif + (void)port; + (void)request; + + // get the input mesh + const_p_teca_table candidates = + std::dynamic_pointer_cast(input_data[0]); + + // in parallel only rank 0 is required to have data + int rank = 0; +#if defined(TECA_HAS_MPI) + int init = 0; + MPI_Initialized(&init); + if (init) + MPI_Comm_rank(this->get_communicator(), &rank); +#endif + if (!candidates) + { + if (rank == 0) + { + TECA_FATAL_ERROR("empty input or not a table") + } + return nullptr; + } + + // check that there are some candidates to work with. + unsigned long n_rows = candidates->get_number_of_rows(); + if (n_rows < 1) + { + TECA_FATAL_ERROR("Failed to form TC tracks because there were no candidates") + return nullptr; + } + + // Check format string for lat/lon + int iLatIndex = (-1); + int iLonIndex = (-1); + + // validate the input table + for (size_t i = 0; i < this->internals->vec_format_strings.size(); ++i) + { + if (!candidates->has_column(this->internals->vec_format_strings[i])) + { + TECA_FATAL_ERROR("Candidate table missing \"" + << this->internals->vec_format_strings[i] << "\"") + return nullptr; + } + + if (this->internals->vec_format_strings[i] == "lat") + { + iLatIndex = i; + } + if (this->internals->vec_format_strings[i] == "lon") + { + iLonIndex = i; + } + } + + if (iLatIndex == (-1)) + { + TECA_ERROR("Latitude \"lat\" must be specified in format (--in_fmt)"); + return nullptr; + } + if (iLonIndex == (-1)) + { + TECA_ERROR("Latitude \"lon\" must be specified in format (--in_fmt)"); + return nullptr; + } + + // Parse the input into mapCandidateInfo, then store iterators to all elements + // of mapCandidateInfo in vecIterCandidates to enable sequential access. + TimeToCandidateInfoMap mapCandidateInfo; + std::vector vecIterCandidates; + + AnnounceStartBlock("Loading candidate data"); + + if (parse_detect_nodes_file(candidates, this->internals->vec_format_strings, + mapCandidateInfo, this->internals->caltype, + this->allow_repeated_times)) + { + TECA_FATAL_ERROR("parse_detect_nodes_file encountered an error") + return nullptr; + } + + vecIterCandidates.reserve(mapCandidateInfo.size()); + for (auto it = mapCandidateInfo.begin(); it != mapCandidateInfo.end(); it++) + { + vecIterCandidates.push_back(it); + } + + const Time &timeFirst = vecIterCandidates[0]->first; + const Time &timeLast = vecIterCandidates[vecIterCandidates.size()-1]->first; + + // Verify total time is larger than mintime + if (this->internals->min_time_seconds > 0.0) + { + double delta_seconds = timeFirst.DeltaSeconds(timeLast); + if (this->internals->min_time_seconds > delta_seconds) + { + TECA_FATAL_ERROR("Duration of record " + << "is shorter than --min_time; no paths possible") + return nullptr; + } + } + + // Verify adjacent times are smaller than maxgap + if (this->internals->max_gap_seconds >= 0.0) + { + for (size_t t = 0; t < vecIterCandidates.size()-1; t++) + { + const Time &timeCurrent = vecIterCandidates[t]->first; + const Time &timeNext = vecIterCandidates[t+1]->first; + + double delta_seconds = timeCurrent.DeltaSeconds(timeNext); + if (delta_seconds > this->internals->max_gap_seconds) + { + TECA_FATAL_ERROR("Discrete times " << t << " and " << t+1 + << "differ by more than --max_gap") + return nullptr; + } + } + } + AnnounceEndBlock("Done"); + + AnnounceStartBlock("Creating KD trees at each time level"); + + // Null pointer + int * noptr = NULL; + + // Vector of lat/lon values + std::vector< std::vector > vecNodes; + vecNodes.resize(mapCandidateInfo.size()); + + // Vector of KD trees + std::vector vecKDTrees; + vecKDTrees.resize(mapCandidateInfo.size()); + + for (size_t t = 0; t < vecIterCandidates.size(); t++) + { + const TimesliceCandidateInformation & tscinfo = vecIterCandidates[t]->second; + + // Create a new kdtree + if (tscinfo.size() == 0) + { + vecKDTrees[t] = NULL; + continue; + } + + vecKDTrees[t] = kd_create(3); + + vecNodes[t].resize(tscinfo.size()); + + // Insert all points at this time level + for (long unsigned int i = 0; i < tscinfo.size(); i++) + { + vecNodes[t][i].lonRad = DegToRad(std::stod(tscinfo[i][iLonIndex])); + vecNodes[t][i].latRad = DegToRad(std::stod(tscinfo[i][iLatIndex])); + + RLLtoXYZ_Rad(vecNodes[t][i].lonRad, vecNodes[t][i].latRad, + vecNodes[t][i].x, vecNodes[t][i].y, vecNodes[t][i].z); + + kd_insert3(vecKDTrees[t], vecNodes[t][i].x, vecNodes[t][i].y, + vecNodes[t][i].z, reinterpret_cast(noptr+i)); + } + } + AnnounceEndBlock("Done"); + + AnnounceStartBlock("Constructing paths"); + + std::vector vecPathSegmentsSet; + + AnnounceStartBlock("Populating set of path segments"); + + // Vector over time of sets of all path segments that start at that time + if (this->prioritize == "") + { + generate_path_segments_set_basic(vecIterCandidates, + vecNodes, + vecKDTrees, + this->internals->max_gap_steps, + this->internals->max_gap_seconds, + this->range, + vecPathSegmentsSet); + } + else + { + int ixPriorityCol = (-1); + for (long unsigned int i = 0; i < this->internals->vec_format_strings.size(); i++) + { + if (this->internals->vec_format_strings[i] == this->prioritize) + { + ixPriorityCol = i; + break; + } + } + + if (ixPriorityCol == (-1)) + { + TECA_ERROR("Format string (--in_fmt) does not contain priority column \"" << this->prioritize << "\"") + } + + generate_path_segments_with_priority(vecIterCandidates, + vecNodes, + vecKDTrees, + ixPriorityCol, + this->internals->max_gap_steps, + this->internals->max_gap_seconds, + this->range, + vecPathSegmentsSet); + } + + AnnounceEndBlock("Done"); + + AnnounceStartBlock("Connecting path segments"); + + std::vector vecPaths; + + generate_paths_basic(vecIterCandidates, + vecPathSegmentsSet, + vecPaths); + + AnnounceEndBlock("Done"); + + AnnounceEndBlock("Done"); + + AnnounceStartBlock("Filtering paths"); + + int nRejectedMinTimePaths = 0; + int nRejectedMinEndpointDistPaths = 0; + int nRejectedMinPathDistPaths = 0; + int nRejectedThresholdPaths = 0; + + std::vector vecPathsUnfiltered = vecPaths; + vecPaths.clear(); + + // Loop through all paths + for (const SimplePath & path : vecPathsUnfiltered) + { + // Reject path due to minimum time + if (path.m_iTimes.size() < this->internals->min_time_steps) + { + nRejectedMinTimePaths++; + continue; + } + if (this->internals->min_time_steps > 0.0) + { + int nT = path.m_iTimes.size(); + int iTime0 = path.m_iTimes[0]; + int iTime1 = path.m_iTimes[nT-1]; + + const Time &timeFirst = vecIterCandidates[iTime0]->first; + const Time &timeLast = vecIterCandidates[iTime1]->first; + + double dDeltaSeconds = timeFirst.DeltaSeconds(timeLast); + if (dDeltaSeconds < this->internals->min_time_seconds) + { + nRejectedMinTimePaths++; + continue; + } + } + + // Reject path due to minimum endpoint distance + if (this->min_endpoint_distance > 0.0) + { + int nT = path.m_iTimes.size(); + + int iTime0 = path.m_iTimes[0]; + int iRes0 = path.m_iCandidates[0]; + + int iTime1 = path.m_iTimes[nT-1]; + int iRes1 = path.m_iCandidates[nT-1]; + + double dLonRad0 = vecNodes[iTime0][iRes0].lonRad; + double dLatRad0 = vecNodes[iTime0][iRes0].latRad; + + double dLonRad1 = vecNodes[iTime1][iRes1].lonRad; + double dLatRad1 = vecNodes[iTime1][iRes1].latRad; + + double dR = GreatCircleDistance_Deg(dLonRad0, dLatRad0, dLonRad1, dLatRad1); + + if (dR < this->min_endpoint_distance) + { + nRejectedMinEndpointDistPaths++; + continue; + } + } + + // Reject path due to minimum total path distance + if (this->min_path_distance > 0.0) + { + double dTotalPathDistance = 0.0; + for (long unsigned int i = 0; i < path.m_iTimes.size() - 1; i++) + { + int iTime0 = path.m_iTimes[i]; + int iRes0 = path.m_iCandidates[i]; + + int iTime1 = path.m_iTimes[i+1]; + int iRes1 = path.m_iCandidates[i+1]; + + double dLonRad0 = vecNodes[iTime0][iRes0].lonRad; + double dLatRad0 = vecNodes[iTime0][iRes0].latRad; + + double dLonRad1 = vecNodes[iTime1][iRes1].lonRad; + double dLatRad1 = vecNodes[iTime1][iRes1].latRad; + + double dR = GreatCircleDistance_Deg(dLonRad0, dLatRad0, dLonRad1, dLatRad1); + + dTotalPathDistance += dR; + } + + if (dTotalPathDistance < this->min_path_distance) + { + nRejectedMinPathDistPaths++; + continue; + } + } + + // Reject path due to threshold + bool fOpResult = true; + for (long unsigned int x = 0; x < this->internals->vec_threshold_op.size(); x++) + { + fOpResult = this->internals->vec_threshold_op[x].Apply(path, vecIterCandidates); + + if (!fOpResult) + { + break; + } + } + if (!fOpResult) + { + nRejectedThresholdPaths++; + continue; + } + + // Add path to array of paths + vecPaths.push_back(path); + } + + Announce("Paths rejected (mintime): %i", nRejectedMinTimePaths); + Announce("Paths rejected (minendpointdist): %i", nRejectedMinEndpointDistPaths); + Announce("Paths rejected (minpathdist): %i", nRejectedMinPathDistPaths); + Announce("Paths rejected (threshold): %i", nRejectedThresholdPaths); + Announce("Total paths found: %i", vecPaths.size()); + AnnounceEndBlock("Done"); + + AnnounceStartBlock("Writing results"); + + // create the table to hold storm tracks + p_teca_table storm_tracks = teca_table::New(); + storm_tracks->copy_metadata(candidates); + + std::string time_units; + storm_tracks->get_time_units(time_units); + + storm_tracks->declare_columns("storm_id", long(), "path_length", long(), + "year", int(), "month", int(), "day", int(), "hour", int()); + + // Copy over column data from candidate file to track file for output + n_rows = 0; + for (long unsigned int i = 0; i < vecPaths.size(); i++) + { + n_rows += vecPaths[i].m_iTimes.size(); + } + + + for (long unsigned int var = 0; var < this->internals->vec_format_strings.size(); var++) + { + int row = 0; + double *d_output_array = new double[n_rows]; + int *i_output_array = new int[n_rows]; + for (long unsigned int i = 0; i < vecPaths.size(); i++) + { + for (long unsigned int t = 0; t < vecPaths[i].m_iTimes.size(); t++) + { + int iTime = vecPaths[i].m_iTimes[t]; + + int iCandidate = vecPaths[i].m_iCandidates[t]; + + if (var == 0) + { + storm_tracks << i << vecPaths[i].m_iTimes.size() + << vecIterCandidates[iTime]->first.GetYear() + << vecIterCandidates[iTime]->first.GetMonth() + << vecIterCandidates[iTime]->first.GetDay() + << ceil(vecIterCandidates[iTime]->first.GetSecond() / 3600.); + } + + if (this->internals->vec_format_strings[var] == "i" || + this->internals->vec_format_strings[var] == "j") + i_output_array[row] = std::stoi(vecIterCandidates[iTime]->second[iCandidate][var]); + else + d_output_array[row] = std::stod(vecIterCandidates[iTime]->second[iCandidate][var]); + row++; + } + } + + + if (this->internals->vec_format_strings[var] == "i" || + this->internals->vec_format_strings[var] == "j") + { + p_teca_variant_array_impl output = + teca_variant_array_impl::New(n_rows, i_output_array); + storm_tracks->append_column(this->internals->vec_format_strings[var], output); + } + else + { + p_teca_variant_array_impl output = + teca_variant_array_impl::New(n_rows, d_output_array); + storm_tracks->append_column(this->internals->vec_format_strings[var], output); + } + + delete [] d_output_array; + delete [] i_output_array; + } + + AnnounceEndBlock("Done"); + + // Cleanup + AnnounceStartBlock("Cleanup"); + + for (long unsigned int t = 0; t < vecKDTrees.size(); t++) + { + kd_free(vecKDTrees[t]); + } + + AnnounceEndBlock("Done"); + + return storm_tracks; +} diff --git a/alg/teca_stitch_nodes.h b/alg/teca_stitch_nodes.h new file mode 100644 index 000000000..a4302a475 --- /dev/null +++ b/alg/teca_stitch_nodes.h @@ -0,0 +1,159 @@ +#ifndef teca_stitch_nodes_h +#define teca_stitch_nodes_h + +#include "teca_algorithm.h" +#include "teca_metadata.h" + +#include +#include + +TECA_SHARED_OBJECT_FORWARD_DECL(teca_stitch_nodes) + +/** + * Porting the StitchNodes function from TempestExtremes to TECA. + * stitch_nodes is used to connect nodal features together in time, + * producing paths associated with singular features. + * Additional filtering of the output of detect_nodes can be applied + * based on the temporal features of these paths. + */ +class TECA_EXPORT teca_stitch_nodes : public teca_algorithm +{ +public: + TECA_ALGORITHM_STATIC_NEW(teca_stitch_nodes) + TECA_ALGORITHM_DELETE_COPY_ASSIGN(teca_stitch_nodes) + TECA_ALGORITHM_CLASS_NAME(teca_stitch_nodes) + ~teca_stitch_nodes(); + + // report/initialize to/from Boost program options + // objects. + TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION() + TECA_SET_ALGORITHM_PROPERTIES() + + int initialize(); + + /** @name in_connet + * Set the connectivity file + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, in_connect) + ///@} + + /** @name in_fmt + * Tracks output commands + * [var,op,dist;...] + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, in_fmt) + ///@} + + /** @name min_time + * Minimum duration of path + * default 10 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, min_time) + ///@} + + /** @name cal_type + * Calendar type + * default standard + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, cal_type) + ///@} + + /** @name max_gap + * Maximum time gap (in time steps or duration) + * default 3 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, max_gap) + ///@} + + /** @name track_threshold_cmd + * Threshold commands for path + *[var,op,value,count;...] + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, track_threshold_cmd) + ///@} + + /** @name prioritize + * Variable to use when prioritizing paths + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, prioritize) + ///@} + + /** @name min_path_length + * Minimum path length + * default 0 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(int, min_path_length) + ///@} + + /** @name range + * Range (in degrees) + * default 8 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, range) + ///@} + + /** @name min_endpoint_distance + * Minimum distance between endpoints of path + * default 0 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, min_endpoint_distance) + ///@} + + /** @name min_path_distance + * Minimum path lengt + * default 0 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, min_path_distance) + ///@} + + /** @name allow_repeated_times + * Allow repeated times + * default false + */ + ///@{ + TECA_ALGORITHM_PROPERTY(bool, allow_repeated_times) + ///@} + +protected: + teca_stitch_nodes(); + +private: + std::vector get_upstream_request( + unsigned int port, + const std::vector &input_md, + const teca_metadata &request) override; + + const_p_teca_dataset execute( + unsigned int port, + const std::vector &input_data, + const teca_metadata &request) override; + +private: + std::string in_connect; + std::string in_fmt; + std::string min_time; + std::string cal_type; + std::string max_gap; + std::string track_threshold_cmd; + std::string prioritize; + int min_path_length; + double range; + double min_endpoint_distance; + double min_path_distance; + bool allow_repeated_times; + + class internals_t; + internals_t *internals; +}; +#endif diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt index 4bab5fdf7..a33a36114 100644 --- a/apps/CMakeLists.txt +++ b/apps/CMakeLists.txt @@ -16,6 +16,11 @@ if (TECA_HAS_BOOST) list(APPEND teca_app_link ${Boost_LIBRARIES}) endif() +if (TECA_HAS_TELITE) + include_directories(SYSTEM ${TELITE_INCLUDE_DIR}) + list(APPEND teca_app_link ${TELITE_LIBRARY}) +endif() + teca_add_app(teca_bayesian_ar_detect LIBS ${teca_app_link} FEATURES ${TECA_HAS_BOOST} ${TECA_HAS_NETCDF} ${TECA_HAS_UDUNITS}) @@ -53,6 +58,9 @@ teca_add_app(teca_regional_moisture_flux LIBS ${teca_app_link} FEATURES ${TECA_HAS_BOOST} ${TECA_HAS_NETCDF} ${TECA_HAS_UDUNITS} ${TECA_HAS_SHAPELIB}) +teca_add_app(teca_tempest_tc_detect LIBS ${teca_app_link} + FEATURES ${TECA_HAS_BOOST} ${TECA_HAS_TELITE} ${TECA_HAS_NETCDF}) + teca_add_python_app(teca_convert_table) teca_add_python_app(teca_dataset_metadata FEATURES ${TECA_HAS_NETCDF}) diff --git a/apps/teca_tempest_tc_detect.cpp b/apps/teca_tempest_tc_detect.cpp new file mode 100644 index 000000000..6343d6ed5 --- /dev/null +++ b/apps/teca_tempest_tc_detect.cpp @@ -0,0 +1,681 @@ +#include "teca_cf_reader.h" +#include "teca_multi_cf_reader.h" +#include "teca_normalize_coordinates.h" +#include "teca_l2_norm.h" +#include "teca_derived_quantity.h" +#include "teca_derived_quantity_numerics.h" +#include "teca_mpi_manager.h" +#include "teca_app_util.h" +#include "teca_table_reduce.h" +#include "teca_table_sort.h" +#include "teca_table_calendar.h" +#include "teca_table_writer.h" +#include "teca_cartesian_mesh_regrid.h" +#include "teca_cartesian_mesh_source.h" +#include "teca_rename_variables.h" +#include "teca_mesh_join.h" +#include "teca_detect_nodes.h" +#include "teca_stitch_nodes.h" + +#include +#include +#include + +using namespace std; +using namespace teca_derived_quantity_numerics; + +using boost::program_options::value; + +// -------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + // initialize mpi + teca_mpi_manager mpi_man(argc, argv); + + // initialize command line options description + // set up some common options to simplify use for most + // common scenarios + int help_width = 100; + options_description basic_opt_defs( + "Basic usage:\n\n" + "The following options are the most commonly used. Information\n" + "on advanced options can be displayed using --advanced_help\n\n" + "Basic command line options", help_width, help_width - 4 + ); + basic_opt_defs.add_options() + ("input_file", value(), "\na teca_multi_cf_reader configuration file" + " identifying the set of NetCDF CF2 files to process. When present data is" + " read using the teca_multi_cf_reader. Use one of either --input_file or" + " --input_regex.\n") + + ("input_regex", value(), "\na teca_cf_reader regex identifying the" + " set of NetCDF CF2 files to process. When present data is read using the" + " teca_cf_reader. Use one of either --input_file or --input_regex.\n") + + ("candidate_file", value()->default_value("candidates.csv"), + "\nfile path to write the storm candidates to. The extension determines" + " the file format. May be one of `.nc`, `.csv`, or `.bin`\n") + + ("track_file", value()->default_value("tracks.csv"), + "\nfile path to write the storm tracks to. The extension determines" + " the file format. May be one of `.nc`, `.csv`, or `.bin`\n") + + ("x_axis_variable", value()->default_value("lon"), + "\nName of the variable to use for x-coordinates.\n") + + ("y_axis_variable", value()->default_value("lat"), + "\nName of the variable to use for y-coordinates.\n") + + ("z_axis_variable", value()->default_value("level"), + "\nName of the variable to use for z-coordinates.\n") + + ("in_connect", value()->default_value(""), + "\nConnectivity file that describes the unstructured grid\n") + ("search_by_min", value()->default_value(""), + "\nThe input variable to use for initially selecting candidate points" + " (defined as local minima). At least one (and at most one)" + " of --search_by_min or --search_by_max must be specified.\n") + ("search_by_max", value()->default_value(""), + "\nThe input variable to use for initially selecting candidate points" + " (defined as local maxima). At least one (and at most one)" + " of --search_by_min or --search_by_max must be specified.\n") + ("closed_contour_cmd", value()->default_value(""), + "\nEliminate candidates if they do not have a closed contour." + " The closed contour is determined by breadth first search:" + " if any paths exist from the candidate point" + " (or nearby minima/maxima if minmaxdist is specified)" + " that reach the specified distance before achieving the specified delta" + " then we say no closed contour is present." + " Closed contour commands are separated by a semicolon (\";;...\")." + " Each closed contour command takes the form \"var,delta,dist,minmaxdist\"." + " These arguments are as follows." + " var is the name of the variable used for the contour search." + " dist is the great-circle distance (in degrees) from the pivot" + " within which the closed_contour criteria must be satisfied." + " delta is the amount by which the field must change from the pivot value." + " If positive (negative) the field must increase (decrease) by this value along the contour." + " minmaxdist is the great-circle distance away from the candidate" + " to search for the minima/maxima. If delta is positive (negative)," + " the pivot is a local minimum (maximum).\n") + ("no_closed_contour_cmd", value()->default_value(""), + "\nAs --closed_contour_cmd," + " except it eliminates candidates if a closed contour is present.\n") + ("candidate_threshold_cmd", value()->default_value(""), + "\nThreshold commands for candidates." + " Eliminate candidates that do not satisfy a threshold criteria" + " (there must exist a point within a given distance of the candidate" + " that satisfies a given equality or inequality)." + " Search is performed by breadth-first search over the grid." + " Threshold commands are separated by a semicolon (\";;...\")." + " Each threshold command takes the form \"var,op,value,dist\"." + " These arguments are as follows." + " var is the name of the variable used for the thresholding." + " op is the operator that must be satisfied for threshold (options include >,>=,<,<=,=,!=)." + " value is the value on the right-hand-side of the comparison." + " dist is the great-circle distance away from the candidate" + " to search for a point that satisfies the threshold.\n") + ("output_cmd", value()->default_value(""), + "\nInclude additional columns in the candidates output file." + " Each output command takes the form \"var,op,dist\"." + " These arguments are as follows." + " var is the name of the variable used for output." + " op is the operator that is applied over all points" + " within the specified distance of the candidate (options include max, min, avg, maxdist, mindist)." + " dist is the great-circle distance away from the candidate wherein the operator is applied.\n") + ("search_by_threshold", value()->default_value(""), + "\nThreshold for search operation in the form \"\"" + " These arguments are as follows." + " op is the operator that must be satisfied for threshold (options include >,>=,<,<=,=,!=)." + " value is the value on the right-hand-side of the comparison.\n") + ("min_lon", value()->default_value(0.0), + "\nMinimum longitude in degrees for detection\n") + ("max_lon", value()->default_value(0.0), + "\nMaximum longitude in degrees for detection" + " As longitude is a periodic dimension," + " when --regional is not specified --min_lon may be larger than --max_lon." + " If --max_lon and --min_lon are equal then these arguments are ignored.\n") + ("min_lat", value()->default_value(0.0), + "\nMinimum latitude in degrees for detection\n") + ("max_lat", value()->default_value(0.0), + "\nMaximum latitude in degrees for detection" + " If --max_lat and --min_lat are equal then these arguments are ignored.\n") + ("min_abs_lat", value()->default_value(0.0), + "\nMinimum absolute value of latitude in degrees for detection" + " This argument has no effect if set to zero.\n") + ("merge_dist", value()->default_value(6.0), + "\nMinimum allowable distance between two candidates in degrees." + " Candidate points with a distance (in degrees great-circle-distance)" + " shorter than the specified value are merged." + " Among two candidates within the merge distance," + " only the candidate with the lowest value of the --search_by_min field" + " or highest value of the --search_by_max field are retained.\n") + ("diag_connect", value()->default_value(false), + "\nWhen the data is on a structured grid," + " consider grid cells to be connected in the diagonal (across the vertex)\n") + ("regional", value()->default_value(true), + "\nUsed to indicate that a given latitude-longitude grid" + " should not be periodic in the longitudinal direction.\n") + ("out_header", value()->default_value(true), + "\nIf present, output a header at the beginning of the output file" + " indicating the columns of the file.\n") + + ("in_fmt", value()->default_value(""), + "\nA comma-separated list of names of the auxiliary columns" + " within the candidates output file." + " (namely, the list must not include the time columns).\n") + ("min_time", value()->default_value("10"), + "\nThe minimum length of a path either in terms of" + " number of discrete times or as a duration, e.g. \"24h\"." + " Note that the duration of a path is computed as" + " the difference between the final time and initial time," + " so a \"24h\" duration correspond to 5 time steps in 6-hourly data" + " (i.e. 0h,6,12,18,24UTC).\n") + ("cal_type", value()->default_value("standard"), + "\nCalendar type\n") + ("max_gap", value()->default_value("3"), + "\nMaximum time gap (in time steps or duration)" + " The number of allowed missing points between spatially proximal candidate nodes" + " while still considering them part of the same path.\n") + ("track_threshold_cmd", value()->default_value(""), + "\nThreshold commands for path" + " Filter paths based on the number of times where a particular threshold is satisfied." + " Threshold commands are separated by a semicolon (\";;...\")." + " Each threshold command takes the form \"col,op,value,count\"." + " These arguments are as follows." + " col is the name of the column to use for thresholding, as specified in --in_fmt." + " op is the operator that must be satisfied for threshold (options include >,>=,<,<=,=,!=,|>=,|<=)." + " value is the value on the right-hand-side of the comparison." + " count is either the minimum number of time slices where the threshold must be satisfied" + " or the instruction \"all\", \"first\", or \"last\"." + " Here \"all\" is used to indicate the threshold must be satisfied at all points along the path," + " \"first\" is used to indicate the threshold must be satisfied only at the first point along the path, and" + " \"last\" is used to indicate the threshold must be satisfied only at the last point along the path.\n") + ("prioritize", value()->default_value(""), + "\nVariable to use when prioritizing paths\n") + ("min_path_length", value()->default_value(0), + "\nMinimum path length\n") + ("range", value()->default_value(8.0), + "\nThe maximum distance between candidates" + " along a path (in great-circle degrees).\n") + ("min_endpoint_distance", value()->default_value(0.0), + "\nThe minimum great-circle distance between the first candidate" + " on a path and the last candidate (in degrees).\n") + ("min_path_distance", value()->default_value(0.0), + "\nThe minimum accumulated great-circle distance between nodes in a path (in degrees).\n") + ("allow_repeated_times", value()->default_value(false), + "\nAllow repeated times\n") + + ("sea_level_pressure", value()->default_value(""), + "\nname of variable with sea level pressure\n") + ("geopotential_at_surface", value()->default_value(""), + "\nname of variable with geopotential at the surface\n") + ("geopotential", value()->default_value(""), + "\nname of variable with geopotential (or geopotential height)" + " for thickness calc\n") + ("500mb_height", value()->default_value(""), + "\nname of variable with 500mb height for thickness calc\n") + ("300mb_height", value()->default_value(""), + "\nname of variable with 300mb height for thickness calc\n") + ("surface_wind_u", value()->default_value(""), + "\nname of variable with surface wind x-component\n") + ("surface_wind_v", value()->default_value(""), + "\nname of variable with surface wind y-component\n") + + ("first_step", value()->default_value(0), + "\nfirst time step to process\n") + ("last_step", value()->default_value(-1), + "\nlast time step to process\n") + ("n_threads", value()->default_value(-1), + "\nSets the thread pool size on each MPI rank." + " When the default value of -1 is used TECA will coordinate the thread" + " pools across ranks such each thread is bound to a unique physical core.\n") + + ("verbose", value()->default_value(0), + "\nUse 1 to enable verbose mode, otherwise 0.\n") + ("help", "\ndisplays documentation for application specific command line options\n") + ("advanced_help", "\ndisplays documentation for algorithm specific command line options\n") + ("full_help", "\ndisplays both basic and advanced documentation together\n") + ; + + // add all options from each pipeline stage for more advanced use + options_description advanced_opt_defs( + "Advanced usage:\n\n" + "The following list contains the full set options giving one full\n" + "control over all runtime modifiable parameters. The basic options\n" + "(see" "--help) map to these, and will override them if both are\n" + "specified.\n\n" + "Advanced command line options", help_width, help_width - 4 + ); + + // create the pipeline stages here, they contain the + // documentation and parse command line. + // objects report all of their properties directly + // set default options here so that command line options override + // them. while we are at it connect the pipeline + p_teca_cf_reader cf_reader = teca_cf_reader::New(); + p_teca_multi_cf_reader mcf_reader = teca_multi_cf_reader::New(); + p_teca_normalize_coordinates sim_coords = teca_normalize_coordinates::New(); + p_teca_cartesian_mesh_source regrid_src_1 = teca_cartesian_mesh_source::New(); + p_teca_cartesian_mesh_source regrid_src_2 = teca_cartesian_mesh_source::New(); + p_teca_cartesian_mesh_regrid regrid_1 = teca_cartesian_mesh_regrid::New(); + p_teca_cartesian_mesh_regrid regrid_2 = teca_cartesian_mesh_regrid::New(); + p_teca_rename_variables rename_1 = teca_rename_variables::New(); + p_teca_rename_variables rename_2 = teca_rename_variables::New(); + p_teca_mesh_join join = teca_mesh_join::New(); + p_teca_derived_quantity thickness = teca_derived_quantity::New(); + p_teca_l2_norm surf_wind = teca_l2_norm::New(); + p_teca_detect_nodes candidates = teca_detect_nodes::New(); + p_teca_table_reduce map_reduce = teca_table_reduce::New(); + p_teca_table_sort sort = teca_table_sort::New(); + p_teca_table_writer candidate_writer = teca_table_writer::New(); + p_teca_stitch_nodes tracks = teca_stitch_nodes::New(); + p_teca_table_calendar calendar = teca_table_calendar::New(); + p_teca_table_writer track_writer = teca_table_writer::New(); + + thickness->set_dependent_variables({"Z500", "Z300"}); + thickness->set_derived_variable("thickness"); + surf_wind->set_component_0_variable("VAR_10U"); + surf_wind->set_component_1_variable("VAR_10V"); + surf_wind->set_l2_norm_variable("surface_wind"); + + cf_reader->get_properties_description("cf_reader", advanced_opt_defs); + mcf_reader->get_properties_description("mcf_reader", advanced_opt_defs); + regrid_src_1->get_properties_description("regrid_source_1", advanced_opt_defs); + regrid_src_2->get_properties_description("regrid_source_2", advanced_opt_defs); + regrid_1->get_properties_description("regrid_1", advanced_opt_defs); + regrid_2->get_properties_description("regrid_2", advanced_opt_defs); + rename_1->get_properties_description("rename_1", advanced_opt_defs); + rename_2->get_properties_description("rename_2", advanced_opt_defs); + join->get_properties_description("join", advanced_opt_defs); + thickness->get_properties_description("thickness", advanced_opt_defs); + surf_wind->get_properties_description("surface_wind_speed", advanced_opt_defs); + candidates->get_properties_description("candidates", advanced_opt_defs); + map_reduce->get_properties_description("map_reduce", advanced_opt_defs); + sort->get_properties_description("sort", advanced_opt_defs); + candidate_writer->get_properties_description("candidate_writer", advanced_opt_defs); + tracks->get_properties_description("tracks", advanced_opt_defs); + calendar->get_properties_description("calendar", advanced_opt_defs); + track_writer->get_properties_description("track_writer", advanced_opt_defs); + + // package basic and advanced options for display + options_description all_opt_defs(help_width, help_width - 4); + all_opt_defs.add(basic_opt_defs).add(advanced_opt_defs); + + // parse the command line + int ierr = 0; + variables_map opt_vals; + if ((ierr = teca_app_util::process_command_line_help( + mpi_man.get_comm_rank(), argc, argv, basic_opt_defs, + advanced_opt_defs, all_opt_defs, opt_vals))) + { + if (ierr == 1) + return 0; + return -1; + } + + // pass command line arguments into the pipeline objects + // advanced options are processed first, so that the basic + // options will override them + cf_reader->set_properties("cf_reader", opt_vals); + mcf_reader->set_properties("mcf_reader", opt_vals); + regrid_src_1->set_properties("regrid_source_1", opt_vals); + regrid_src_2->set_properties("regrid_source_2", opt_vals); + regrid_1->set_properties("regrid_1", opt_vals); + regrid_2->set_properties("regrid_2", opt_vals); + rename_1->set_properties("rename_1", opt_vals); + rename_2->set_properties("rename_2", opt_vals); + join->set_properties("join", opt_vals); + thickness->set_properties("thickness", opt_vals); + surf_wind->set_properties("surface_wind_speed", opt_vals); + candidates->set_properties("candidates", opt_vals); + map_reduce->set_properties("map_reduce", opt_vals); + sort->set_properties("sort", opt_vals); + candidate_writer->set_properties("candidate_writer", opt_vals); + tracks->set_properties("tracks", opt_vals); + calendar->set_properties("calendar", opt_vals); + track_writer->set_properties("track_writer", opt_vals); + + // now pass in the basic options, these are processed + // last so that they will take precedence + + // configure the reader + bool have_file = opt_vals.count("input_file"); + bool have_regex = opt_vals.count("input_regex"); + // some minimal check for missing options + if ((have_file && have_regex) || !(have_file || have_regex)) + { + if (mpi_man.get_comm_rank() == 0) + { + TECA_FATAL_ERROR("Exactly one of --input_file or --input_regex can be specified. " + "Use --input_file to activate the multi_cf_reader (HighResMIP datasets) " + "and --input_regex to activate the cf_reader (CAM like datasets)") + } + return -1; + } + + if (!opt_vals["x_axis_variable"].defaulted()) + { + cf_reader->set_x_axis_variable(opt_vals["x_axis_variable"].as()); + mcf_reader->set_x_axis_variable(opt_vals["x_axis_variable"].as()); + } + + if (!opt_vals["y_axis_variable"].defaulted()) + { + cf_reader->set_y_axis_variable(opt_vals["y_axis_variable"].as()); + mcf_reader->set_y_axis_variable(opt_vals["y_axis_variable"].as()); + } + + if (!opt_vals["z_axis_variable"].defaulted()) + { + cf_reader->set_z_axis_variable(opt_vals["z_axis_variable"].as()); + mcf_reader->set_z_axis_variable(opt_vals["z_axis_variable"].as()); + } + + if (opt_vals.count("input_file")) + { + mcf_reader->set_input_file(opt_vals["input_file"].as()); + mcf_reader->set_validate_time_axis(0); + sim_coords->set_input_connection(mcf_reader->get_output_port()); + } + else if (opt_vals.count("input_regex")) + { + cf_reader->set_files_regex(opt_vals["input_regex"].as()); + sim_coords->set_input_connection(cf_reader->get_output_port()); + } + p_teca_algorithm head = sim_coords; + + if (!opt_vals["in_connect"].defaulted()) + { + candidates->set_in_connect(opt_vals["in_connect"].as()); + tracks->set_in_connect(opt_vals["in_connect"].as()); + } + + if (!opt_vals["search_by_min"].defaulted()) + { + candidates->set_search_by_min(opt_vals["search_by_min"].as()); + } + + if (!opt_vals["search_by_max"].defaulted()) + { + candidates->set_search_by_max(opt_vals["search_by_max"].as()); + } + + if (opt_vals["search_by_min"].as() == "" && opt_vals["search_by_max"].as() == "") + { + if (opt_vals["sea_level_pressure"].as() == "") + TECA_FATAL_ERROR("Missing name of variable with sea level pressure") + else + candidates->set_search_by_min(opt_vals["sea_level_pressure"].as()); + } + + if (!opt_vals["closed_contour_cmd"].defaulted()) + { + candidates->set_closed_contour_cmd(opt_vals["closed_contour_cmd"].as()); + } + else + { + if (opt_vals["geopotential"].as() == "") + { + if (opt_vals["300mb_height"].as() == "") + { + TECA_FATAL_ERROR("Missing name of variable with 500mb height" + " or with geopotential for thickness calc") + } + else + { + thickness->set_dependent_variable(1, opt_vals["300mb_height"].as()); + } + if (opt_vals["500mb_height"].as() == "") + { + TECA_FATAL_ERROR("Missing name of variable with 300mb height" + " or with geopotential for thickness calc") + } + else + { + thickness->set_dependent_variable(0, opt_vals["500mb_height"].as()); + } + } + else + { + teca_metadata md = sim_coords->update_metadata(); + teca_metadata coords; + md.get("coordinates", coords); + const_p_teca_variant_array x = coords.get("x"); + const_p_teca_variant_array y = coords.get("y"); + double bounds[6] = {0.0}; + md.get("bounds", bounds, 6); + + //slice variables at pressure level of 300mb + regrid_src_1->set_bounds({bounds[0], bounds[1], bounds[2], bounds[3], 300, 300, 0, 0}); + regrid_src_1->set_whole_extents({0lu, x->size() - 1lu, 0lu, y->size() - 1lu, 0lu, 0lu, 0lu, 0lu}); + regrid_src_1->set_t_axis_variable(md); + regrid_src_1->set_t_axis(md); + + regrid_1->set_interpolation_mode_linear(); + regrid_1->set_input_connection(0, regrid_src_1->get_output_port()); + regrid_1->set_input_connection(1, sim_coords->get_output_port()); + + rename_1->set_original_variable_names({opt_vals["geopotential"].as()}); + rename_1->set_new_variable_names({"Z300"}); + rename_1->set_input_connection(regrid_1->get_output_port()); + + //slice variables at pressure level of 500mb + regrid_src_2->set_bounds({bounds[0], bounds[1], bounds[2], bounds[3], 500, 500, 0, 0}); + regrid_src_2->set_whole_extents({0lu, x->size() - 1lu, 0lu, y->size() - 1lu, 0lu, 0lu, 0lu, 0lu}); + regrid_src_2->set_t_axis_variable(md); + regrid_src_2->set_t_axis(md); + + regrid_2->set_interpolation_mode_linear(); + regrid_2->set_input_connection(0, regrid_src_2->get_output_port()); + regrid_2->set_input_connection(1, sim_coords->get_output_port()); + + rename_2->set_original_variable_names({opt_vals["geopotential"].as()}); + rename_2->set_new_variable_names({"Z500"}); + rename_2->set_input_connection(regrid_2->get_output_port()); + + //join the two new meshes + join->set_number_of_input_connections(2); + join->set_input_connection(0, rename_1->get_output_port()); + join->set_input_connection(1, rename_2->get_output_port()); + head = join; + } + + size_t n_var = thickness->get_number_of_dependent_variables(); + if (n_var != 2) + { + TECA_FATAL_ERROR("thickness calculation requires 2 " + "variables. given " << n_var) + return -1; + } + thickness->set_execute_callback(point_wise_difference(thickness->get_dependent_variable(0), + thickness->get_dependent_variable(1), + thickness->get_derived_variable(0))); + std::string text = opt_vals["sea_level_pressure"].as()+",200.0,5.5,0;"+thickness->get_derived_variable(0)+",-6.0,6.5,1.0"; + candidates->set_closed_contour_cmd(text); + } + + if (!opt_vals["no_closed_contour_cmd"].defaulted()) + { + candidates->set_no_closed_contour_cmd(opt_vals["no_closed_contour_cmd"].as()); + } + + if (!opt_vals["candidate_threshold_cmd"].defaulted()) + { + candidates->set_candidate_threshold_cmd(opt_vals["candidate_threshold_cmd"].as()); + } + + if (!opt_vals["output_cmd"].defaulted()) + { + candidates->set_output_cmd(opt_vals["output_cmd"].as()); + } + else + { + if (opt_vals["surface_wind_u"].as() == "") + TECA_FATAL_ERROR("Missing name of variable with surface wind x-component") + else + surf_wind->set_component_0_variable(opt_vals["surface_wind_u"].as()); + + if (opt_vals["surface_wind_v"].as() == "") + TECA_FATAL_ERROR("Missing name of variable with surface wind y-component") + else + surf_wind->set_component_1_variable(opt_vals["surface_wind_v"].as()); + + if (opt_vals["geopotential_at_surface"].as() == "") + TECA_FATAL_ERROR("Missing name of variable with geopotential at the surface") + + std::string text = opt_vals["sea_level_pressure"].as()+",min,0;"+surf_wind->get_l2_norm_variable()+",max,2;"+opt_vals["geopotential_at_surface"].as()+",min,0"; + candidates->set_output_cmd(text); + } + if (!opt_vals["in_fmt"].defaulted()) + { + tracks->set_in_fmt(opt_vals["in_fmt"].as()); + } + else + { + std::string text = "i,j,lat,lon,"+opt_vals["sea_level_pressure"].as()+","+surf_wind->get_l2_norm_variable()+","+opt_vals["geopotential_at_surface"].as(); + tracks->set_in_fmt(text); + } + + if (!opt_vals["search_by_threshold"].defaulted()) + { + candidates->set_search_by_threshold(opt_vals["search_by_threshold"].as()); + } + + if (!opt_vals["min_lon"].defaulted()) + { + candidates->set_min_lon(opt_vals["min_lon"].as()); + } + + if (!opt_vals["max_lon"].defaulted()) + { + candidates->set_max_lon(opt_vals["max_lon"].as()); + } + + if (!opt_vals["min_lat"].defaulted()) + { + candidates->set_min_lat(opt_vals["min_lat"].as()); + } + + if (!opt_vals["max_lat"].defaulted()) + { + candidates->set_max_lat(opt_vals["max_lat"].as()); + } + + if (!opt_vals["min_abs_lat"].defaulted()) + { + candidates->set_min_abs_lat(opt_vals["min_abs_lat"].as()); + } + + if (!opt_vals["merge_dist"].defaulted()) + { + candidates->set_merge_dist(opt_vals["merge_dist"].as()); + } + + if (!opt_vals["diag_connect"].defaulted()) + { + candidates->set_diag_connect(opt_vals["diag_connect"].as()); + } + + if (!opt_vals["regional"].defaulted()) + { + candidates->set_regional(opt_vals["regional"].as()); + } + + if (!opt_vals["out_header"].defaulted()) + { + candidates->set_out_header(opt_vals["out_header"].as()); + } + candidates->set_verbose(opt_vals["verbose"].as()); + candidates->initialize(); + + if (!opt_vals["first_step"].defaulted()) + map_reduce->set_start_index(opt_vals["first_step"].as()); + + if (!opt_vals["last_step"].defaulted()) + map_reduce->set_end_index(opt_vals["last_step"].as()); + + if (!opt_vals["n_threads"].defaulted()) + map_reduce->set_thread_pool_size(opt_vals["n_threads"].as()); + else + map_reduce->set_thread_pool_size(-1); + + candidate_writer->set_file_name(opt_vals["candidate_file"].as()); + candidate_writer->set_output_format_auto(); + + sort->set_index_column("step"); + + if (!opt_vals["min_time"].defaulted()) + { + tracks->set_min_time(opt_vals["min_time"].as()); + } + + if (!opt_vals["cal_type"].defaulted()) + { + tracks->set_cal_type(opt_vals["cal_type"].as()); + } + + if (!opt_vals["max_gap"].defaulted()) + { + tracks->set_max_gap(opt_vals["max_gap"].as()); + } + + if (!opt_vals["track_threshold_cmd"].defaulted()) + { + tracks->set_track_threshold_cmd(opt_vals["track_threshold_cmd"].as()); + } + else + { + std::string text = surf_wind->get_l2_norm_variable()+",>=,10.0,10;lat,<=,50.0,10;lat,>=,-50.0,10;"+opt_vals["geopotential_at_surface"].as()+",<=,15.0,10"; + tracks->set_track_threshold_cmd(text); + } + + if (!opt_vals["prioritize"].defaulted()) + { + tracks->set_prioritize(opt_vals["prioritize"].as()); + } + + if (!opt_vals["min_path_length"].defaulted()) + { + tracks->set_min_path_length(opt_vals["min_path_length"].as()); + } + + if (!opt_vals["range"].defaulted()) + { + tracks->set_range(opt_vals["range"].as()); + } + + if (!opt_vals["min_endpoint_distance"].defaulted()) + { + tracks->set_min_endpoint_distance(opt_vals["min_endpoint_distance"].as()); + } + + if (!opt_vals["min_path_distance"].defaulted()) + { + tracks->set_min_path_distance(opt_vals["min_path_distance"].as()); + } + + if (!opt_vals["allow_repeated_times"].defaulted()) + { + tracks->set_allow_repeated_times(opt_vals["allow_repeated_times"].as()); + } + tracks->initialize(); + + track_writer->set_file_name(opt_vals["track_file"].as()); + track_writer->set_output_format_auto(); + + // connect all the stages + thickness->set_input_connection(head->get_output_port()); + surf_wind->set_input_connection(thickness->get_output_port()); + candidates->set_input_connection(surf_wind->get_output_port()); + map_reduce->set_input_connection(candidates->get_output_port()); + candidate_writer->set_input_connection(map_reduce->get_output_port()); + sort->set_input_connection(candidate_writer->get_output_port()); + calendar->set_input_connection(sort->get_output_port()); + tracks->set_input_connection(calendar->get_output_port()); + track_writer->set_input_connection(tracks->get_output_port()); + + // run the pipeline + track_writer->update(); + + return 0; +} diff --git a/doc/rtd/applications.rst b/doc/rtd/applications.rst index e35879d7d..afb66d561 100644 --- a/doc/rtd/applications.rst +++ b/doc/rtd/applications.rst @@ -30,6 +30,8 @@ batch script rather than in your shell. | :ref:`teca_temporal_reduction` | Computes reductions (min, max, average, | | | summation) over the time dimension | +----------------------------------------+--------------------------------------------------+ +| :ref:`teca_tempest_tc_detect` | TC detection using TempestExtremes | ++----------------------------------------+--------------------------------------------------+ | :ref:`teca_tc_detect` | TC detection using GFDL algorithm | +----------------------------------------+--------------------------------------------------+ | :ref:`teca_tc_trajectory` | Computes TC tracks from a set of candidates | @@ -57,6 +59,8 @@ batch script rather than in your shell. | :ref:`teca_cf_restripe` | Convert the internal layout of a dataset on disk | | | with optional subsetting and/or regridding. | +----------------------------------------+--------------------------------------------------+ +| :ref:`teca_lapse_rate` | Calculates the mean lapse rate | ++----------------------------------------+--------------------------------------------------+ Applying the Command Line Applications at Scale ----------------------------------------------- @@ -1557,6 +1561,310 @@ Command Line Arguments displays both basic and advanced documentation together +.. _teca_tempest_tc_detect: + +teca_tempest_tc_detect +------------------------------ +The cyclone detecton algorithm is based on two TempestExtremes codes :cite:`tempestextremes`: +DetectNodes and StitchNodes. + +Inputs +~~~~~~ +A Cartesian mesh stored in a collection of NetCDF CF2 files. The detector requires on +the following fields. + +1. Sea level pressure +2. Surface wind x-component +3. Surface wind y-component +4. Geopotential height (300 and 500 mB) +5. Surface geopotential height + +Outputs +~~~~~~~ +1. Cyclone candidate table +2. Cyclone track table + +Command Line Arguments +~~~~~~~~~~~~~~~~~~~~~~ + +--input_file INPUT_FILE + a teca_multi_cf_reader configuration file identifying the set of NetCDF CF2 files to process. + When present data is read using the teca_multi_cf_reader. Use one of either `--input_file` or + `--input_regex`. + +--input_regex INPUT_REGEX + a teca_cf_reader regex identifying the set of NetCDF CF2 files to process. When present data is + read using the teca_cf_reader. Use one of either `--input_file` or `--input_regex`. + +--candidate_file CANDIDATE_FILE + file path to write the storm candidates to. The extension determines the file format. May be one of + `.nc`, `.csv`, or `.bin`. (default: candidates.csv) + +--track_file TRACK_FILE + file path to write the storm tracks to. The extension determines the file format. May be one of + `.nc`, `.csv`, or `.bin`. (default: tracks.csv) + +--x_axis_variable X_AXIS_VARIABLE + name of x coordinate variable (default: lon) + +--y_axis_variable Y_AXIS_VARIABLE + name of y coordinate variable (default: lat) + +--z_axis_variable Z_AXIS_VARIABLE + name of z coordinate variable (default: level) + +--sea_level_pressure SEA_LEVEL_PRESSURE + name of variable with sea level pressure + +--surface_wind_u SURFACE_WIND_U + name of variable with surface wind x-component + +--surface_wind_v SURFACE_WIND_Y + name of variable with surface wind y-component + +--geopotential_at_surface GEOPOTENTIAL_AT_SURFACE + name of variable with geopotential at the surface + +--geopotential GEOPOTENTIAL + name of variable with geopotential height + for thickness calc (300 and 500 mB) + +--300mb_height 300_HEIGHT + name of variable with 300mb height for thickness calc + +--500mb_height 500_HEIGHT + name of variable with 500mb height for thickness calc + +--in_connect CONNECTIVITY_FILE + a connectivity file that describes the unstructured grid. + Not supported yet. + +--diag_connect true/false + When the data is on a structured grid, consider grid cells + to be connected in the diagonal (across the vertex). + (default: false) + +--regional true/false + Used to indicate that a given latitude-longitude grid + should not be periodic in the longitudinal direction. + (default: true) + +--out_header true/false + If present, output a header at the beginning of the output file + indicating the columns of the file. + (default: true) + +--search_by_min VARIABLE_NAME + The input variable to use for initially selecting candidate points (defined as local minima). + At least one (and at most one) of --search_by_min or --search_by_max must be specified. + +--search_by_max VARIABLE_NAME + The input variable to use for initially selecting candidate points (defined as local maxima). + At least one (and at most one) of --search_by_min or --search_by_max must be specified. + +--search_by_threshold SEARCH_BY_THRESHOLD + Threshold for search operation in the form "" + These arguments are as follows. + op is the operator that must be satisfied for threshold (options include >,>=,<,<=,=,!=). + value is the value on the right-hand-side of the comparison. + +--min_lon MIN_LON + The minimum longitude for candidate points. (default: 0.0) + +--max_lon MAX_LON + The maximum longitude for candidate points. + As longitude is a periodic dimension, + when --regional is not specified --min_lon may be larger than --max_lon. + If --max_lon and --min_lon are equal then these arguments are ignored. + (default: 0.0) + +--min_lat MIN_LAT + The minimum latitude for candidate points. (default: 0.0) + +--max_lat MAX_LAT + The maximum latitude for candidate points. + If --max_lat and --min_lat are equal then these arguments are ignored. + (default: 0.0) + +--min_abs_lat MIN_ABS_LAT + The minimum absolute value of the latitude for candidate points. + This argument has no effect if set to zero. + (default: 0.0) + +--merge_dist MERGE_DIST + Minimum allowable distance between two candidates in degrees. + Candidate points with a distance (in degrees great-circle-distance) + shorter than the specified value are merged. + Among two candidates within the merge distance, + only the candidate with the lowest value of the --search_by_min field + or highest value of the --search_by_max field are retained. + (default: 6.0) + +--closed_contour_cmd CLOSED_CONTOUR_CMD + Eliminate candidates if they do not have a closed contour. + The closed contour is determined by breadth first search: + if any paths exist from the candidate point + (or nearby minima/maxima if minmaxdist is specified) + that reach the specified distance before achieving the specified delta + then we say no closed contour is present. + Closed contour commands are separated by a semicolon (";;..."). + Each closed contour command takes the form "var,delta,dist,minmaxdist". + These arguments are as follows. + var is the name of the variable used for the contour search. + dist is the great-circle distance (in degrees) from the pivot + within which the closed_contour criteria must be satisfied. + delta is the amount by which the field must change from the pivot value. + If positive (negative) the field must increase (decrease) by this value along the contour. + minmaxdist is the great-circle distance away from the candidate + to search for the minima/maxima. If delta is positive (negative), + the pivot is a local minimum (maximum). + (default: SEA_LEVEL_PRESSURE,200.0,5.5,0;thickness,-6.0,6.5,1.0) + +--no_closed_contour_cmd NO_CLOSED_CONTOUR_CMD + As --closed_contour_cmd, + except it eliminates candidates if a closed contour is present. + +--candidate_threshold_cmd CANDIDATE_THRESHOLD_CMD + Threshold commands for candidates. + Eliminate candidates that do not satisfy a threshold criteria + (there must exist a point within a given distance of the candidate + that satisfies a given equality or inequality). + Search is performed by breadth-first search over the grid. + Threshold commands are separated by a semicolon (";;..."). + Each threshold command takes the form "var,op,value,dist". + These arguments are as follows. + var is the name of the variable used for the thresholding. + op is the operator that must be satisfied for threshold (options include >,>=,<,<=,=,!=). + value is the value on the right-hand-side of the comparison. + dist is the great-circle distance away from the candidate + to search for a point that satisfies the threshold. + +--output_cmd OUTPUT_CMD + Include additional columns in the candidates output file. + Each output command takes the form "var,op,dist". + These arguments are as follows. + var is the name of the variable used for output. + op is the operator that is applied over all points + within the specified distance of the candidate (options include max, min, avg, maxdist, mindist). + dist is the great-circle distance away from the candidate wherein the operator is applied. + (default: SEA_LEVEL_PRESSURE,min,0;surface_wind_speed,max,2;GEOPOTENTIAL_AT_SURFACE,min,0) + +--in_fmt IN_FMT + A comma-separated list of names of the auxiliary columns within the candidates output file. + (namely, the list must not include the time columns). + (default: i,j,lat,lon,SEA_LEVEL_PRESSURE,surface_wind_speed,GEOPOTENTIAL_AT_SURFACE) + +--range RANGE + The maximum distance between candidates along a path (in great-circle degrees). + (default: 8.0) + +--min_time MIN_TIME + The minimum length of a path either in terms of number of discrete times or as a duration, e.g. "24h". + Note that the duration of a path is computed as the difference between the final time and initial time, + so a "24h" duration correspond to 5 time steps in 6-hourly data (i.e. 0h,6,12,18,24UTC). + (default: 10) + +--min_endpoint_distance MIN_ENDPOINT_DIST + The minimum great-circle distance between the first candidate on a path and the last candidate (in degrees). + (default: 0.0) + +--min_path_distance MIN_PATH_DISTANCE + The minimum accumulated great-circle distance between nodes in a path (in degrees). + (default: 0.0) + +--min_path_length MIN_PATH_LENGTH + Minimum path length (default: 0.0) + +--max_gap MAX_GAP + The number of allowed missing points between spatially proximal candidate nodes + while still considering them part of the same path. + (default: 3) + +--track_threshold_cmd TRACK_THRESHOLD_CMD + Filter paths based on the number of times where a particular threshold is satisfied. + Threshold commands are separated by a semicolon (";;..."). + Each threshold command takes the form "col,op,value,count". + These arguments are as follows. + col is the name of the column to use for thresholding, as specified in --in_fmt. + op is the operator that must be satisfied for threshold (options include >,>=,<,<=,=,!=,|>=,|<=). + value is the value on the right-hand-side of the comparison. + count is either the minimum number of time slices where the threshold must be satisfied + or the instruction "all", "first", or "last". + Here "all" is used to indicate the threshold must be satisfied at all points along the path, + "first" is used to indicate the threshold must be satisfied only at the first point along the path, and + "last" is used to indicate the threshold must be satisfied only at the last point along the path. + (default: surface_wind_speed,>=,10.0,10;lat,<=,50.0,10;lat,>=,-50.0,10;GEOPOTENTIAL_AT_SURFACE,<=,15.0,10) + +--prioritize VARIABLE_NAME + Variable to use when prioritizing path + +--allow_repeated_times true/false + Allow repeated times (default: false) + +--cal_type CAL_TYPE + Calendar type (default: standard) + +--first_step FIRST_STEP + first time step to process (default: 0) + +--last_step LAST_STEP + last time step to process (default: -1) + +--n_threads NUMBER_OF_THREADS + Sets the thread pool size on each MPI rank. When the default value of -1 is used TECA will + coordinate the thread pools across ranks such each thread is bound to a unique physical core. + (default: -1) + +--help + displays documentation for application specific command line options + +--advanced_help + displays documentation for algorithm specific command line options + +--full_help + displays both basic and advanced documentation together + +--verbose VERBOSE + Enable verbose output (default: 0) + +Examples +~~~~~~~~ + +ERA5 data +^^^^^^^^^ + +.. code-block:: bash + + #!/bin/bash + #SBATCH -C cpu + #SBATCH -q debug + #SBATCH -A m1517 + #SBATCH -N 1 + #SBATCH -n 64 + #SBATCH -t 00:30:00 + #SBATCH -J teca_tempest_tc_detect.out + #SBATCH --output=%x-%j.out + + # load the TECA module + module use /global/common/software/m1517/teca/perlmutter_cpu/develop/modulefiles + module load teca + + srun -n 64 teca_tempest_tc_detect \ + --input_file era5_2017.mcf \ + --sea_level_pressure MSL \ + --surface_wind_u VAR_10U \ + --surface_wind_v VAR_10V \ + --geopotential Z \ + --geopotential_at_surface ZS \ + --x_axis_variable longitude \ + --y_axis_variable latitude \ + --max_lat 90 \ + --min_lat -90 \ + --max_lon 359.75 \ + --min_lon 0 \ + --search_by_min MSL + + .. _teca_tc_detect: teca_tc_detect @@ -1595,7 +1903,7 @@ the following fields. Outputs ~~~~~~~ -1. Cyclone andidate table +1. Cyclone candidate table 2. Cyclone track table Command Line Arguments diff --git a/doc/rtd/bibliography.bib b/doc/rtd/bibliography.bib index c929cdb8d..7746c6c3d 100644 --- a/doc/rtd/bibliography.bib +++ b/doc/rtd/bibliography.bib @@ -81,3 +81,15 @@ @article{tc_props url = {https://doi.org/10.1175/MWR3420.1}, eprint = {https://journals.ametsoc.org/mwr/article-pdf/135/7/2568/4230462/mwr3420\_1.pdf}, } + +@article{tempestextremes, + author = {Ullrich, P. A. and Zarzycki, C. M. and McClenny, E. E. and Pinheiro, M. C. and Stansfield, A. M. and Reed, K. A.}, + title = {TempestExtremes v2.1: a community framework for feature detection, tracking, and analysis in large datasets}, + journaL = {Geoscientific Model Development}, + volume = {14}, + year = {2021}, + number = {8}, + pages = {5023--5048}, + url = {https://gmd.copernicus.org/articles/14/5023/2021/}, + doi = {10.5194/gmd-14-5023-2021} +} diff --git a/teca_config.h.in b/teca_config.h.in index 0fee0866b..b4146943f 100644 --- a/teca_config.h.in +++ b/teca_config.h.in @@ -19,6 +19,7 @@ #cmakedefine TECA_HAS_OPENSSL #cmakedefine TECA_HAS_SHAPELIB #cmakedefine TECA_HAS_CUDA +#cmakedefine TECA_HAS_TELITE #cmakedefine TECA_VERSION_DESCR "@TECA_VERSION_DESCR@" #cmakedefine TECA_PYTHON_VERSION @TECA_PYTHON_VERSION@ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index dc517624a..a676c7bce 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -36,6 +36,9 @@ if (TECA_HAS_DATA) configure_file(${TECA_DATA_ROOT}/ERA5_lapse_rate_group_test.mcf.in ${CMAKE_CURRENT_BINARY_DIR}/ERA5_lapse_rate_group_test.mcf) + + configure_file(${TECA_DATA_ROOT}/ERA5/ERA5_TC_test.mcf.in + ${CMAKE_CURRENT_BINARY_DIR}/ERA5_TC_test.mcf) endif() teca_add_test(test_pipeline_time_average @@ -1064,3 +1067,15 @@ teca_add_test(test_regional_surface_area 512 "${TECA_DATA_ROOT}/vc965bq8111" 17814000.0 3e3 FEATURES ${TECA_HAS_NETCDF} ${TECA_HAS_SHAPELIB} REQ_TECA_DATA) + +teca_add_test(test_mesh_join + EXEC_NAME test_mesh_join + SOURCES test_mesh_join.cpp + LIBS teca_core teca_data teca_io teca_alg ${teca_test_link} + COMMAND test_mesh_join + -i "${TECA_DATA_ROOT}/HighResMIP/TC_test/PSL/PSL.*\.nc" + "${TECA_DATA_ROOT}/HighResMIP/TC_test/U850/U850.*\.nc" + "${TECA_DATA_ROOT}/HighResMIP/TC_test/V850/V850.*\.nc" + -o "test_mesh_join.nc" -s 0,-1 -x lon -y lat -t time -c 4 -n 1 -p PSL U850 V850 + FEATURES ${TECA_HAS_NETCDF} + REQ_TECA_DATA) diff --git a/test/apps/CMakeLists.txt b/test/apps/CMakeLists.txt index 7dc3234b5..44abe2f2c 100644 --- a/test/apps/CMakeLists.txt +++ b/test/apps/CMakeLists.txt @@ -134,6 +134,12 @@ teca_add_app_test(test_tc_trajectory_app ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT} REQ_TECA_DATA) +teca_add_app_test(test_tempest_tc_detect_app_mcf teca_tempest_tc_detect + COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_tempest_tc_detect_app_mcf.sh + ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT} + FEATURES ${TECA_HAS_TELITE} ${TECA_HAS_NETCDF} + REQ_TECA_DATA) + teca_add_app_test(test_tc_wind_radii_app_serial teca_tc_wiind_radii COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_tc_wind_radii_app.sh @@ -792,7 +798,8 @@ set(app_names teca_tc_wind_radii teca_tc_wind_radii_stats teca_temporal_reduction - teca_cpp_temporal_reduction) + teca_cpp_temporal_reduction + teca_tempest_tc_detect) foreach (app_name ${app_names}) diff --git a/test/apps/test_tempest_tc_detect_app_mcf.sh b/test/apps/test_tempest_tc_detect_app_mcf.sh new file mode 100755 index 000000000..7cf5fd22d --- /dev/null +++ b/test/apps/test_tempest_tc_detect_app_mcf.sh @@ -0,0 +1,45 @@ +#!/bin/bash + +if [[ $# < 2 ]] +then + echo "usage: test_tempest_tc_detect_app_mcf.sh [app prefix] [data root] " \ + "[mpi exec] [test cores] $#" + exit -1 +fi + +app_prefix=${1} +data_root=${2} + +if [[ $# -eq 4 ]] +then + mpi_exec=${3} + test_cores=${4} + launcher="${mpi_exec} -n ${test_cores}" +fi + +set -x + +# run the app +${launcher} ${app_prefix}/teca_tempest_tc_detect \ + --input_file "${app_prefix}/../test/ERA5_TC_test.mcf" \ + --sea_level_pressure MSL \ + --geopotential Z \ + --surface_wind_u VAR_10U \ + --surface_wind_v VAR_10V \ + --geopotential_at_surface ZS \ + --x_axis_variable longitude \ + --y_axis_variable latitude \ + --max_lat 90 \ + --min_lat -90 \ + --max_lon 359.5 \ + --min_lon 0 \ + --last_step 1 \ + --candidate_file test_te_candidates_app_output.csv + +# run the diff +${app_prefix}/teca_table_diff \ + "${data_root}/test_te_candidates_app_ref.csv" \ + test_te_candidates_app_output.csv + +# clean up +rm test_te_candidates_app_output.csv diff --git a/test/test_mesh_join.cpp b/test/test_mesh_join.cpp new file mode 100644 index 000000000..6220060cc --- /dev/null +++ b/test/test_mesh_join.cpp @@ -0,0 +1,190 @@ +#include "teca_config.h" +#include "teca_cf_reader.h" +#include "teca_normalize_coordinates.h" +#include "teca_mesh_join.h" +#include "teca_cf_writer.h" +#include "teca_index_executive.h" +#include "teca_mpi_manager.h" +#include "teca_system_interface.h" +#include "teca_mpi.h" + +#include +#include +#include +using namespace std; + +int parse_command_line( + int argc, + char **argv, + int rank, + const p_teca_cf_reader &cf_reader_1, + const p_teca_cf_reader &cf_reader_2, + const p_teca_cf_reader &cf_reader_3, + const p_teca_cf_writer &cf_writer, + const p_teca_index_executive exec); + + +int main(int argc, char **argv) +{ + teca_mpi_manager mpi_man(argc, argv); + int rank = mpi_man.get_comm_rank(); + + teca_system_interface::set_stack_trace_on_error(); + teca_system_interface::set_stack_trace_on_mpi_error(); + + // create the pipeline objects + p_teca_cf_reader cf_reader_1 = teca_cf_reader::New(); + p_teca_cf_reader cf_reader_2 = teca_cf_reader::New(); + p_teca_cf_reader cf_reader_3 = teca_cf_reader::New(); + p_teca_cf_writer cf_writer = teca_cf_writer::New(); + p_teca_index_executive exec = teca_index_executive::New(); + + // initialize them from command line options + if (parse_command_line(argc, argv, rank, cf_reader_1, + cf_reader_2, cf_reader_3, cf_writer, exec)) + return -1; + + p_teca_mesh_join join = teca_mesh_join::New(); + join->set_number_of_input_connections(3); + join->set_verbose(1); + join->set_input_connection(0, cf_reader_1->get_output_port()); + join->set_input_connection(1, cf_reader_2->get_output_port()); + join->set_input_connection(2, cf_reader_3->get_output_port()); + + p_teca_normalize_coordinates coords = teca_normalize_coordinates::New(); + coords->set_input_connection(join->get_output_port()); + + cf_writer->set_input_connection(coords->get_output_port()); + cf_writer->set_executive(exec); + + // run the pipeline + cf_writer->update(); + + return 0; +} + + +// -------------------------------------------------------------------------- +int parse_command_line(int argc, char **argv, int rank, + const p_teca_cf_reader &cf_reader_1, + const p_teca_cf_reader &cf_reader_2, + const p_teca_cf_reader &cf_reader_3, + const p_teca_cf_writer &cf_writer, + const p_teca_index_executive exec) +{ + if (argc < 3) + { + if (rank == 0) + { + cerr << endl << "Usage error:" << endl + << "test_mesh_join [-i input1 input2 input3] [-o output] [-s first step,last step] " + << "[-x x axis variable] [-y y axis variable] [-z z axis variable] " + << "[-t t axis variable] [-c steps per file] [-n num threads] " + << "[-p var0 var1 ...]" + << endl << endl; + } + return -1; + } + + int n_files = 3; + vector input; + string output; + string x_ax = "lon"; + string y_ax = "lat"; + string z_ax = ""; + string t_ax = ""; + int n_threads = -1; + long first_step = 0; + long last_step = -1; + long steps_per_file = 1; + vector arrays; + + int j = 0; + for (int i = 1; i < argc; ++i) + { + if (!strcmp("-i", argv[i])) + { + for (int k = 0; k < n_files; ++k) + input.push_back(argv[++i]); + ++j; + } + else if (!strcmp("-o", argv[i])) + { + output = argv[++i]; + ++j; + } + else if (!strcmp("-x", argv[i])) + { + x_ax = argv[++i]; + ++j; + } + else if (!strcmp("-y", argv[i])) + { + y_ax = argv[++i]; + ++j; + } + else if (!strcmp("-z", argv[i])) + { + z_ax = argv[++i]; + ++j; + } + else if (!strcmp("-t", argv[i])) + { + t_ax = argv[++i]; + ++j; + } + else if (!strcmp("-s", argv[i])) + { + sscanf(argv[++i], "%li,%li", + &first_step, &last_step); + ++j; + } + else if (!strcmp("-n", argv[i])) + { + n_threads = atoi(argv[++i]); + ++j; + } + else if (!strcmp("-c", argv[i])) + { + steps_per_file = atoi(argv[++i]); + ++j; + } + if (!strcmp("-p", argv[i])) + { + for (int k = i; k < argc-1; ++k) + arrays.push_back(argv[++i]); + ++j; + } + } + + // pass the command line options + cf_reader_1->set_x_axis_variable(x_ax); + cf_reader_1->set_y_axis_variable(y_ax); + cf_reader_1->set_z_axis_variable(z_ax); + cf_reader_1->set_t_axis_variable(t_ax); + cf_reader_1->set_files_regex(input[0]); + + cf_reader_2->set_x_axis_variable(x_ax); + cf_reader_2->set_y_axis_variable(y_ax); + cf_reader_2->set_z_axis_variable(z_ax); + cf_reader_2->set_t_axis_variable(t_ax); + cf_reader_2->set_files_regex(input[1]); + + cf_reader_3->set_x_axis_variable(x_ax); + cf_reader_3->set_y_axis_variable(y_ax); + cf_reader_3->set_z_axis_variable(z_ax); + cf_reader_3->set_t_axis_variable(t_ax); + cf_reader_3->set_files_regex(input[2]); + + cf_writer->set_file_name(output); + cf_writer->set_thread_pool_size(n_threads); + cf_writer->set_first_step(first_step); + cf_writer->set_last_step(last_step); + cf_writer->set_layout(teca_cf_writer::number_of_steps); + cf_writer->set_steps_per_file(steps_per_file); + cf_writer->set_point_arrays(arrays); + + exec->set_arrays(arrays); + + return 0; +} diff --git a/test/travis_ci/ctest_linux.cmake b/test/travis_ci/ctest_linux.cmake index 63ea5bf67..fa848a9af 100644 --- a/test/travis_ci/ctest_linux.cmake +++ b/test/travis_ci/ctest_linux.cmake @@ -29,7 +29,9 @@ REQUIRE_NETCDF_MPI=$ENV{REQUIRE_NETCDF_MPI} REQUIRE_UDUNITS=TRUE REQUIRE_MPI=TRUE REQUIRE_PYTHON=TRUE -REQUIRE_TECA_DATA=TRUE") +REQUIRE_TECA_DATA=TRUE +REQUIRE_TELITE=TRUE +TELite_DIR=/usr/local/lib") file(WRITE "${CTEST_BINARY_DIRECTORY}/CMakeCache.txt" ${INITIAL_CACHE}) ctest_start("${DASHBOARD_TRACK}") ctest_read_custom_files("${CTEST_BINARY_DIRECTORY}") diff --git a/test/travis_ci/ctest_linux.sh b/test/travis_ci/ctest_linux.sh index 5aecb1344..d8babca5d 100755 --- a/test/travis_ci/ctest_linux.sh +++ b/test/travis_ci/ctest_linux.sh @@ -27,7 +27,7 @@ export LD_LIBRARY_PATH=${DASHROOT}/build/lib export MPLBACKEND=Agg mkdir build -ctest -S ${DASHROOT}/test/travis_ci/ctest_linux.cmake --output-on-failure --timeout 400 & +ctest -S ${DASHROOT}/test/travis_ci/ctest_linux.cmake -V --timeout 400 & ctest_pid=$! # this loop prevents travis from killing the job diff --git a/test/travis_ci/install_fedora_31.sh b/test/travis_ci/install_fedora_31.sh index 1e97f266b..a41074d23 100755 --- a/test/travis_ci/install_fedora_31.sh +++ b/test/travis_ci/install_fedora_31.sh @@ -37,3 +37,12 @@ pip3 install "numpy<2.0" mpi4py matplotlib torch # install data files. svn co svn://svn.code.sf.net/p/teca/TECA_data@${TECA_DATA_REVISION} TECA_data + +# install TELite library +git clone https://github.com/LBL-EESA/TELite.git +cd TELite +git checkout ${TECA_TELITE_REVISION} +mkdir build +cd build +cmake .. +make -j2 && make -j2 install diff --git a/test/travis_ci/install_fedora_32.sh b/test/travis_ci/install_fedora_32.sh index 1e97f266b..f754e70ee 100755 --- a/test/travis_ci/install_fedora_32.sh +++ b/test/travis_ci/install_fedora_32.sh @@ -28,6 +28,7 @@ echo ${DOCKER_IMAGE} echo ${IMAGE_VERSION} echo ${TECA_PYTHON_VERSION} echo ${TECA_DATA_REVISION} +echo ${TECA_TELITE_REVISION} python3 -mvenv `pwd`/../tci set +x @@ -37,3 +38,12 @@ pip3 install "numpy<2.0" mpi4py matplotlib torch # install data files. svn co svn://svn.code.sf.net/p/teca/TECA_data@${TECA_DATA_REVISION} TECA_data + +# install TELite library +git clone https://github.com/LBL-EESA/TELite.git +cd TELite +git checkout ${TECA_TELITE_REVISION} +mkdir build +cd build +cmake .. +make -j2 && make -j2 install diff --git a/test/travis_ci/install_fedora_33.sh b/test/travis_ci/install_fedora_33.sh index 1e97f266b..f754e70ee 100755 --- a/test/travis_ci/install_fedora_33.sh +++ b/test/travis_ci/install_fedora_33.sh @@ -28,6 +28,7 @@ echo ${DOCKER_IMAGE} echo ${IMAGE_VERSION} echo ${TECA_PYTHON_VERSION} echo ${TECA_DATA_REVISION} +echo ${TECA_TELITE_REVISION} python3 -mvenv `pwd`/../tci set +x @@ -37,3 +38,12 @@ pip3 install "numpy<2.0" mpi4py matplotlib torch # install data files. svn co svn://svn.code.sf.net/p/teca/TECA_data@${TECA_DATA_REVISION} TECA_data + +# install TELite library +git clone https://github.com/LBL-EESA/TELite.git +cd TELite +git checkout ${TECA_TELITE_REVISION} +mkdir build +cd build +cmake .. +make -j2 && make -j2 install diff --git a/test/travis_ci/install_fedora_37.sh b/test/travis_ci/install_fedora_37.sh index 5800b2ea6..fc4cb457e 100755 --- a/test/travis_ci/install_fedora_37.sh +++ b/test/travis_ci/install_fedora_37.sh @@ -40,3 +40,12 @@ pip3 install "numpy<2.0" mpi4py matplotlib torch # install data files. svn co svn://svn.code.sf.net/p/teca/TECA_data@${TECA_DATA_REVISION} TECA_data + +# install TELite library +git clone https://github.com/LBL-EESA/TELite.git +cd TELite +git checkout ${TECA_TELITE_REVISION} +mkdir build +cd build +cmake .. +make -j2 && make -j2 install diff --git a/test/travis_ci/install_ubuntu_20_04.sh b/test/travis_ci/install_ubuntu_20_04.sh index 9dfe548c4..eb52e518b 100755 --- a/test/travis_ci/install_ubuntu_20_04.sh +++ b/test/travis_ci/install_ubuntu_20_04.sh @@ -36,6 +36,7 @@ echo ${DOCKER_IMAGE} echo ${IMAGE_VERSION} echo ${TECA_PYTHON_VERSION} echo ${TECA_DATA_REVISION} +echo ${TECA_TELITE_REVISION} python3 -mvenv `pwd`/../tci set +x @@ -45,3 +46,12 @@ pip3 install "numpy<2.0" mpi4py matplotlib torch # install data files. svn co svn://svn.code.sf.net/p/teca/TECA_data@${TECA_DATA_REVISION} TECA_data + +# install TELite library +git clone https://github.com/LBL-EESA/TELite.git +cd TELite +git checkout ${TECA_TELITE_REVISION} +mkdir build +cd build +cmake .. +make -j2 && make -j2 install diff --git a/test/travis_ci/install_ubuntu_22_04.sh b/test/travis_ci/install_ubuntu_22_04.sh index 9dfe548c4..c39b6dcc7 100755 --- a/test/travis_ci/install_ubuntu_22_04.sh +++ b/test/travis_ci/install_ubuntu_22_04.sh @@ -45,3 +45,12 @@ pip3 install "numpy<2.0" mpi4py matplotlib torch # install data files. svn co svn://svn.code.sf.net/p/teca/TECA_data@${TECA_DATA_REVISION} TECA_data + +# install TELite library +git clone https://github.com/LBL-EESA/TELite.git +cd TELite +git checkout ${TECA_TELITE_REVISION} +mkdir build +cd build +cmake .. +make -j2 && make -j2 install