diff --git a/.gitignore b/.gitignore
index bf6a22fc6..957306f03 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,6 @@
*.patch
_build
*.pt
+*.vscode*
+.DS_Store
+generated_rtd*
diff --git a/.readthedocs.yml b/.readthedocs.yaml
similarity index 86%
rename from .readthedocs.yml
rename to .readthedocs.yaml
index e2c25ef5e..16f057784 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yaml
@@ -4,10 +4,13 @@
version: 2
+build:
+ image: latest
+
sphinx:
configuration: doc/rtd/conf.py
python:
- version: 3.7
+ version: 3.8
install:
- requirements: doc/rtd/requirements.txt
diff --git a/.travis.yml b/.travis.yml
index 909717b1d..72b7b9552 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -17,7 +17,7 @@ env:
- BUILD_TYPE=Debug
- TECA_DIR=/travis_teca_dir
- TECA_PYTHON_VERSION=3
- - TECA_DATA_REVISION=101
+ - TECA_DATA_REVISION=117
jobs:
- DOCKER_IMAGE=ubuntu IMAGE_VERSION=20.04 IMAGE_NAME=ubuntu_20_04 REQUIRE_NETCDF_MPI=TRUE
- DOCKER_IMAGE=ubuntu IMAGE_VERSION=20.04 IMAGE_NAME=ubuntu_20_04 REQUIRE_NETCDF_MPI=FALSE
diff --git a/CMake/teca_app.cmake b/CMake/teca_app.cmake
index fbaa0e923..82bb06234 100644
--- a/CMake/teca_app.cmake
+++ b/CMake/teca_app.cmake
@@ -31,6 +31,7 @@ function (teca_add_app app_name)
teca_system teca_core teca_data teca_io teca_alg
${APP_LIBS})
endif()
+ set_target_properties(${app_name} PROPERTIES APP_TYPE C++)
install(TARGETS ${app_name} RUNTIME DESTINATION ${BIN_PREFIX})
else()
message(STATUS "command line application ${app_name} -- disabled")
diff --git a/CMake/teca_python.cmake b/CMake/teca_python.cmake
index 6ee3f7397..b94dcb80c 100644
--- a/CMake/teca_python.cmake
+++ b/CMake/teca_python.cmake
@@ -60,6 +60,8 @@ function (teca_add_python_app app_name)
if (NOT APP_SOURCES)
set(APP_SOURCES "${app_name}.in")
endif()
+ add_custom_target(${app_name})
+ set_target_properties(${app_name} PROPERTIES APP_TYPE Python)
teca_py_install_apps(${APP_SOURCES})
else()
message(STATUS "command line application ${app_name} -- disabled")
diff --git a/CMake/teca_test.cmake b/CMake/teca_test.cmake
index 375c3524f..346eaf319 100644
--- a/CMake/teca_test.cmake
+++ b/CMake/teca_test.cmake
@@ -51,3 +51,9 @@ function (teca_add_test T_NAME)
endif()
endif()
endfunction()
+
+function (teca_add_app_test T_NAME T_TARGET)
+ if (TARGET ${T_TARGET})
+ teca_add_test(${T_NAME} ${ARGV})
+ endif()
+endfunction()
diff --git a/CMakeLists.txt b/CMakeLists.txt
index c1ecb1fbb..4eb5ff0ff 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -387,26 +387,50 @@ if (BUILD_TESTING)
# figure out how many cores we can use for parallel tests
set(TECA_TEST_CORES 0 CACHE STRING
"Max number of cores for use in parallel tests")
+
+ # by default assume 2 hyperthreads per core, if this is not
+ # the case override here
+ set(HYPERTHREADS_PER_CORE 2 CACHE STRING
+ "The number of hyperthreads per core.")
+
+ # use CMake to get the number of logical cores. includes hyperthreads
+ # in the count.
if (TECA_TEST_CORES LESS 1)
ProcessorCount(LOGICAL_CORES)
if (LOGICAL_CORES EQUAL 0)
- set(LOGICAL_CORES 4)
+ message(FATAL_ERROR "Failed to detect the number of cores. "
+ "Set TECA_TEST_CORES")
endif()
else()
- math(EXPR LOGICAL_CORES "${TECA_TEST_CORES}*2")
+ math(EXPR LOGICAL_CORES "${TECA_TEST_CORES}*${HYPERTHREADS_PER_CORE}")
endif()
- math(EXPR PHYSICAL_CORES "${LOGICAL_CORES}/2")
- if (PHYSICAL_CORES LESS 3)
- set(TEST_CORES 2)
- set(HALF_TEST_CORES 2)
- set(TWICE_TEST_CORES 4)
- else()
- set(TEST_CORES ${PHYSICAL_CORES})
- math(EXPR HALF_TEST_CORES "${TEST_CORES}/2")
- set(TWICE_TEST_CORES ${LOGICAL_CORES})
+
+ # adjust count for hyperthreads.
+ math(EXPR PHYSICAL_CORES "${LOGICAL_CORES}/${HYPERTHREADS_PER_CORE}")
+ if (PHYSICAL_CORES LESS 1)
+ message(FATAL_ERROR "Invalid CPU configuration. "
+ "LOGICAL_CORES=${LOGICAL_CORES} HYPERTHREADS_PER_CORE="
+ "${HYPERTHREADS_PER_CORE}")
endif()
+
+ # set the number of cores to use for pure MPI or purely threaded tests
+ set(TEST_CORES ${PHYSICAL_CORES})
message(STATUS "regression testing -- enabled (${TEST_CORES} cores).")
+ # set the number of cores to use for MPI + threads tests. if there are too
+ # few physical cores then disable hybrid parallel tests
+ math(EXPR HALF_TEST_CORES "${TEST_CORES}/2")
+ if (HALF_TEST_CORES LESS 2)
+ message(STATUS "Hybrid parallel tests -- disabled.")
+ set(TEST_MPI_THREADS OFF)
+ else()
+ message(STATUS "Hybrid parallel tests -- enabled.")
+ set(TEST_MPI_THREADS ON)
+ endif()
+
+ # set the number of cores for oversubscription/streaming tests
+ math(EXPR TWICE_TEST_CORES "${TEST_CORES}*2")
+
add_subdirectory(test)
else()
message(STATUS "regression testing -- disbaled")
diff --git a/README.md b/README.md
index 8e0485626..a10c5ad87 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,8 @@
-
+
+
+
+ |
+
@@ -8,7 +12,7 @@
TECA is a collection of climate analysis algorithms geared toward extreme event detection and tracking implemented in a scalable parallel framework. The code has been successfully deployed and run at massive scales on current DOE supercomputers. TECA's core is written in modern C++ and exploits MPI + X parallelism where X is one of threads, OpenMP, or GPUs. The framework supports a number of parallel design patterns including distributed data parallelism and map-reduce. While modern C++ delivers the highest performance, Python bindings make the code approachable and easy to use.
### Documentation
-The [TECA User's Guide](https://teca.readthedocs.io/en/latest/) is the authorotative source for documentation on topics such as [installing TECA](https://teca.readthedocs.io/en/latest/installation.html), running TECA's [command line applications](https://teca.readthedocs.io/en/latest/applications.html), and [Python development](https://teca.readthedocs.io/en/latest/python.html).
+The [TECA User's Guide](https://teca.readthedocs.io/en/latest/) is the authorotative source for documentation on topics such as [installing TECA](https://teca.readthedocs.io/en/latest/installation.html), running TECA's [command line applications](https://teca.readthedocs.io/en/latest/applications.html), and [Python development](https://teca.readthedocs.io/en/latest/python.html). The TECA source code is documented on our [Doxygen site](https://teca.readthedocs.io/en/latest/doxygen/index.html).
### Tutorials
The [TECA tutorials](https://sourceforge.net/p/teca/TECA_tutorials) subversion repository contains slides from previous tutorials.
diff --git a/alg/CMakeLists.txt b/alg/CMakeLists.txt
index 48038978d..d8f7e1734 100644
--- a/alg/CMakeLists.txt
+++ b/alg/CMakeLists.txt
@@ -21,16 +21,20 @@ set(teca_alg_cxx_srcs
teca_component_statistics.cxx
teca_derived_quantity.cxx
teca_descriptive_statistics.cxx
+ teca_elevation_mask.cxx
teca_evaluate_expression.cxx
teca_face_to_cell_centering.cxx
teca_geography.cxx
+ teca_indexed_dataset_cache.cxx
teca_integrated_vapor_transport.cxx
+ teca_integrated_water_vapor.cxx
teca_l2_norm.cxx
teca_latitude_damper.cxx
teca_laplacian.cxx
teca_mask.cxx
teca_normalize_coordinates.cxx
teca_parser.cxx
+ teca_rename_variables.cxx
teca_table_calendar.cxx
teca_table_reduce.cxx
teca_table_region_mask.cxx
@@ -41,7 +45,8 @@ set(teca_alg_cxx_srcs
teca_tc_classify.cxx
teca_tc_wind_radii.cxx
teca_tc_trajectory.cxx
- teca_temporal_average.cxx
+ teca_simple_moving_average.cxx
+ teca_unpack_data.cxx
teca_valid_value_mask.cxx
teca_variant_array_operand.cxx
teca_vertical_coordinate_transform.cxx
diff --git a/alg/teca_2d_component_area.cxx b/alg/teca_2d_component_area.cxx
index 285d40f52..9bdcbac60 100644
--- a/alg/teca_2d_component_area.cxx
+++ b/alg/teca_2d_component_area.cxx
@@ -132,11 +132,13 @@ void teca_2d_component_area::get_properties_description(
"name of the varibale containing region labels")
TECA_POPTS_GET(int, prefix, contiguous_component_ids,
"when the region label ids start at 0 and are consecutive "
- "this flag enables use of an optimization (0)")
+ "this flag enables use of an optimization")
TECA_POPTS_GET(long, prefix, background_id,
- "the label id that corresponds to the background (-1)")
+ "the label id that corresponds to the background")
;
+ this->teca_algorithm::get_properties_description(prefix, opts);
+
global_opts.add(opts);
}
@@ -144,6 +146,8 @@ void teca_2d_component_area::get_properties_description(
void teca_2d_component_area::set_properties(const std::string &prefix,
variables_map &opts)
{
+ this->teca_algorithm::set_properties(prefix, opts);
+
TECA_POPTS_SET(opts, std::string, prefix, component_variable)
TECA_POPTS_SET(opts, int, prefix, contiguous_component_ids)
TECA_POPTS_SET(opts, long, prefix, background_id)
diff --git a/alg/teca_2d_component_area.h b/alg/teca_2d_component_area.h
index 5a9d57e6a..54f4a2e42 100644
--- a/alg/teca_2d_component_area.h
+++ b/alg/teca_2d_component_area.h
@@ -10,40 +10,40 @@
TECA_SHARED_OBJECT_FORWARD_DECL(teca_2d_component_area)
-/// an algorithm that computes the area of labeled regions
+/// An algorithm that computes the areas of labeled regions
/**
-Given a set of labels on a Cartesian mesh, the algorithm computes the area of
-each region. Regions are identified by assigning a unique integer value to each
-mesh point that belongs in the region. The component_variable property names
-the variable containing the region labels.
-
-if the region labels start at 0 and are contiguous then an optimization can be
-used. Set contiguous_component_ids property to enable the optimization. Note that
-TECA's connected component labeler assigns the background (i.e. cells not inside
-the segmentation) the label 0. One can identify the background region and area
-via this label. When processing data generated outside of TECA it might be
-necessary to supply the background label. Use -2 if there is no background.
-
-the input dataset is passed through and the results of the calculations are
-stored in the output dataset metadata in the following keys:
-
- number_of_components - number of component ids for which area was
- computed. Note that this can include a background
- component i.e. for cells outside of the segmentation.
-
- component_ids - a vector containing the label of each component. This is
- always starts with 0, where the label 0 identifies cells
- out side of the segmentation, and ranges up to
- number_of_components - 1, where the labels from 1 up to
- number_of_components - 1 identify connected regions of
- cells inside the segmentation.
-
- component_area - a vector containing the area for the corresponding entry
- in the component_ids array.
-
- background_id - the label used for cells outside of the segmentation,
- i.e. the background. This can be used to skip processing
- of the background when desirable.
+ * Given a set of labels on a Cartesian mesh, the algorithm computes the area
+ * of each region. Regions are identified by assigning a unique integer value
+ * to each mesh point that belongs in the region. The component_variable
+ * property names the variable containing the region labels.
+ *
+ * if the region labels start at 0 and are contiguous then an optimization can
+ * be used. Set contiguous_component_ids property to enable the optimization.
+ * Note that TECA's connected component labeler assigns the background (i.e.
+ * cells not inside the segmentation) the label 0. One can identify the
+ * background region and area via this label. When processing data generated
+ * outside of TECA it might be necessary to supply the background label. Use -2
+ * if there is no background.
+ *
+ * the input dataset is passed through and the results of the calculations are
+ * stored in the output dataset metadata in the following keys:
+ *
+ * | name | description |
+ * | ---- | ----------- |
+ * | number_of_components | number of component ids for which area was |
+ * | | computed. Note that this can include a background |
+ * | | component i.e. for cells outside of the segmentation. |
+ * | component_ids | a vector containing the label of each component. This is |
+ * | | always starts with 0, where the label 0 identifies |
+ * | | cells out side of the segmentation, and ranges up |
+ * | | to number_of_components - 1, where the labels from |
+ * | | 1 up to number_of_components - 1 identify |
+ * | | connected regions of cells inside the segmentation. |
+ * | component_area | a vector containing the area for the corresponding |
+ * | | entry in the component_ids array. |
+ * | background_id | the label used for cells outside of the segmentation, |
+ * | | i.e. the background. This can be used to skip processing |
+ * | | of the background when desirable. |
*/
class teca_2d_component_area : public teca_algorithm
{
@@ -53,23 +53,40 @@ class teca_2d_component_area : public teca_algorithm
TECA_ALGORITHM_CLASS_NAME(teca_2d_component_area)
~teca_2d_component_area();
- // report/initialize to/from Boost program options objects.
+ /** @name program_options
+ * report/initialize to/from Boost program options objects.
+ */
+ ///@{
TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION()
TECA_SET_ALGORITHM_PROPERTIES()
+ ///@}
- // set the name of the input array
+ /** @name component_variable
+ * Sets the name of the array containing component labels to compute the
+ * area of.
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(std::string, component_variable)
+ ///@}
- // set this only if you know for certain that label ids are contiguous and
- // start at 0. this enables use of a faster implementation.
+ /** @name contiguous_component_ids
+ * set this only if you know for certain that label ids are contiguous and
+ * start at 0. this enables use of a faster implementation.
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(int, contiguous_component_ids)
-
- // set this to override the component label used for background. By default
- // this is set to -1 to indicate that the value should be obtained from the
- // metadata key `background_id`. Note that TECA's connected component
- // labeler uses the id 0 for the background and passes this in a metadata
- // key and as a result no action is required.
+ ///@}
+
+ /** @name background_id
+ * set this to override the component label used for background. By default
+ * this is set to -1 to indicate that the value should be obtained from the
+ * metadata key `background_id`. Note that TECA's connected component
+ * labeler uses the id 0 for the background and passes this in a metadata
+ * key and as a result no action is required.
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(long, background_id)
+ ///@}
protected:
teca_2d_component_area();
diff --git a/alg/teca_apply_binary_mask.cxx b/alg/teca_apply_binary_mask.cxx
index 4ec27d4df..0e2d9aa8b 100644
--- a/alg/teca_apply_binary_mask.cxx
+++ b/alg/teca_apply_binary_mask.cxx
@@ -1,26 +1,44 @@
#include "teca_apply_binary_mask.h"
-#include "teca_mesh.h"
+#include "teca_cartesian_mesh.h"
#include "teca_array_collection.h"
#include "teca_variant_array.h"
#include "teca_metadata.h"
-#include "teca_mesh.h"
+#include "teca_array_attributes.h"
+#include "teca_mpi_util.h"
#include
#include
-#include
#include
-using std::deque;
-using std::vector;
-using std::set;
+#include
+
+#if defined(TECA_HAS_BOOST)
+#include
+#endif
+
+//#define TECA_DEBUG
+
using std::cerr;
using std::endl;
-//#define TECA_DEBUG
+namespace internal
+{
+// output = mask*input
+template
+void apply_mask(var_t * __restrict__ output,
+ const mask_t * __restrict__ mask,
+ const var_t * __restrict__ input,
+ unsigned long n)
+{
+ for (size_t i = 0; i < n; ++i)
+ output[i] = mask[i]*input[i];
+}
+}
// --------------------------------------------------------------------------
-teca_apply_binary_mask::teca_apply_binary_mask() : mask_variable("")
+teca_apply_binary_mask::teca_apply_binary_mask() :
+ mask_variable(""), output_variable_prefix("masked_")
{
this->set_number_of_input_connections(1);
this->set_number_of_output_ports(1);
@@ -30,6 +48,125 @@ teca_apply_binary_mask::teca_apply_binary_mask() : mask_variable("")
teca_apply_binary_mask::~teca_apply_binary_mask()
{}
+#if defined(TECA_HAS_BOOST)
+// --------------------------------------------------------------------------
+void teca_apply_binary_mask::get_properties_description(
+ const std::string &prefix, options_description &global_opts)
+{
+ options_description opts("Options for "
+ + (prefix.empty()?"teca_apply_binary_mask":prefix));
+
+ opts.add_options()
+ TECA_POPTS_MULTI_GET(std::vector, prefix, masked_variables,
+ "A list of variables to apply the mask to.")
+ TECA_POPTS_GET(std::string, prefix, mask_variable,
+ "The name of the variable containing the mask values.")
+ TECA_POPTS_GET(std::string, prefix, output_variable_prefix,
+ "A string prepended to the output variable names. If empty the"
+ " input variables will be replaced by their masked results")
+ ;
+
+ this->teca_algorithm::get_properties_description(prefix, opts);
+
+ global_opts.add(opts);
+}
+
+// --------------------------------------------------------------------------
+void teca_apply_binary_mask::set_properties(
+ const std::string &prefix, variables_map &opts)
+{
+ this->teca_algorithm::set_properties(prefix, opts);
+
+ TECA_POPTS_SET(opts, std::vector, prefix, masked_variables)
+ TECA_POPTS_SET(opts, std::string, prefix, mask_variable)
+ TECA_POPTS_SET(opts, std::string, prefix, output_variable_prefix)
+}
+#endif
+
+// --------------------------------------------------------------------------
+std::string teca_apply_binary_mask::get_output_variable_name(std::string input_var)
+{
+ return this->output_variable_prefix + input_var;
+}
+
+// --------------------------------------------------------------------------
+void teca_apply_binary_mask::get_output_variable_names(
+ std::vector &names)
+{
+ int n_inputs = this->masked_variables.size();
+ for (int i = 0; i < n_inputs; ++i)
+ {
+ names.push_back(
+ this->get_output_variable_name(this->masked_variables[i]));
+ }
+}
+
+// --------------------------------------------------------------------------
+teca_metadata teca_apply_binary_mask::get_output_metadata(
+ unsigned int port,
+ const std::vector &input_md)
+{
+#ifdef TECA_DEBUG
+ cerr << teca_parallel_id()
+ << "teca_apply_binary_mask::get_output_metadata" << endl;
+#endif
+ (void)port;
+
+ // check that the input variables have been specified.
+ // this is likely a user error.
+ if (this->masked_variables.empty() &&
+ teca_mpi_util::mpi_rank_0(this->get_communicator()))
+ {
+ TECA_WARNING("Nothing to do, masked_variables have not"
+ " been specified.")
+ }
+
+ // add in the array we will generate
+ teca_metadata out_md(input_md[0]);
+
+ // get the attributes
+ teca_metadata attributes;
+ out_md.get("attributes", attributes);
+
+ // construct the list of output variable names
+ for (auto& input_var : masked_variables)
+ {
+ std::string output_var = this->get_output_variable_name(input_var);
+
+ // add the varible to the list of output variables
+ out_md.append("variables", output_var);
+
+ // insert attributes to enable this variable to be written by the CF writer
+ teca_metadata input_atts;
+ if (attributes.get(input_var, input_atts))
+ {
+ TECA_WARNING("Failed to get attributes for \"" << input_var
+ << "\". Writing the result will not be possible")
+ }
+ else
+ {
+ // copy the attributes from the input. this will capture the
+ // data type, size, units, etc.
+ teca_array_attributes output_atts(input_atts);
+
+ // update description and long name
+ output_atts.description = input_var +
+ " multiplied by " + this->mask_variable;
+
+ output_atts.long_name.clear();
+
+ // update the array attributes
+ attributes.set(output_var, (teca_metadata)output_atts);
+ }
+
+ }
+
+ // update the attributes
+ out_md.set("attributes", attributes);
+
+ return out_md;
+}
+
// --------------------------------------------------------------------------
std::vector teca_apply_binary_mask::get_upstream_request(
unsigned int port, const std::vector &input_md,
@@ -42,9 +179,9 @@ std::vector teca_apply_binary_mask::get_upstream_request(
(void) port;
(void) input_md;
- vector up_reqs;
+ std::vector up_reqs;
- // get the name of the array to request
+ // get the name of the mask array
if (this->mask_variable.empty())
{
TECA_ERROR("A mask variable was not specified")
@@ -55,11 +192,36 @@ std::vector teca_apply_binary_mask::get_upstream_request(
// add in what we need
teca_metadata req(request);
std::set arrays;
+
if (req.has("arrays"))
req.get("arrays", arrays);
+
arrays.insert(this->mask_variable);
- if (!this->mask_arrays.empty())
- arrays.insert(this->mask_arrays.begin(), this->mask_arrays.end());
+
+ // check that the input variables have been specified.
+ // this is likely a user error.
+ if (this->masked_variables.empty() &&
+ teca_mpi_util::mpi_rank_0(this->get_communicator()))
+ {
+ TECA_WARNING("Nothing to do, masked_variables have not"
+ " been specified.")
+ }
+
+ // request the arrays to mask
+ for (auto& input_var : masked_variables)
+ {
+ // request the needed variable
+ arrays.insert(input_var);
+
+ // intercept request for our output if the variable will have a new name
+ std::string out_var = this->get_output_variable_name(input_var);
+ if (out_var != input_var)
+ {
+ arrays.erase(out_var);
+ }
+ }
+
+ // update the list of arrays to request
req.set("arrays", arrays);
// send up
@@ -79,68 +241,80 @@ const_p_teca_dataset teca_apply_binary_mask::execute(
(void)request;
// get the input
- const_p_teca_mesh in_mesh =
- std::dynamic_pointer_cast(input_data[0]);
+ const_p_teca_mesh in_mesh
+ = std::dynamic_pointer_cast(input_data[0]);
if (!in_mesh)
{
- TECA_ERROR("empty input, or not a mesh")
+ TECA_ERROR("Failed to apply mask. Dataset is not a teca_mesh")
return nullptr;
}
- // create output and copy metadata, coordinates, etc
+ // create the output mesh, pass everything through
+ // masked arrays are added or replaced below
p_teca_mesh out_mesh =
- std::dynamic_pointer_cast(in_mesh->new_instance());
- out_mesh->copy(in_mesh);
+ std::static_pointer_cast
+ (std::const_pointer_cast(in_mesh)->new_shallow_copy());
- // get the mask array
+ // check that a masking variable has been provided
if (this->mask_variable.empty())
{
- TECA_ERROR("A mask variable was not specified")
+ TECA_ERROR("The mask_variable name was not specified")
return nullptr;
}
- p_teca_array_collection arrays = out_mesh->get_point_arrays();
-
- p_teca_variant_array mask_array = arrays->get(this->mask_variable);
+ // get the mask array
+ const_p_teca_variant_array mask_array
+ = in_mesh->get_point_arrays()->get(this->mask_variable);
if (!mask_array)
{
- TECA_ERROR("mask variable \"" << this->mask_variable
- << "\" is not in the input")
+ TECA_ERROR("The mask_variable \"" << this->mask_variable
+ << "\" was requested but is not present in the input data.")
return nullptr;
}
// apply the mask
- unsigned long nelem = mask_array->size();
+ NESTED_TEMPLATE_DISPATCH(const teca_variant_array_impl,
+ mask_array.get(), _MASK,
- NESTED_TEMPLATE_DISPATCH(teca_variant_array_impl,
- mask_array.get(), _1,
+ // loop over input variables
+ for (auto& input_var : masked_variables)
+ {
+ std::string output_var = this->get_output_variable_name(input_var);
- NT_1 *pmask = static_cast(mask_array.get())->get();
+ // get the input array
+ const_p_teca_variant_array input_array
+ = in_mesh->get_point_arrays()->get(input_var);
+ if (!input_array)
+ {
+ TECA_ERROR("The masked_variable \"" << input_var
+ << "\" was requested but is not present in the input data.")
+ return nullptr;
+ }
- unsigned int narrays = arrays->size();
- for (unsigned int i = 0; i < narrays; ++i)
- {
- // if the user provided a list, restrict masking to that
- // list. and if not, mask everything
- if (!this->mask_arrays.empty() &&
- !std::count(this->mask_arrays.begin(),
- this->mask_arrays.end(), arrays->get_name(i)))
- continue;
+ // allocate the output array
+ size_t n = input_array->size();
- p_teca_variant_array array = arrays->get(i);
+ p_teca_variant_array output_array = input_array->new_instance(n);
- NESTED_TEMPLATE_DISPATCH(teca_variant_array_impl,
- array.get(), _2,
+ //output_array->resize(n);
- NT_2 *parray = static_cast(array.get())->get();
+ // do the mask calculation
+ NESTED_TEMPLATE_DISPATCH(
+ teca_variant_array_impl,
+ output_array.get(), _VAR,
- for (unsigned long q = 0; q < nelem; ++q)
- {
- parray[q] *= static_cast(pmask[q]);
- }
+ internal::apply_mask(
+ dynamic_cast(output_array.get())->get(),
+ static_cast(mask_array.get())->get(),
+ static_cast(input_array.get())->get(),
+ n);
)
+
+ out_mesh->get_point_arrays()->set(
+ output_var, output_array);
}
- )
+ )
+
return out_mesh;
}
diff --git a/alg/teca_apply_binary_mask.h b/alg/teca_apply_binary_mask.h
index f30a1081e..ceeed4040 100644
--- a/alg/teca_apply_binary_mask.h
+++ b/alg/teca_apply_binary_mask.h
@@ -10,12 +10,26 @@
TECA_SHARED_OBJECT_FORWARD_DECL(teca_apply_binary_mask)
-/// an algorithm that applies a binary mask multiplicatively
+/// Applies a mask to a given list of variables
/**
-an algorithm that applies a binary mask multiplicatively to all
-arrays in the input dataset. where mask is 1 values are passed
-through, where mask is 0 values are removed.
-*/
+ * Given a mask variable, this routine applies the mask to a list of input
+ * variables.
+ *
+ * The mask variable can either be binary, or it can represent a probability
+ * ranging from 0 to 1. For mask variable `mask` and input variable `var`, this
+ * algorithm computes `mask * var` and sends the resulting array downstream; this
+ * masking operation is applied for all variables in the input list.
+ *
+ * A potential use-case for this algorithm is masking quantities like
+ * precipitation by the probability of atmospheric river presence; the average
+ * of this masked precipitation variable gives the average precipitation due to
+ * atmospheric rivers.
+ *
+ * The output variable names are given a prefix to distinguish them from the
+ * upstream versions. E.g., if the algorithm property `output_variable_prefix` is set
+ * to 'ar_', and the variable being masked is 'precip', then the output array
+ * name is 'ar_precip'.
+ */
class teca_apply_binary_mask : public teca_algorithm
{
public:
@@ -24,19 +38,48 @@ class teca_apply_binary_mask : public teca_algorithm
TECA_ALGORITHM_CLASS_NAME(teca_apply_binary_mask)
~teca_apply_binary_mask();
- // set the name of the output array
+ // report/initialize to/from Boost program options
+ // objects.
+ TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION()
+ TECA_SET_ALGORITHM_PROPERTIES()
+
+ /** @name mask_variable
+ * set the name of the variable containing the mask values
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(std::string, mask_variable)
+ ///@}
+
+ /** @name masked_variable
+ * A list of of variables to apply the mask to. If empty no arrays will be
+ * requested, and no variables will be masked
+ */
+ ///@{
+ TECA_ALGORITHM_VECTOR_PROPERTY(std::string, masked_variable)
+ ///@}
- // set the arrays to mask. if empty no arrays will be
- // requested, but all present will be masked
- TECA_ALGORITHM_VECTOR_PROPERTY(std::string, mask_array)
+ /** @name output_variable_prefix
+ * A prefix for the names of the variables that have been masked. If this
+ * is empty masked data replaces its input, otherwise input data is
+ * preserved and masked data is added.
+ */
+ ///@{
+ TECA_ALGORITHM_PROPERTY(std::string, output_variable_prefix)
+ ///@}
+
+ /** helper that constructs and returns the result variable names taking
+ * into account he list of masked_variables and the output_variable_prefix.
+ * use this to know what variables will be produced.
+ */
+ void get_output_variable_names(std::vector &names);
protected:
teca_apply_binary_mask();
private:
- //teca_metadata get_output_metadata(unsigned int port,
- // const std::vector &input_md) override;
+ 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,
@@ -46,9 +89,14 @@ class teca_apply_binary_mask : public teca_algorithm
const std::vector &input_data,
const teca_metadata &request) override;
+ // helper that given and input variable name constructs the result variable
+ // name taking into account the output_variable_prefix
+ std::string get_output_variable_name(std::string input_var);
+
private:
std::string mask_variable;
- std::vector mask_arrays;
+ std::vector masked_variables;
+ std::string output_variable_prefix;
};
#endif
diff --git a/alg/teca_ar_detect.cxx b/alg/teca_ar_detect.cxx
deleted file mode 100644
index 23b848c3b..000000000
--- a/alg/teca_ar_detect.cxx
+++ /dev/null
@@ -1,1347 +0,0 @@
-#include "teca_ar_detect.h"
-
-#include "teca_cartesian_mesh.h"
-#include "teca_variant_array.h"
-#include "teca_table.h"
-#include "teca_calendar.h"
-#include "teca_coordinate_util.h"
-
-#include
-#include
-#include
-#include
-#include
-
-#if defined(TECA_HAS_BOOST)
-#include
-#endif
-
-using std::cerr;
-using std::endl;
-
-#define TECA_DEBUG 0
-#if TECA_DEBUG > 0
-#include "teca_cartesian_mesh_writer.h"
-#include "teca_programmable_algorithm.h"
-int write_mesh(
- const const_p_teca_cartesian_mesh &mesh,
- const const_p_teca_variant_array &vapor,
- const const_p_teca_variant_array &thres,
- const const_p_teca_variant_array &ccomp,
- const const_p_teca_variant_array &lsmask,
- const std::string &base_name);
-#endif
-
-// a description of the atmospheric river
-struct atmospheric_river
-{
- atmospheric_river() :
- pe(false), length(0.0),
- min_width(0.0), max_width(0.0),
- end_top_lat(0.0), end_top_lon(0.0),
- end_bot_lat(0.0), end_bot_lon(0.0)
- {}
-
- bool pe;
- double length;
- double min_width;
- double max_width;
- double end_top_lat;
- double end_top_lon;
- double end_bot_lat;
- double end_bot_lon;
-};
-
-std::ostream &operator<<(std::ostream &os, const atmospheric_river &ar)
-{
- os << " type=" << (ar.pe ? "PE" : "AR")
- << " length=" << ar.length
- << " width=" << ar.min_width << ", " << ar.max_width
- << " bounds=" << ar.end_bot_lon << ", " << ar.end_bot_lat << ", "
- << ar.end_top_lon << ", " << ar.end_top_lat;
- return os;
-}
-
-unsigned sauf(const unsigned nrow, const unsigned ncol, unsigned int *image);
-
-bool ar_detect(
- const_p_teca_variant_array lat,
- const_p_teca_variant_array lon,
- const_p_teca_variant_array land_sea_mask,
- p_teca_unsigned_int_array con_comp,
- unsigned long n_comp,
- double river_start_lat,
- double river_start_lon,
- double river_end_lat_low,
- double river_end_lon_low,
- double river_end_lat_high,
- double river_end_lon_high,
- double percent_in_mesh,
- double river_width,
- double river_length,
- double land_threshold_low,
- double land_threshold_high,
- atmospheric_river &ar);
-
-// set locations in the output where the input array
-// has values within the low high range.
-template
-void threshold(
- const T *input, unsigned int *output,
- size_t n_vals, T low, T high)
-{
- for (size_t i = 0; i < n_vals; ++i)
- output[i] = ((input[i] >= low) && (input[i] <= high)) ? 1 : 0;
-}
-
-
-
-// --------------------------------------------------------------------------
-teca_ar_detect::teca_ar_detect() :
- water_vapor_variable("prw"),
- land_sea_mask_variable(""),
- low_water_vapor_threshold(20),
- high_water_vapor_threshold(75),
- search_lat_low(19.0),
- search_lon_low(180.0),
- search_lat_high(56.0),
- search_lon_high(250.0),
- river_start_lat_low(18.0),
- river_start_lon_low(180.0),
- river_end_lat_low(29.0),
- river_end_lon_low(233.0),
- river_end_lat_high(56.0),
- river_end_lon_high(238.0),
- percent_in_mesh(5.0),
- river_width(1250.0),
- river_length(2000.0),
- land_threshold_low(1.0),
- land_threshold_high(std::numeric_limits::max())
-{
- this->set_number_of_input_connections(1);
- this->set_number_of_output_ports(1);
-}
-
-// --------------------------------------------------------------------------
-teca_ar_detect::~teca_ar_detect()
-{}
-
-#if defined(TECA_HAS_BOOST)
-// --------------------------------------------------------------------------
-void teca_ar_detect::get_properties_description(
- const std::string &prefix, options_description &opts)
-{
- options_description ard_opts("Options for "
- + (prefix.empty()?"teca_ar_detect":prefix));
-
- ard_opts.add_options()
- TECA_POPTS_GET(std::string, prefix, water_vapor_variable,
- "name of variable containing water vapor values")
- TECA_POPTS_GET(double, prefix, low_water_vapor_threshold,
- "low water vapor threshold")
- TECA_POPTS_GET(double, prefix, high_water_vapor_threshold,
- "high water vapor threshold")
- TECA_POPTS_GET(double, prefix, search_lat_low,
- "search space low latitude")
- TECA_POPTS_GET(double, prefix, search_lon_low,
- "search space low longitude")
- TECA_POPTS_GET(double, prefix, search_lat_high,
- "search space high latitude")
- TECA_POPTS_GET(double, prefix, search_lon_high,
- "search space high longitude")
- TECA_POPTS_GET(double, prefix, river_start_lat_low,
- "latitude used to classify as AR or PE")
- TECA_POPTS_GET(double, prefix, river_start_lon_low,
- "longitude used to classify as AR or PE")
- TECA_POPTS_GET(double, prefix, river_end_lat_low,
- "CA coastal region low latitude")
- TECA_POPTS_GET(double, prefix, river_end_lon_low,
- "CA coastal region low longitude")
- TECA_POPTS_GET(double, prefix, river_end_lat_high,
- "CA coastal region high latitude")
- TECA_POPTS_GET(double, prefix, river_end_lon_high,
- "CA coastal region high longitude")
- TECA_POPTS_GET(double, prefix, percent_in_mesh,
- "size of river in relation to search space area")
- TECA_POPTS_GET(double, prefix, river_width,
- "minimum river width")
- TECA_POPTS_GET(double, prefix, river_length,
- "minimum river length")
- TECA_POPTS_GET(std::string, prefix, land_sea_mask_variable,
- "name of variable containing land-sea mask values")
- TECA_POPTS_GET(double, prefix, land_threshold_low,
- "low land value")
- TECA_POPTS_GET(double, prefix, land_threshold_high,
- "high land value")
- ;
-
- opts.add(ard_opts);
-}
-
-// --------------------------------------------------------------------------
-void teca_ar_detect::set_properties(
- const std::string &prefix, variables_map &opts)
-{
- TECA_POPTS_SET(opts, std::string, prefix, water_vapor_variable)
- TECA_POPTS_SET(opts, double, prefix, low_water_vapor_threshold)
- TECA_POPTS_SET(opts, double, prefix, high_water_vapor_threshold)
- TECA_POPTS_SET(opts, double, prefix, search_lat_low)
- TECA_POPTS_SET(opts, double, prefix, search_lon_low)
- TECA_POPTS_SET(opts, double, prefix, search_lat_high)
- TECA_POPTS_SET(opts, double, prefix, search_lon_high)
- TECA_POPTS_SET(opts, double, prefix, river_start_lat_low)
- TECA_POPTS_SET(opts, double, prefix, river_start_lon_low)
- TECA_POPTS_SET(opts, double, prefix, river_end_lat_low)
- TECA_POPTS_SET(opts, double, prefix, river_end_lon_low)
- TECA_POPTS_SET(opts, double, prefix, river_end_lat_high)
- TECA_POPTS_SET(opts, double, prefix, river_end_lon_high)
- TECA_POPTS_SET(opts, double, prefix, percent_in_mesh)
- TECA_POPTS_SET(opts, double, prefix, river_width)
- TECA_POPTS_SET(opts, double, prefix, river_length)
- TECA_POPTS_SET(opts, std::string, prefix, land_sea_mask_variable)
- TECA_POPTS_SET(opts, double, prefix, land_threshold_low)
- TECA_POPTS_SET(opts, double, prefix, land_threshold_high)
-}
-
-#endif
-
-// --------------------------------------------------------------------------
-teca_metadata teca_ar_detect::get_output_metadata(
- unsigned int port,
- const std::vector &input_md)
-{
-#if TECA_DEBUG > 1
- cerr << teca_parallel_id()
- << "teca_ar_detect::get_output_metadata" << endl;
-#endif
- (void)port;
-
- teca_metadata output_md(input_md[0]);
- return output_md;
-}
-
-// --------------------------------------------------------------------------
-std::vector teca_ar_detect::get_upstream_request(
- unsigned int port,
- const std::vector &input_md,
- const teca_metadata &request)
-{
-#if TECA_DEBUG > 1
- cerr << teca_parallel_id()
- << "teca_ar_detect::get_upstream_request" << endl;
-#endif
- (void)port;
-
- std::vector up_reqs;
-
- teca_metadata md = input_md[0];
-
- // locate the extents of the user supplied region of
- // interest
- teca_metadata coords;
- if (md.get("coordinates", coords))
- {
- TECA_ERROR("metadata is missing \"coordinates\"")
- return up_reqs;
- }
-
- p_teca_variant_array lat;
- p_teca_variant_array lon;
- if (!(lat = coords.get("y")) || !(lon = coords.get("x")))
- {
- TECA_ERROR("metadata missing lat lon coordinates")
- return up_reqs;
- }
-
- std::vector bounds = {this->search_lon_low,
- this->search_lon_high, this->search_lat_low,
- this->search_lat_high, 0.0, 0.0};
-
- // build the request
- std::vector arrays;
- request.get("arrays", arrays);
- arrays.push_back(this->water_vapor_variable);
- if (!this->land_sea_mask_variable.empty())
- arrays.push_back(this->land_sea_mask_variable);
-
- teca_metadata up_req(request);
- up_req.set("arrays", arrays);
- up_req.set("bounds", bounds);
-
- up_reqs.push_back(up_req);
- return up_reqs;
-}
-
-// --------------------------------------------------------------------------
-const_p_teca_dataset teca_ar_detect::execute(
- unsigned int port,
- const std::vector &input_data,
- const teca_metadata &request)
-{
-#if TECA_DEBUG > 1
- cerr << teca_parallel_id() << "teca_ar_detect::execute";
- this->to_stream(cerr);
- cerr << endl;
-#endif
- (void)port;
- (void)request;
-
- // get the input dataset
- const_p_teca_cartesian_mesh mesh
- = std::dynamic_pointer_cast(input_data[0]);
- if (!mesh)
- {
- TECA_ERROR("invalid input. teca_cartesian_mesh is required")
- return nullptr;
- }
-
- // get coordinate arrays
- const_p_teca_variant_array lat = mesh->get_y_coordinates();
- const_p_teca_variant_array lon = mesh->get_x_coordinates();
-
- if (!lon || !lat)
- {
- TECA_ERROR("invalid mesh. missing lat lon coordinates")
- return nullptr;
- }
-
- // get land sea mask
- const_p_teca_variant_array land_sea_mask;
- if (this->land_sea_mask_variable.empty() ||
- !(land_sea_mask = mesh->get_point_arrays()->get(this->land_sea_mask_variable)))
- {
- // input doesn't have it, generate a stand in such
- // that land fall criteria will evaluate true
- size_t n = lat->size()*lon->size();
- p_teca_double_array lsm = teca_double_array::New(n, this->land_threshold_low);
- land_sea_mask = lsm;
- }
-
- // get the mesh extents
- std::vector extent;
- mesh->get_extent(extent);
-
- unsigned long num_rows = extent[3] - extent[2] + 1;
- unsigned long num_cols = extent[1] - extent[0] + 1;
- unsigned long num_rc = num_rows*num_cols;
-
- // get water vapor data
- const_p_teca_variant_array water_vapor
- = mesh->get_point_arrays()->get(this->water_vapor_variable);
-
- if (!water_vapor)
- {
- TECA_ERROR(
- << "Dataset missing water vapor variable \""
- << this->water_vapor_variable << "\"")
- return nullptr;
- }
-
- p_teca_table event = teca_table::New();
-
- event->declare_columns(
- "time", double(), "time_step", long(),
- "length", double(), "min width", double(),
- "max width", double(), "end_top_lat", double(),
- "end_top_lon", double(), "end_bot_lat", double(),
- "end_bot_lon", double(), "type", std::string());
-
- // get calendar
- std::string calendar;
- mesh->get_calendar(calendar);
- event->set_calendar(calendar);
-
- // get units
- std::string units;
- mesh->get_time_units(units);
- event->set_time_units(units);
-
- // get time step
- unsigned long time_step;
- mesh->get_time_step(time_step);
-
- // get offset of the current timestep
- double time = 0.0;
- mesh->get_time(time);
-
- TEMPLATE_DISPATCH(
- const teca_variant_array_impl,
- water_vapor.get(),
-
- const NT *p_wv = dynamic_cast(water_vapor.get())->get();
-
- // threshold
- p_teca_unsigned_int_array con_comp
- = teca_unsigned_int_array::New(num_rc, 0);
-
- unsigned int *p_con_comp = con_comp->get();
-
- threshold(p_wv, p_con_comp, num_rc,
- static_cast(this->low_water_vapor_threshold),
- static_cast(this->high_water_vapor_threshold));
-
-#if TECA_DEBUG > 0
- p_teca_variant_array thresh = con_comp->new_copy();
-#endif
-
- // label
- int num_comp = sauf(num_rows, num_cols, p_con_comp);
-
-#if TECA_DEBUG > 0
- write_mesh(mesh, water_vapor, thresh, con_comp,
- land_sea_mask, "ar_mesh_%t%.%e%");
-#endif
-
- // detect ar
- atmospheric_river ar;
- if (num_comp &&
- ar_detect(lat, lon, land_sea_mask, con_comp, num_comp,
- this->river_start_lat_low, this->river_start_lon_low,
- this->river_end_lat_low, this->river_end_lon_low,
- this->river_end_lat_high, this->river_end_lon_high,
- this->percent_in_mesh, this->river_width,
- this->river_length, this->land_threshold_low,
- this->land_threshold_high, ar))
- {
-#if TECA_DEBUG > 0
- cerr << teca_parallel_id() << " event detected " << time_step << endl;
-#endif
- event << time << time_step
- << ar.length << ar.min_width << ar.max_width
- << ar.end_top_lat << ar.end_top_lon
- << ar.end_bot_lat << ar.end_bot_lon
- << std::string(ar.pe ? "PE" : "AR");
- }
- )
-
- return event;
-}
-
-// --------------------------------------------------------------------------
-void teca_ar_detect::to_stream(std::ostream &os) const
-{
- os << " water_vapor_variable=" << this->water_vapor_variable
- << " land_sea_mask_variable=" << this->land_sea_mask_variable
- << " low_water_vapor_threshold=" << this->low_water_vapor_threshold
- << " high_water_vapor_threshold=" << this->high_water_vapor_threshold
- << " river_start_lon_low=" << this->river_start_lon_low
- << " river_start_lat_low=" << this->river_start_lat_low
- << " river_end_lon_low=" << this->river_end_lon_low
- << " river_end_lat_low=" << this->river_end_lat_low
- << " river_end_lon_high=" << this->river_end_lon_high
- << " river_end_lat_high=" << this->river_end_lat_high
- << " percent_in_mesh=" << this->percent_in_mesh
- << " river_width=" << this->river_width
- << " river_length=" << this->river_length
- << " land_threshodl_low=" << this->land_threshold_low
- << " land_threshodl_high=" << this->land_threshold_high;
-}
-
-
-// Code borrowed from John Wu's sauf.cpp
-// Find the minimal value starting @arg ind.
-inline unsigned afind(const std::vector& equiv,
- const unsigned ind)
-{
- unsigned ret = ind;
- while (equiv[ret] < ret)
- {
- ret = equiv[ret];
- }
- return ret;
-}
-
-// Set the values starting with @arg ind.
-inline void aset(std::vector& equiv,
- const unsigned ind, const unsigned val)
-{
- unsigned i = ind;
- while (equiv[i] < i)
- {
- unsigned j = equiv[i];
- equiv[i] = val;
- i = j;
- }
- equiv[i] = val;
-}
-
-/*
-* Purpose: Scan with Array-based Union-Find
-* Return vals: Number of connected components
-* Description: SAUF -- Scan with Array-based Union-Find.
-* This is an implementation that follows the decision try to minimize
-* number of neighbors visited and uses the array-based union-find
-* algorithms to minimize work on the union-find data structure. It works
-* with each pixel/cell of the 2D binary image individually.
-* The 2D binary image is passed to sauf as a unsigned*. On input, the
-* zero value is treated as the background, and non-zero is treated as
-* object. On successful completion of this function, the non-zero values
-* in array image is replaced by its label.
-* The return value is the number of components found.
-*/
-unsigned sauf(const unsigned nrow, const unsigned ncol, unsigned *image)
-{
- const unsigned ncells = ncol * nrow;
- const unsigned ncp1 = ncol + 1;
- const unsigned ncm1 = ncol - 1;
- std::vector equiv; // equivalence array
- unsigned nextLabel = 1;
-
- equiv.reserve(ncol);
- equiv.push_back(0);
-
- // the first cell of the first line
- if (*image != 0)
- {
- *image = nextLabel;
- equiv.push_back(nextLabel);
- ++ nextLabel;
- }
- // first row of cells
- for (unsigned i = 1; i < ncol; ++ i)
- {
- if (image[i] != 0)
- {
- if (image[i-1] != 0)
- {
- image[i] = image[i-1];
- }
- else
- {
- equiv.push_back(nextLabel);
- image[i] = nextLabel;
- ++ nextLabel;
- }
- }
- }
-
- // scan the rest of lines, check neighbor b first
- for (unsigned j = ncol; j < ncells; j += ncol)
- {
- unsigned nc, nd, k, l;
-
- // the first point of the line has two neighbors, and the two
- // neighbors must have at most one label (recorded as nc)
- if (image[j] != 0)
- {
- if (image[j-ncm1] != 0)
- nc = image[j-ncm1];
- else if (image[j-ncol] != 0)
- nc = image[j-ncol];
- else
- nc = nextLabel;
- if (nc != nextLabel) { // use existing label
- nc = equiv[nc];
- image[j] = nc;
- }
- else { // need a new label
- equiv.push_back(nc);
- image[j] = nc;
- ++ nextLabel;
- }
- }
-
- // the rest of the line
- for (unsigned i = j+1; i < j+ncol; ++ i)
- {
- if (image[i] != 0) {
- if (image[i-ncol] != 0) {
- nc = image[i-ncol];
- l = afind(equiv, nc);
- aset(equiv, nc, l);
- image[i] = l;
- }
- else if (i-ncm1 "
- << equiv[equiv[i]] << std::endl;
-#endif
- equiv[i] = equiv[equiv[i]];
- }
- else { // change to the next smallest unused label
-#if defined(_DEBUG) || defined(DEBUG)
- std::cout << i << " final " << equiv[i] << " ==> "
- << nextLabel << std::endl;
-#endif
- equiv[i] = nextLabel;
- ++ nextLabel;
- }
- }
-
- if (nextLabel < nequiv) {// relabel all cells to their minimal labels
- for (unsigned i = 0; i < ncells; ++ i)
- image[i] = equiv[image[i]];
- }
-
-#if defined(_DEBUG) || defined(DEBUG)
- std::cout << "sauf(" << nrow << ", " << ncol << ") assigned "
- << nextLabel-1 << " label" << (nextLabel>2?"s":"")
- << ", used " << nequiv << " provisional labels"
- << std::endl;
-#endif
- return nextLabel-1;
-}
-
-// do any of the detected points meet the river start
-// criteria. retrun true if so.
-template
-bool river_start_criteria_lat(
- const std::vector &con_comp_r,
- const T *p_lat,
- T river_start_lat)
-{
- unsigned long n = con_comp_r.size();
- for (unsigned long q = 0; q < n; ++q)
- {
- if (p_lat[con_comp_r[q]] >= river_start_lat)
- return true;
- }
- return false;
-}
-
-// do any of the detected points meet the river start
-// criteria. retrun true if so.
-template
-bool river_start_criteria_lon(
- const std::vector &con_comp_c,
- const T *p_lon,
- T river_start_lon)
-{
- unsigned long n = con_comp_c.size();
- for (unsigned long q = 0; q < n; ++q)
- {
- if (p_lon[con_comp_c[q]] >= river_start_lon)
- return true;
- }
- return false;
-}
-
-// helper return true if the start criteria is
-// met, and classifies the ar as PE if it starts
-// in the bottom boundary.
-template
-bool river_start_criteria(
- const std::vector &con_comp_r,
- const std::vector &con_comp_c,
- const T *p_lat,
- const T *p_lon,
- T start_lat,
- T start_lon,
- atmospheric_river &ar)
-{
- return
- ((ar.pe = river_start_criteria_lat(con_comp_r, p_lat, start_lat))
- || river_start_criteria_lon(con_comp_c, p_lon, start_lon));
-}
-
-// do any of the detected points meet the river end
-// criteria? (ie. does it hit the west coasts?) if so
-// store a bounding box covering the river and return
-// true.
-template
-bool river_end_criteria(
- const std::vector &con_comp_r,
- const std::vector &con_comp_c,
- const T *p_lat,
- const T *p_lon,
- T river_end_lat_low,
- T river_end_lon_low,
- T river_end_lat_high,
- T river_end_lon_high,
- atmospheric_river &ar)
-{
- bool end_criteria = false;
-
- std::vector end_col_idx;
-
- unsigned int count = con_comp_r.size();
- for (unsigned int i = 0; i < count; ++i)
- {
- // approximate land mask boundaries for the western coast of the US,
- T lon_val = p_lon[con_comp_c[i]];
- if ((lon_val >= river_end_lon_low) && (lon_val <= river_end_lon_high))
- end_col_idx.push_back(i);
- }
-
- // look for rows that fall between lat boundaries
- T top_lat = T();
- T top_lon = T();
- T bot_lat = T();
- T bot_lon = T();
-
- bool top_touch = false;
- unsigned int end_col_count = end_col_idx.size();
- for (unsigned int i = 0; i < end_col_count; ++i)
- {
- // approximate land mask boundaries for the western coast of the US,
- T lat_val = p_lat[con_comp_r[end_col_idx[i]]];
- if ((lat_val >= river_end_lat_low) && (lat_val <= river_end_lat_high))
- {
- T lon_val = p_lon[con_comp_c[end_col_idx[i]]];
- end_criteria = true;
- if (!top_touch)
- {
- top_touch = true;
- top_lat = lat_val;
- top_lon = lon_val;
- }
- bot_lat = lat_val;
- bot_lon = lon_val;
- }
- }
-
- ar.end_top_lat = top_lat;
- ar.end_top_lon = top_lon;
- ar.end_bot_lat = bot_lat;
- ar.end_bot_lon = bot_lon;
-
- return end_criteria;
-}
-
-/*
-* Calculate geodesic distance between two lat, long pairs
-* CODE borrowed from: from http://www.geodatasource.com/developers/c
-* from http://osiris.tuwien.ac.at/~wgarn/gis-gps/latlong.html
-* from http://www.codeproject.com/KB/cpp/Distancecplusplus.aspx
-*/
-template
-T geodesic_distance(T lat1, T lon1, T lat2, T lon2)
-{
- T deg_to_rad = T(M_PI/180.0);
-
- T dlat1 = lat1*deg_to_rad;
- T dlon1 = lon1*deg_to_rad;
- T dlat2 = lat2*deg_to_rad;
- T dlon2 = lon2*deg_to_rad;
-
- T dLon = dlon1 - dlon2;
- T dLat = dlat1 - dlat2;
-
- T sin_dLat_half_sq = sin(dLat/T(2.0));
- sin_dLat_half_sq *= sin_dLat_half_sq;
-
- T sin_dLon_half_sq = sin(dLon/T(2.0));
- sin_dLon_half_sq *= sin_dLon_half_sq;
-
- T aHarv = sin_dLat_half_sq
- + cos(dlat1)*cos(dlat2)*sin_dLon_half_sq;
-
- T cHarv = T(2.0)*atan2(sqrt(aHarv), sqrt(T(1.0) - aHarv));
-
- T R = T(6371.0);
- T distance = R*cHarv;
-
- return distance;
-}
-
-/*
-* This function calculates the average geodesic width
-* As each pixel represents certain area, the total area
-* is the product of the number of pixels and the area of
-* one pixel. The average width is: the total area divided
-* by the medial axis length
-* We are calculating the average width, since we are not
-* determining where exactly to cut off the tropical region
-* to calculate the real width of an atmospheric river
-*/
-template
-T avg_width(
- const std::vector &con_comp_r,
- const std::vector &con_comp_c,
- T ar_len,
- const T *p_lat,
- const T *p_lon)
-{
-/*
- // TODO -- need bounds checking when doing things like
- // p_lat[con_comp_r[0] + 1]. also because it's potentially
- // a stretched cartesian mesh need to compute area of
- // individual cells
-
- // length of cell in lat direction
- T lat_val[2] = {p_lat[con_comp_r[0]], p_lat[con_comp_r[0] + 1]};
- T lon_val[2] = {p_lon[con_comp_c[0]], p_lon[con_comp_c[0]]};
- T dlat = geodesic_distance(lat_val[0], lon_val[0], lat_val[1], lon_val[1]);
-
- // length of cell in lon direction
- lat_val[1] = lat_val[0];
- lon_val[1] = p_lon[con_comp_c[0] + 1];
- T dlon = geodesic_distance(lat_val[0], lon_val[0], lat_val[1], lon_val[1]);
-*/
- (void)con_comp_c;
- // compute area of the first cell in the input mesh
- // length of cell in lat direction
- T lat_val[2] = {p_lat[0], p_lat[1]};
- T lon_val[2] = {p_lon[0], p_lon[0]};
- T dlat = geodesic_distance(lat_val[0], lon_val[0], lat_val[1], lon_val[1]);
-
- // length of cell in lon direction
- lat_val[1] = lat_val[0];
- lon_val[1] = p_lon[1];
- T dlon = geodesic_distance(lat_val[0], lon_val[0], lat_val[1], lon_val[1]);
-
- // area
- T pixel_area = dlat*dlon;
- T total_area = pixel_area*con_comp_r.size();
-
- // avg width
- T avg_width = total_area/ar_len;
- return avg_width;
-}
-
-/*
-* Find the middle point between two pairs of lat and lon values
-* http://stackoverflow.com/questions/4164830/geographic-midpoint-between-two-coordinates
-*/
-template
-void geodesic_midpoint(T lat1, T lon1, T lat2, T lon2, T &mid_lat, T &mid_lon)
-{
- T deg_to_rad = T(M_PI/180.0);
- T dLon = (lon2 - lon1) * deg_to_rad;
- T dLat1 = lat1 * deg_to_rad;
- T dLat2 = lat2 * deg_to_rad;
- T dLon1 = lon1 * deg_to_rad;
-
- T Bx = cos(dLat2) * cos(dLon);
- T By = cos(dLat2) * sin(dLon);
-
- mid_lat = atan2(sin(dLat1)+sin(dLat2),
- sqrt((cos(dLat1)+Bx)*(cos(dLat1)+Bx)+By*By));
-
- mid_lon = dLon1 + atan2(By, (cos(dLat1)+Bx));
-
- T rad_to_deg = T(180.0/M_PI);
- mid_lat *= rad_to_deg;
- mid_lon *= rad_to_deg;
-}
-
-/*
-* Find the length along the medial axis of a connected component
-* Medial length is the sum of the distances between the medial
-* points in the connected component
-*/
-template
-T medial_length(
- const std::vector &con_comp_r,
- const std::vector &con_comp_c,
- const T *p_lat,
- const T *p_lon)
-{
- std::vector jb_r1;
- std::vector jb_c1;
- std::vector jb_c2;
-
- long row_track = -1;
-
- unsigned long count = con_comp_r.size();
- for (unsigned long i = 0; i < count; ++i)
- {
- if (row_track != con_comp_r[i])
- {
- jb_r1.push_back(con_comp_r[i]);
- jb_c1.push_back(con_comp_c[i]);
-
- jb_c2.push_back(con_comp_c[i]);
-
- row_track = con_comp_r[i];
- }
- else
- {
- jb_c2.back() = con_comp_c[i];
- }
- }
-
- T total_dist = T();
-
- long b_count = jb_r1.size() - 1;
- for (long i = 0; i < b_count; ++i)
- {
- T lat_val[2];
- T lon_val[2];
-
- lat_val[0] = p_lat[jb_r1[i]];
- lat_val[1] = p_lat[jb_r1[i]];
-
- lon_val[0] = p_lon[jb_c1[i]];
- lon_val[1] = p_lon[jb_c2[i]];
-
- T mid_lat1;
- T mid_lon1;
-
- geodesic_midpoint(
- lat_val[0], lon_val[0], lat_val[1], lon_val[1],
- mid_lat1, mid_lon1);
-
- lat_val[0] = p_lat[jb_r1[i+1]];
- lat_val[1] = p_lat[jb_r1[i+1]];
-
- lon_val[0] = p_lon[jb_c1[i+1]];
- lon_val[1] = p_lon[jb_c2[i+1]];
-
- T mid_lat2;
- T mid_lon2;
-
- geodesic_midpoint(
- lat_val[0], lon_val[0], lat_val[1], lon_val[1],
- mid_lat2, mid_lon2);
-
- total_dist
- += geodesic_distance(mid_lat1, mid_lon1, mid_lat2, mid_lon2);
- }
-
- return total_dist;
-}
-
-/*
-// Suren's function
-// helper return true if the geometric conditions
-// on an ar are satisfied. also stores the length
-// and width of the river.
-template
-bool river_geometric_criteria(
- const std::vector &con_comp_r,
- const std::vector &con_comp_c,
- const T *p_lat,
- const T *p_lon,
- double river_length,
- double river_width,
- atmospheric_river &ar)
-{
- ar.length = medial_length(con_comp_r, con_comp_c, p_lat, p_lon);
-
- ar.width = avg_width(con_comp_r, con_comp_c,
- static_cast(ar.length), p_lat, p_lon);
-
- return (ar.length >= river_length) && (ar.width <= river_width);
-}
-*/
-
-// Junmin's function for height of a triangle
-template
-T triangle_height(T base, T s1, T s2)
-{
- // area from Heron's fomula
- T p = (base + s1 + s2)/T(2);
- T area = p*(p - base)*(p - s1)*(p - s2);
- // detect impossible triangle
- if (area < T())
- return std::min(s1, s2);
- // height from A = 1/2 b h
- return T(2)*sqrt(area)/base;
-}
-
-// TDataProcessor::check_geodesic_width_top_down
-// Junmin's function for detecting river based on
-// it's geometric properties
-template
-bool river_geometric_criteria(
- const std::vector &con_comp_r,
- const std::vector &con_comp_c,
- const T *p_lat,
- const T *p_lon,
- double river_length,
- double river_width,
- atmospheric_river &ar)
-{
- std::vector distinct_rows;
- std::vector leftmost_col;
- std::vector rightmost_col;
-
- int row_track = -1;
- size_t count = con_comp_r.size();
- for (size_t i = 0; i < count; ++i)
- {
- if (row_track != con_comp_r[i])
- {
- row_track = con_comp_r[i];
-
- distinct_rows.push_back(con_comp_r[i]);
- leftmost_col.push_back(con_comp_c[i]);
- rightmost_col.push_back(con_comp_c[i]);
- }
- else
- {
- rightmost_col.back() = con_comp_c[i];
- }
- }
-
- // river metrics
- T length_from_top = T();
- T min_width = std::numeric_limits::max();
- T max_width = std::numeric_limits::lowest();
-
- for (long i = distinct_rows.size() - 2; i >= 0; --i)
- {
- // for each row w respect to row above it. triangulate
- // a quadrilateral composed of left and right most points
- // in this and the above rows. ccw ordering from lower
- // left corner is A,B,D,C.
-
- // low left-right distance
- T AB = geodesic_distance(
- p_lat[distinct_rows[i]], p_lon[leftmost_col[i]],
- p_lat[distinct_rows[i]], p_lon[rightmost_col[i]]);
-
- // left side bottom-top distance
- T AC = geodesic_distance(
- p_lat[distinct_rows[i]], p_lon[leftmost_col[i]],
- p_lat[distinct_rows[i+1]], p_lon[leftmost_col[i]]);
-
- // distance from top left to bottom right, across
- T BC = geodesic_distance(
- p_lat[distinct_rows[i]], p_lon[rightmost_col[i]],
- p_lat[distinct_rows[i+1]], p_lon[leftmost_col[i+1]]);
-
- // high left-right distance
- T CD = geodesic_distance(
- p_lat[distinct_rows[i+1]], p_lon[leftmost_col[i+1]],
- p_lat[distinct_rows[i+1]], p_lon[rightmost_col[i+1]]);
-
- // right side bottom-top distance
- T BD = geodesic_distance(
- p_lat[distinct_rows[i]], p_lon[rightmost_col[i]],
- p_lat[distinct_rows[i+1]], p_lon[rightmost_col[i+1]]);
-
- T height_from_b = triangle_height(AC, AB, BC);
- T height_from_c = triangle_height(BD, BC, CD);
-
- T curr_min = std::min(height_from_b, height_from_c);
-
- // test width criteria
- if (curr_min > river_width)
- {
- // TODO -- first time through the loop length == 0. is it intentional
- // to discard the detection or should length calc take place before this test?
- // note: first time through loop none of the event details have been recoreded
- if (length_from_top <= river_length)
- {
- // too short to be a river
- return false;
- }
- else
- {
- // part of a connected region is AR
- ar.min_width = static_cast(min_width);
- ar.max_width = static_cast(max_width);
- ar.length = static_cast(length_from_top);
- return true;
- }
- }
-
- // update width
- min_width = std::min(min_width, curr_min);
- max_width = std::max(max_width, curr_min);
-
- // update length
- T mid_bot_lat;
- T mid_bot_lon;
- geodesic_midpoint(
- p_lat[distinct_rows[i]], p_lon[leftmost_col[i]],
- p_lat[distinct_rows[i]], p_lon[rightmost_col[i]],
- mid_bot_lat, mid_bot_lon);
-
- T mid_top_lat;
- T mid_top_lon;
- geodesic_midpoint(
- p_lat[distinct_rows[i+1]], p_lon[leftmost_col[i+1]],
- p_lat[distinct_rows[i+1]], p_lon[rightmost_col[i+1]],
- mid_top_lat, mid_top_lon);
-
- length_from_top += geodesic_distance(
- mid_bot_lat, mid_bot_lon, mid_top_lat, mid_top_lon);
- }
-
- // check the length criteria.
- // TODO: if we are here the widtrh critera was not met
- // so the following detection is based solely on the length?
- if (length_from_top > river_length)
- {
- // AR
- ar.min_width = static_cast(min_width);
- ar.max_width = static_cast(max_width);
- ar.length = static_cast(length_from_top);
- return true;
- }
-
- return false;
-}
-
-
-
-// Junmin's function checkRightBoundary
-// note: if land sea mask is not available land array
-// must all be true.
-template
-bool river_end_criteria(
- const std::vector &con_comp_r,
- const std::vector &con_comp_c,
- const std::vector &land,
- const T *p_lat,
- const T *p_lon,
- T river_end_lat_low,
- T river_end_lon_low,
- T river_end_lat_high,
- T river_end_lon_high,
- atmospheric_river &ar)
-{
- // locate component points within shoreline
- // box
- bool first_crossing = false;
- bool event_detected = false;
-
- T top_lat = T();
- T top_lon = T();
- T bot_lat = T();
- T bot_lon = T();
-
- std::vector right_bound_col_idx;
- size_t count = con_comp_c.size();
- for (size_t i = 0; i < count; ++i)
- {
- T lat = p_lat[con_comp_r[i]];
- T lon = p_lon[con_comp_c[i]];
-
- if ((lat >= river_end_lat_low) && (lat <= river_end_lat_high)
- && (lon >= river_end_lon_low) && (lon <= river_end_lon_high))
- {
- if (!event_detected)
- event_detected = land[i];
-
- if (!first_crossing)
- {
- first_crossing = true;
- top_lat = lat;
- top_lon = lon;
- }
- bot_lat = lat;
- bot_lon = lon;
- }
- }
-
- ar.end_top_lat = top_lat;
- ar.end_top_lon = top_lon;
- ar.end_bot_lat = bot_lat;
- ar.end_bot_lon = bot_lon;
-
- return event_detected;
-}
-
-// Junmin's function
-template
-void classify_event(
- const std::vector &con_comp_r,
- const std::vector &con_comp_c,
- const T *p_lat,
- const T *p_lon,
- T start_lat,
- T start_lon,
- atmospheric_river &ar)
-{
- // classification determined by first detected point in event
- // is closer to left or to bottom
- T lat = p_lat[con_comp_r[0]];
- T lon = p_lon[con_comp_c[0]];
-
- ar.pe = false;
- if ((lat - start_lat) < (lon - start_lon))
- ar.pe = true; // PE
-}
-
-/*
-* The main function that checks whether an AR event exists in
-* given sub-plane of data. This currently applies only to the Western
-* coast of the USA. return true if an ar is found.
-*/
-bool ar_detect(
- const_p_teca_variant_array lat,
- const_p_teca_variant_array lon,
- const_p_teca_variant_array land_sea_mask,
- p_teca_unsigned_int_array con_comp,
- unsigned long n_comp,
- double river_start_lat,
- double river_start_lon,
- double river_end_lat_low,
- double river_end_lon_low,
- double river_end_lat_high,
- double river_end_lon_high,
- double percent_in_mesh,
- double river_width,
- double river_length,
- double land_threshold_low,
- double land_threshold_high,
- atmospheric_river &ar)
-{
- NESTED_TEMPLATE_DISPATCH_FP(
- const teca_variant_array_impl,
- lat.get(),
- 1,
-
- NESTED_TEMPLATE_DISPATCH(
- const teca_variant_array_impl,
- land_sea_mask.get(),
- 2,
-
- const NT1 *p_lat = dynamic_cast(lat.get())->get();
- const NT1 *p_lon = dynamic_cast(lon.get())->get();
-
- const NT2 *p_land_sea_mask
- = dynamic_cast(land_sea_mask.get())->get();
-
- NT1 start_lat = static_cast(river_start_lat);
- NT1 start_lon = static_cast(river_start_lon);
- NT1 end_lat_low = static_cast(river_end_lat_low);
- NT1 end_lon_low = static_cast(river_end_lon_low);
- NT1 end_lat_high = static_cast(river_end_lat_high);
- NT1 end_lon_high = static_cast(river_end_lon_high);
-
- unsigned long num_rows = lat->size();
- unsigned long num_cols = lon->size();
-
- // # in PE is % of points in regious mesh
- unsigned long num_rc = num_rows*num_cols;
- unsigned long thr_count = num_rc*percent_in_mesh/100.0;
-
- unsigned int *p_labels = con_comp->get();
- for (unsigned int i = 1; i <= n_comp; ++i)
- {
- // for all discrete connected component labels
- // verify if there exists an AR
- std::vector con_comp_r;
- std::vector con_comp_c;
- std::vector land;
-
- for (unsigned long r = 0, q = 0; r < num_rows; ++r)
- {
- for (unsigned long c = 0; c < num_cols; ++c, ++q)
- {
- if (p_labels[q] == i)
- {
- // gather points of this connected component
- con_comp_r.push_back(r);
- con_comp_c.push_back(c);
-
- // identify them as land or not
- land.push_back(
- (p_land_sea_mask[q] >= land_threshold_low)
- && (p_land_sea_mask[q] < land_threshold_high));
- }
- }
- }
-
- // check for ar criteria
- unsigned long count = con_comp_r.size();
- if ((count > thr_count)
- && river_end_criteria(
- con_comp_r, con_comp_c, land,
- p_lat, p_lon,
- end_lat_low, end_lon_low,
- end_lat_high, end_lon_high,
- ar)
- && river_geometric_criteria(
- con_comp_r, con_comp_c, p_lat, p_lon,
- river_length, river_width, ar))
- {
- // determine if PE or AR
- classify_event(
- con_comp_r, con_comp_c, p_lat, p_lon,
- start_lat, start_lon, ar);
- return true;
- }
- }
- )
- )
- return false;
-}
-
-#if TECA_DEBUG > 0
-// helper to dump a dataset for debugging
-int write_mesh(
- const const_p_teca_cartesian_mesh &mesh,
- const const_p_teca_variant_array &vapor,
- const const_p_teca_variant_array &thresh,
- const const_p_teca_variant_array &ccomp,
- const const_p_teca_variant_array &lsmask,
- const std::string &file_name)
-{
- p_teca_cartesian_mesh m = teca_cartesian_mesh::New();
- m->copy_metadata(mesh);
-
- p_teca_array_collection pac = m->get_point_arrays();
- pac->append("vapor", std::const_pointer_cast(vapor));
- pac->append("thresh", std::const_pointer_cast(thresh));
- pac->append("ccomp", std::const_pointer_cast(ccomp));
- pac->append("lsmask", std::const_pointer_cast(lsmask));
-
- p_teca_programmable_algorithm s = teca_programmable_algorithm::New();
- s->set_name("serve_mesh");
- s->set_number_of_input_connections(0);
- s->set_number_of_output_ports(1);
- s->set_execute_callback(
- [m] (unsigned int, const std::vector &,
- const teca_metadata &) -> const_p_teca_dataset { return m; }
- );
-
- p_teca_cartesian_mesh_writer w
- = teca_cartesian_mesh_writer::New();
-
- w->set_file_name(file_name);
- w->set_input_connection(s->get_output_port());
- w->update();
-
- return 0;
-}
-#endif
diff --git a/alg/teca_ar_detect.h b/alg/teca_ar_detect.h
deleted file mode 100644
index 11ef0f199..000000000
--- a/alg/teca_ar_detect.h
+++ /dev/null
@@ -1,136 +0,0 @@
-#ifndef teca_ar_detect_h
-#define teca_ar_detect_h
-
-#include "teca_shared_object.h"
-#include "teca_algorithm.h"
-#include "teca_metadata.h"
-
-#include
-#include
-
-TECA_SHARED_OBJECT_FORWARD_DECL(teca_ar_detect)
-
-/**
-Suren and Junmin's atmospheric river detector.
-
-The algorithm searches for atmospheric rivers that
-end on the California coast in water vapor data over
-a specific subset of the input data. A river is detected
-based on it's length, width, and percent area of the
-search space. The algorithm can optionally use a
-land-sea mask to increase accuracy of the California
-coast. Without the land-sea mask a box is used.
-*/
-class teca_ar_detect : public teca_algorithm
-{
-public:
- TECA_ALGORITHM_STATIC_NEW(teca_ar_detect)
- TECA_ALGORITHM_DELETE_COPY_ASSIGN(teca_ar_detect)
- TECA_ALGORITHM_CLASS_NAME(teca_ar_detect)
- ~teca_ar_detect();
-
- // report/initialize to/from Boost program options
- // objects.
- TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION()
- TECA_SET_ALGORITHM_PROPERTIES()
-
- // set/get the name of the integrated water vapor variable
- TECA_ALGORITHM_PROPERTY(std::string, water_vapor_variable)
-
- // set/get threshold on water vapor variable used
- // to segment the data
- TECA_ALGORITHM_PROPERTY(double, low_water_vapor_threshold)
- TECA_ALGORITHM_PROPERTY(double, high_water_vapor_threshold)
-
- // set/get the region of interest in lat lon coordinate system
- // defaults are 19 56 180 250
- TECA_ALGORITHM_PROPERTY(double, search_lat_low)
- TECA_ALGORITHM_PROPERTY(double, search_lon_low)
- TECA_ALGORITHM_PROPERTY(double, search_lat_high)
- TECA_ALGORITHM_PROPERTY(double, search_lon_high)
-
- // set/get the river source region in lat lon coordinate system
- // defaults are 18 180
- TECA_ALGORITHM_PROPERTY(double, river_start_lat_low)
- TECA_ALGORITHM_PROPERTY(double, river_start_lon_low)
-
- // set/get the river ladfall region in lat lon coordinate system
- // defaults are 29 233 56 238
- TECA_ALGORITHM_PROPERTY(double, river_end_lat_low)
- TECA_ALGORITHM_PROPERTY(double, river_end_lon_low)
- TECA_ALGORITHM_PROPERTY(double, river_end_lat_high)
- TECA_ALGORITHM_PROPERTY(double, river_end_lon_high)
-
- // set/get the area as a percent of the search space that
- // a potential river must occupy
- TECA_ALGORITHM_PROPERTY(double, percent_in_mesh)
-
- // set/get the minimum river width and length. defaults
- // are 1250 2000
- TECA_ALGORITHM_PROPERTY(double, river_width)
- TECA_ALGORITHM_PROPERTY(double, river_length)
-
- // set/get the land-sea mask variable. this array
- // will be used to identify land from ocean using
- // land_threshold properties.
- TECA_ALGORITHM_PROPERTY(std::string, land_sea_mask_variable)
-
- // set/get the land classification range [low high). defaults
- // are [1.0 DOUBLE_MAX)
- TECA_ALGORITHM_PROPERTY(double, land_threshold_low)
- TECA_ALGORITHM_PROPERTY(double, land_threshold_high)
-
- // send humand readable representation to the
- // stream
- virtual void to_stream(std::ostream &os) const override;
-
-protected:
- teca_ar_detect();
-
- // helper that computes the output extent
- int get_active_extent(
- p_teca_variant_array lat,
- p_teca_variant_array lon,
- std::vector &extent) const;
-
-private:
- virtual
- teca_metadata get_output_metadata(
- unsigned int port,
- const std::vector &input_md) override;
-
- virtual
- std::vector get_upstream_request(
- unsigned int port,
- const std::vector &input_md,
- const teca_metadata &request) override;
-
- virtual
- const_p_teca_dataset execute(
- unsigned int port,
- const std::vector &input_data,
- const teca_metadata &request) override;
-
-private:
- std::string water_vapor_variable;
- std::string land_sea_mask_variable;
- double low_water_vapor_threshold;
- double high_water_vapor_threshold;
- double search_lat_low;
- double search_lon_low;
- double search_lat_high;
- double search_lon_high;
- double river_start_lat_low;
- double river_start_lon_low;
- double river_end_lat_low;
- double river_end_lon_low;
- double river_end_lat_high;
- double river_end_lon_high;
- double percent_in_mesh;
- double river_width;
- double river_length;
- double land_threshold_low;
- double land_threshold_high;
-};
-
-#endif
diff --git a/alg/teca_bayesian_ar_detect.cxx b/alg/teca_bayesian_ar_detect.cxx
index 8a7705a09..c14dd6f1f 100644
--- a/alg/teca_bayesian_ar_detect.cxx
+++ b/alg/teca_bayesian_ar_detect.cxx
@@ -608,8 +608,10 @@ void teca_bayesian_ar_detect::internals_t::clear()
teca_bayesian_ar_detect::teca_bayesian_ar_detect() :
min_component_area_variable("min_component_area"),
min_ivt_variable("min_water_vapor"),
- hwhm_latitude_variable("hwhm_latitude"), thread_pool_size(1),
- verbose(0), internals(new internals_t)
+ hwhm_latitude_variable("hwhm_latitude"),
+ ar_probability_variable("ar_probability"),
+ thread_pool_size(1),
+ internals(new internals_t)
{
this->set_number_of_input_connections(1);
this->set_number_of_output_ports(1);
@@ -631,23 +633,25 @@ void teca_bayesian_ar_detect::get_properties_description(
opts.add_options()
TECA_POPTS_GET(std::string, prefix, ivt_variable,
- "name of the water vapor variable (\"\")")
+ "Set the name of the integrated vaopr transport(IVT) variable to"
+ " compute AR probability from.")
TECA_POPTS_GET(std::string, prefix, min_component_area_variable,
- "name of the column in the parameter table containing the "
- "component area threshold (\"min_component_area\")")
+ "Set the name of the column in the parameter table containing the "
+ "minimum feature area threshold.")
TECA_POPTS_GET(std::string, prefix, min_ivt_variable,
- "name of the column in the parameter table containing the "
- "water vapor threshold (\"min_water_vapor\")")
+ "Set the name of the column in the parameter table containing the "
+ "minimum percentile IVT threshold.")
TECA_POPTS_GET(std::string, prefix, hwhm_latitude_variable,
- "name of the column in the parameter table containing the "
- "half width at half max latitude (\"hwhm_latitude\")")
+ "Set the name of the column in the parameter table containing the "
+ "half width at half max latitude mask value.")
+ TECA_POPTS_GET(std::string, prefix, ar_probability_variable,
+ "Set the name of the variable to store the computed AR probability in.")
TECA_POPTS_GET(int, prefix, thread_pool_size,
- "number of threads to parallelize execution over (1)")
- TECA_POPTS_GET(int, prefix, verbose,
- "flag indicating diagnostic info should be displayed in "
- "the terminal (0)")
+ "Set the number of threads to parallelize execution over.")
;
+ this->teca_algorithm::get_properties_description(prefix, opts);
+
global_opts.add(opts);
}
@@ -655,10 +659,13 @@ void teca_bayesian_ar_detect::get_properties_description(
void teca_bayesian_ar_detect::set_properties(const std::string &prefix,
variables_map &opts)
{
+ this->teca_algorithm::set_properties(prefix, opts);
+
TECA_POPTS_SET(opts, std::string, prefix, ivt_variable)
TECA_POPTS_SET(opts, std::string, prefix, min_component_area_variable)
TECA_POPTS_SET(opts, std::string, prefix, min_ivt_variable)
TECA_POPTS_SET(opts, std::string, prefix, hwhm_latitude_variable)
+ TECA_POPTS_SET(opts, std::string, prefix, ar_probability_variable)
TECA_POPTS_SET(opts, int, prefix, thread_pool_size)
TECA_POPTS_SET(opts, int, prefix, verbose)
}
@@ -693,7 +700,10 @@ void teca_bayesian_ar_detect::set_thread_pool_size(int n)
// --------------------------------------------------------------------------
unsigned int teca_bayesian_ar_detect::get_thread_pool_size() const noexcept
{
- return this->internals->queue->size();
+ unsigned int n_threads = 0;
+ if (this->internals->queue)
+ n_threads = this->internals->queue->size();
+ return n_threads;
}
// --------------------------------------------------------------------------
@@ -805,7 +815,7 @@ teca_metadata teca_bayesian_ar_detect::get_output_metadata(
// report the variable that we compute, for each timestep from the
// parameter tables.
teca_metadata md(input_md[0]);
- md.append("variables", std::string("ar_probability"));
+ md.append("variables", std::string(this->ar_probability_variable));
// add attributes to enable CF I/O
teca_metadata atts;
@@ -816,7 +826,7 @@ teca_metadata teca_bayesian_ar_detect::get_output_metadata(
0, "unitless", "posterior AR flag",
"the posterior probability of the presence of an atmospheric river");
- atts.set("ar_probability", (teca_metadata)prob_atts);
+ atts.set(this->ar_probability_variable, (teca_metadata)prob_atts);
unsigned long num_params =
this->internals->parameter_table->get_number_of_rows();
@@ -875,7 +885,7 @@ std::vector teca_bayesian_ar_detect::get_upstream_request(
arrays.insert(this->ivt_variable);
// remove what we produce
- arrays.erase("ar_probability");
+ arrays.erase(this->ar_probability_variable);
arrays.erase("ar_count");
arrays.erase("parameter_table_row");
@@ -1039,7 +1049,7 @@ const_p_teca_dataset teca_bayesian_ar_detect::execute(
// set up the reduction which computes the average over runs of all control
// parameter combinations provided in the parameter table
::parameter_table_reduction reduce(parameter_table_size,
- "wv_cc", "ar_probability");
+ "wv_cc", this->ar_probability_variable);
p_teca_programmable_reduce pr = teca_programmable_reduce::New();
pr->set_name("parameter_table_reduce");
diff --git a/alg/teca_bayesian_ar_detect.h b/alg/teca_bayesian_ar_detect.h
index c2805cf99..aecbd9411 100644
--- a/alg/teca_bayesian_ar_detect.h
+++ b/alg/teca_bayesian_ar_detect.h
@@ -10,31 +10,35 @@
TECA_SHARED_OBJECT_FORWARD_DECL(teca_bayesian_ar_detect)
-/// CASCADE BARD atmospheric river detector
+/// The TECA BARD atmospheric river detector.
/**
-Given a point wise IVT (integrated vapor transport) field and a training
-parameter table computes the point wise probability of an atmospheric river
-using the CASCADE BARD algorithm.
-
-Required inputs:
-
- 1. IVT (integrated vapor transport) array on a Cartesian nesh.
- 2. a compatible parameter table. columns of which are : min IVT,
- component area, HWHM lattitude
-
-The names of the input varibale and columns can be specified at run time
-through algorithm properties.
-
-Produces:
-
- A Cartesian mesh with probability of an AR stored in the point centered
- array named "ar_probability". The diagnostic quantites "ar_count" amd
- "parameter_table_row" are stored in information arrays.
-
-For more information see:
-
-Detection of Atmospheric Rivers with Inline Uncertainty Quantification: TECA-BARD v1.0
-O'Brien, T. A et al. Geoscientific Model Development, 2020
+ * Given a point wise IVT (integrated vapor transport) field and a training
+ * parameter table computes the point wise probability of an atmospheric river
+ * using the TECA BARD algorithm.
+ *
+ * Required inputs:
+ *
+ * 1. IVT (integrated vapor transport) array on a Cartesian nesh.
+ * 2. a compatible parameter table. columns of which are : min IVT,
+ * component area, HWHM lattitude
+ *
+ * The names of the input varibale and columns can be specified at run time
+ * through algorithm properties.
+ *
+ * Produces:
+ *
+ * A Cartesian mesh with probability of an AR stored in the point centered
+ * array named "ar_probability". The diagnostic quantites "ar_count" amd
+ * "parameter_table_row" are stored in information arrays.
+ *
+ * For more information see:
+ *
+ * O’Brien, T. A., Risser, M. D., Loring, B., Elbashandy, A. A., Krishnan, H.,
+ * Johnson, J., Patricola, C. M., O’Brien, J. P., Mahesh, A., Arriaga Ramirez,
+ * S., Rhoades, A. M., Charn, A., Inda DÃaz, H., & Collins, W. D. (2020).
+ * Detection of atmospheric rivers with inline uncertainty quantification:
+ * TECA-BARD v1.0.1. Geoscientific Model Development, 13(12), 6131–6148.
+ * https://doi.org/10.5194/gmd-13-6131-2020
*/
class teca_bayesian_ar_detect : public teca_algorithm
{
@@ -49,28 +53,54 @@ class teca_bayesian_ar_detect : public teca_algorithm
TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION()
TECA_SET_ALGORITHM_PROPERTIES()
- // set the name of the input array
+ /** @name ivt_variable
+ * Sets the name of the array containing the IVT field to detect ARs in.
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(std::string, ivt_variable)
+ ///@}
- // set the names of columns in the parameter table.
+ /** @name min_ivt_variable
+ * Set the names of the minimum IVT column in the parameter table.
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(std::string, min_ivt_variable)
- TECA_ALGORITHM_PROPERTY(std::string, min_component_area_variable)
- TECA_ALGORITHM_PROPERTY(std::string, hwhm_latitude_variable)
+ ///@}
- // flag indicating verbose terminal output is desired.
- // default is 0
- TECA_ALGORITHM_PROPERTY(int, verbose)
+ /** @name min_component_area_variable
+ * Set the names of the minimum area column in the parameter table.
+ */
+ ///@{
+ TECA_ALGORITHM_PROPERTY(std::string, min_component_area_variable)
+ ///@}
- // set/get the number of threads in the pool. setting
- // to -1 results in a thread per core factoring in all MPI
- // ranks running on the node. the default is -1.
+ /** @name hwhm_latitude_variable
+ * Set the names of the HWHM column in the parameter table.
+ */
+ ///@{
+ TECA_ALGORITHM_PROPERTY(std::string, hwhm_latitude_variable)
+ ///@}
+
+ /** @name probability variable
+ * Set the name of the variable to store output probability as.
+ */
+ ///@{
+ TECA_ALGORITHM_PROPERTY(std::string, ar_probability_variable)
+ ///@}
+
+ /** Set the number of threads in the pool. Setting to -1 results in a
+ * thread per core factoring in all MPI ranks running on the node. the
+ * default is -1.
+ */
void set_thread_pool_size(int n_threads);
+
+ /// Get the number of threads in the pool.
unsigned int get_thread_pool_size() const noexcept;
- // override the input connections because we are going to
- // take the first input and use it to generate metadata.
- // the second input then becomes the only one the pipeline
- // knows about.
+ /** override the input connections because we are going to take the first
+ * input and use it to generate metadata. the second input then becomes
+ * the only one the pipeline knows about.
+ */
void set_input_connection(unsigned int id,
const teca_algorithm_output_port &port) override;
@@ -98,8 +128,8 @@ class teca_bayesian_ar_detect : public teca_algorithm
std::string min_component_area_variable;
std::string min_ivt_variable;
std::string hwhm_latitude_variable;
+ std::string ar_probability_variable;
int thread_pool_size;
- int verbose;
struct internals_t;
internals_t *internals;
diff --git a/alg/teca_bayesian_ar_detect_parameters.cxx b/alg/teca_bayesian_ar_detect_parameters.cxx
index ca58d5c85..23511a112 100644
--- a/alg/teca_bayesian_ar_detect_parameters.cxx
+++ b/alg/teca_bayesian_ar_detect_parameters.cxx
@@ -973,9 +973,11 @@ void teca_bayesian_ar_detect_parameters::get_properties_description(
opts.add_options()
TECA_POPTS_GET(long, prefix, number_of_rows,
- "the number of parameter table rows to serve (-1)")
+ "the number of parameter table rows to serve")
;
+ this->teca_algorithm::get_properties_description(prefix, opts);
+
global_opts.add(opts);
}
@@ -983,6 +985,8 @@ void teca_bayesian_ar_detect_parameters::get_properties_description(
void teca_bayesian_ar_detect_parameters::set_properties(const std::string &prefix,
variables_map &opts)
{
+ this->teca_algorithm::set_properties(prefix, opts);
+
TECA_POPTS_SET(opts, long, prefix, number_of_rows)
}
#endif
diff --git a/alg/teca_bayesian_ar_detect_parameters.h b/alg/teca_bayesian_ar_detect_parameters.h
index 7f5b1dc56..518305301 100644
--- a/alg/teca_bayesian_ar_detect_parameters.h
+++ b/alg/teca_bayesian_ar_detect_parameters.h
@@ -5,10 +5,10 @@
TECA_SHARED_OBJECT_FORWARD_DECL(teca_bayesian_ar_detect_parameters)
-/**
-An algorithm that constructs and serves up the parameter
-table needed to run the Bayesain AR detector
-*/
+/** @brief
+ * An algorithm that constructs and serves up the parameter
+ * table needed to run the Bayesian AR detector.
+ */
class teca_bayesian_ar_detect_parameters : public teca_algorithm
{
public:
@@ -22,12 +22,16 @@ class teca_bayesian_ar_detect_parameters : public teca_algorithm
TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION()
TECA_SET_ALGORITHM_PROPERTIES()
- // control the number of rows coppied into the table. The rows are
- // copppied in sequential order starting from row zero. The default value
- // of -1 is used to serve all rows. See also get_parameter_table_size.
+ /** @name number_of_rows
+ * control the number of rows copied into the table. The rows are copied
+ * in sequential order starting from row zero. The default value of -1 is
+ * used to serve all rows. See also get_parameter_table_size.
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(long, number_of_rows)
+ ///@}
- // return the number of rows in the internal parameter table.
+ /// return the number of rows in the internal parameter table.
unsigned long get_parameter_table_size();
protected:
diff --git a/alg/teca_binary_segmentation.h b/alg/teca_binary_segmentation.h
index 983a497fb..138a28b03 100644
--- a/alg/teca_binary_segmentation.h
+++ b/alg/teca_binary_segmentation.h
@@ -10,17 +10,17 @@
TECA_SHARED_OBJECT_FORWARD_DECL(teca_binary_segmentation)
-/// an algorithm that computes a binary segmentation
+/// An algorithm that computes a binary segmentation.
/**
-an algorithm that computes a binary segmentation for 1D, 2D, and 3D data. The
-segmentation is computed using threshold operation where values in a range
-(low, high] are in the segmentation (assigned 1). Values outside the range
-are outside of the segmentation (assigned 0).
-
-The algorithm has 2 modes, BY_VALUE and BY_PERCENTILE. In the BY_VALUE mode,
-the test for inclusion is applied on the raw data. In the BY_PERCENTILE mode
-the range is given in percentiles and each data point is converted to a
-percentile before applying the test for inclusion.
+ * an algorithm that computes a binary segmentation for 1D, 2D, and 3D data.
+ * The segmentation is computed using threshold operation where values in a
+ * range (low, high] are in the segmentation (assigned 1). Values outside the
+ * range are outside of the segmentation (assigned 0).
+ *
+ * The algorithm has 2 modes, BY_VALUE and BY_PERCENTILE. In the BY_VALUE mode,
+ * the test for inclusion is applied on the raw data. In the BY_PERCENTILE mode
+ * the range is given in percentiles and each data point is converted to a
+ * percentile before applying the test for inclusion.
*/
class teca_binary_segmentation : public teca_algorithm
{
@@ -31,26 +31,59 @@ class teca_binary_segmentation : public teca_algorithm
~teca_binary_segmentation();
// set the name of the output array to store the resulting segmentation in
+ /** @name segmentation_variable
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(std::string, segmentation_variable)
+ ///@}
+
// set extra metadata for the segmentation variable
+ /** @name segmentation_variable_attributes
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(teca_metadata, segmentation_variable_attributes)
+ ///@}
+
- // set the name of the input array to segment
+ /** @name threshold_variable
+ * set the name of the input array to segment
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(std::string, threshold_variable)
+ ///@}
- // Set the threshold range. The defaults are (-infinity, infinity].
+ /** @name low_threshold_value
+ * Set the threshold range. The defaults are (-infinity, infinity].
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(double, low_threshold_value)
- TECA_ALGORITHM_PROPERTY(double, high_threshold_value)
+ ///@}
- // Set the threshold mode. In BY_PERCENTILE mode low and high thresholds
- // define the percentiles (0 to 100) between which data is in the
- // segmentation. default is BY_VALUE.
+ /** @name high_threshold_value
+ * Set the threshold range. The defaults are (-infinity, infinity].
+ */
+ ///@{
+ TECA_ALGORITHM_PROPERTY(double, high_threshold_value)
+ ///@}
+
+ /** @name threshold_mode
+ * Set the threshold mode. In BY_PERCENTILE mode low and high thresholds
+ * define the percentiles (0 to 100) between which data is in the
+ * segmentation. default is BY_VALUE.
+ */
+ ///@{
+ /// The threshold_mode modes
enum {BY_VALUE=0, BY_PERCENTILE=1};
- TECA_ALGORITHM_PROPERTY(int, threshold_mode);
+ TECA_ALGORITHM_PROPERTY(int, threshold_mode)
+
+ /// Set the threshold_mode to percentile.
void set_threshold_by_percentile() { set_threshold_mode(BY_PERCENTILE); }
+
+ /// Set the threshold_mode to value.
void set_threshold_by_value() { set_threshold_mode(BY_VALUE); }
+ ///@}
protected:
teca_binary_segmentation();
diff --git a/alg/teca_cartesian_mesh_regrid.cxx b/alg/teca_cartesian_mesh_regrid.cxx
index 632738efa..b8f90aa35 100644
--- a/alg/teca_cartesian_mesh_regrid.cxx
+++ b/alg/teca_cartesian_mesh_regrid.cxx
@@ -18,6 +18,28 @@ using std::endl;
//#define TECA_DEBUG
+// always use nearest neighbor interpolation for integers
+// to avoid truncation errors. an alternative would be to
+// implement rounding in the interpolator for integer types
+template
+int get_interpolation_mode(int desired_mode,
+ typename std::enable_if::value>::type* = 0)
+{
+ (void)desired_mode;
+ return teca_cartesian_mesh_regrid::nearest;
+}
+
+// use the requested interpolation mode for floating point
+// data
+template
+int get_interpolation_mode(int desired_mode,
+ typename std::enable_if::value>::type* = 0)
+{
+ return desired_mode;
+}
+
+
+// 3D
template
int interpolate(unsigned long target_nx, unsigned long target_ny,
unsigned long target_nz, const NT1 *p_target_xc, const NT1 *p_target_yc,
@@ -53,31 +75,100 @@ int interpolate(unsigned long target_nx, unsigned long target_ny,
return 0;
}
+// 2D - x-y
+template
+int interpolate(unsigned long target_nx, unsigned long target_ny,
+ const NT1 *p_target_xc, const NT1 *p_target_yc,
+ NT3 *p_target_a, const NT2 *p_source_xc,
+ const NT2 *p_source_yc, const NT3 *p_source_a,
+ unsigned long source_ihi, unsigned long source_jhi,
+ unsigned long source_nx)
+{
+ interp_t f;
+ unsigned long q = 0;
+ for (unsigned long j = 0; j < target_ny; ++j)
+ {
+ NT2 ty = static_cast(p_target_yc[j]);
+ for (unsigned long i = 0; i < target_nx; ++i, ++q)
+ {
+ NT2 tx = static_cast(p_target_xc[i]);
+ if (f(tx,ty,
+ p_source_xc, p_source_yc,
+ p_source_a, source_ihi, source_jhi,
+ source_nx, p_target_a[q]))
+ {
+ TECA_ERROR("failed to interpolate i=(" << i << ", " << j
+ << ") x=(" << tx << ", " << ty << ", " << ")")
+ return -1;
+ }
+ }
+ }
+ return 0;
+}
+
template
int interpolate(int mode, unsigned long target_nx, unsigned long target_ny,
- unsigned long target_nz, const taget_coord_t *p_target_xc, const taget_coord_t *p_target_yc,
- const taget_coord_t *p_target_zc, array_t *p_target_a, const source_coord_t *p_source_xc,
- const source_coord_t *p_source_yc, const source_coord_t *p_source_zc, const array_t *p_source_a,
- unsigned long source_ihi, unsigned long source_jhi, unsigned long source_khi,
- unsigned long source_nx, unsigned long source_nxy)
+ unsigned long target_nz, const taget_coord_t *p_target_xc,
+ const taget_coord_t *p_target_yc, const taget_coord_t *p_target_zc,
+ array_t *p_target_a, const source_coord_t *p_source_xc,
+ const source_coord_t *p_source_yc, const source_coord_t *p_source_zc,
+ const array_t *p_source_a, unsigned long source_ihi, unsigned long source_jhi,
+ unsigned long source_khi, unsigned long source_nx, unsigned long source_ny,
+ unsigned long source_nz)
{
using nearest_interp_t = teca_coordinate_util::interpolate_t<0>;
using linear_interp_t = teca_coordinate_util::interpolate_t<1>;
- switch (mode)
+ unsigned long source_nxy = source_nx*source_ny;
+
+ switch (get_interpolation_mode(mode))
{
case teca_cartesian_mesh_regrid::nearest:
- return interpolate(
- target_nx, target_ny, target_nz, p_target_xc, p_target_yc, p_target_zc,
- p_target_a, p_source_xc, p_source_yc, p_source_zc, p_source_a,
- source_ihi, source_jhi, source_khi, source_nx, source_nxy);
+ {
+ if ((target_nz == 1) && (source_nz == 1))
+ {
+ // 2D in the x-y plane
+ return interpolate(
+ target_nx, target_ny, p_target_xc, p_target_yc,
+ p_target_a, p_source_xc, p_source_yc, p_source_a,
+ source_ihi, source_jhi, source_nx);
+ }
+ else
+ {
+ // 3D
+ return interpolate(
+ target_nx, target_ny, target_nz, p_target_xc,
+ p_target_yc, p_target_zc, p_target_a, p_source_xc,
+ p_source_yc, p_source_zc, p_source_a, source_ihi,
+ source_jhi, source_khi, source_nx, source_nxy);
+ }
break;
+ }
case teca_cartesian_mesh_regrid::linear:
- return interpolate(
- target_nx, target_ny, target_nz, p_target_xc, p_target_yc, p_target_zc,
- p_target_a, p_source_xc, p_source_yc, p_source_zc, p_source_a,
- source_ihi, source_jhi, source_khi, source_nx, source_nxy);
+ {
+ if ((target_nz == 1) && (source_nz == 1))
+ {
+ // 2D in the x-y plane
+ return interpolate(
+ target_nx, target_ny, p_target_xc, p_target_yc,
+ p_target_a, p_source_xc, p_source_yc, p_source_a,
+ source_ihi, source_jhi, source_nx);
+ }
+ else
+ {
+ // 3D
+ return interpolate(
+ target_nx, target_ny, target_nz, p_target_xc,
+ p_target_yc, p_target_zc, p_target_a, p_source_xc,
+ p_source_yc, p_source_zc, p_source_a, source_ihi,
+ source_jhi, source_khi, source_nx, source_nxy);
+ }
break;
+ }
}
TECA_ERROR("invalid interpolation mode \"" << mode << "\"")
@@ -86,6 +177,7 @@ int interpolate(int mode, unsigned long target_nx, unsigned long target_ny,
+
// --------------------------------------------------------------------------
teca_cartesian_mesh_regrid::teca_cartesian_mesh_regrid()
: target_input(0), interpolation_mode(nearest)
@@ -109,13 +201,15 @@ void teca_cartesian_mesh_regrid::get_properties_description(
opts.add_options()
TECA_POPTS_GET(int, prefix, target_input,
- "select input connection that contains metadata (0)")
- TECA_POPTS_GET(std::vector, prefix, arrays,
- "list of arrays to move from source to target mesh ("")")
+ "select input connection that contains metadata")
+ TECA_POPTS_MULTI_GET(std::vector, prefix, arrays,
+ "list of arrays to move from source to target mesh")
TECA_POPTS_GET(int, prefix, interpolation_mode,
- "linear or nearest interpolation (1)")
+ "linear or nearest interpolation")
;
+ this->teca_algorithm::get_properties_description(prefix, opts);
+
global_opts.add(opts);
}
@@ -123,6 +217,8 @@ void teca_cartesian_mesh_regrid::get_properties_description(
void teca_cartesian_mesh_regrid::set_properties(
const std::string &prefix, variables_map &opts)
{
+ this->teca_algorithm::set_properties(prefix, opts);
+
TECA_POPTS_SET(opts, int, prefix, target_input)
TECA_POPTS_SET(opts, std::vector, prefix, arrays)
TECA_POPTS_SET(opts, int, prefix, interpolation_mode)
@@ -276,7 +372,9 @@ std::vector teca_cartesian_mesh_regrid::get_upstream_request(
else
{
if (teca_coordinate_util::bounds_to_extent(request_bounds,
- target_x, target_y, target_z, target_extent))
+ target_x, target_y, target_z, target_extent) ||
+ teca_coordinate_util::validate_extent(target_x->size(),
+ target_y->size(), target_z->size(), target_extent, true))
{
TECA_ERROR("invalid bounds requested [" << request_bounds[0] << ", "
<< request_bounds[1] << ", " << request_bounds[2] << ", "
@@ -294,6 +392,24 @@ std::vector teca_cartesian_mesh_regrid::get_upstream_request(
target_z->get(target_extent[4], target_bounds[4]);
target_z->get(target_extent[5], target_bounds[5]);
+ // if the source is 2D, the cf_reader may have faked the vertical dimension.
+ // in that case, use the source's vertical coordinate in the requested bounds
+ teca_metadata source_coords;
+ p_teca_variant_array source_z;
+
+ if (input_md[md_src].get("coordinates", source_coords)
+ || !(source_z = source_coords.get("z")))
+ {
+ TECA_ERROR("failed to locate source mesh coordinates")
+ return up_reqs;
+ }
+
+ if (source_z->size() == 1)
+ {
+ source_z->get(0, target_bounds[4]);
+ source_z->get(0, target_bounds[5]);
+ }
+
// send the target bounds to the source as well
source_req.set("bounds", target_bounds, 6);
@@ -301,10 +417,19 @@ std::vector teca_cartesian_mesh_regrid::get_upstream_request(
up_reqs[md_tgt] = target_req;
up_reqs[md_src] = source_req;
+#ifdef TECA_DEBUG
+ std::cerr << "source request = ";
+ source_req.to_stream(std::cerr);
+ std::cerr << std::endl;
+
+ std::cerr << "target request = ";
+ target_req.to_stream(std::cerr);
+ std::cerr << std::endl;
+#endif
+
return up_reqs;
}
-
// --------------------------------------------------------------------------
const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
unsigned int port, const std::vector &input_data,
@@ -347,7 +472,7 @@ const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
// get the list of arrays to move
std::vector req_arrays;
- request.get("regrid_arrays", req_arrays);
+ request.get("arrays", req_arrays);
// add any explicitly named
std::copy(this->arrays.begin(), this->arrays.end(),
@@ -362,7 +487,25 @@ const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
for (; it != end; ++it)
{
if (!target->get_point_arrays()->has(*it))
- source_arrays.push_back(*it);
+ {
+ if (source->get_point_arrays()->has(*it))
+ {
+ source_arrays.push_back(*it);
+ }
+ else
+ {
+ TECA_ERROR("Array \"" << *it
+ << "\" is neither present in source or target mesh")
+ return nullptr;
+ }
+ }
+ }
+
+ // catch a user error
+ if (!source_arrays.size() &&
+ teca_mpi_util::mpi_rank_0(this->get_communicator()))
+ {
+ TECA_WARNING("No arrays will be interpolated")
}
// move the arrays
@@ -384,7 +527,6 @@ const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
unsigned long source_nx = source_xc->size();
unsigned long source_ny = source_yc->size();
unsigned long source_nz = source_zc->size();
- unsigned long source_nxy = source_nx*source_ny;
unsigned long source_ihi = source_nx - 1;
unsigned long source_jhi = source_ny - 1;
unsigned long source_khi = source_nz - 1;
@@ -392,20 +534,20 @@ const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
NESTED_TEMPLATE_DISPATCH_FP(
const teca_variant_array_impl,
target_xc.get(),
- 1,
+ _TGT,
- const NT1 *p_target_xc = std::dynamic_pointer_cast(target_xc)->get();
- const NT1 *p_target_yc = std::dynamic_pointer_cast(target_yc)->get();
- const NT1 *p_target_zc = std::dynamic_pointer_cast(target_zc)->get();
+ const NT_TGT *p_target_xc = std::dynamic_pointer_cast(target_xc)->get();
+ const NT_TGT *p_target_yc = std::dynamic_pointer_cast(target_yc)->get();
+ const NT_TGT *p_target_zc = std::dynamic_pointer_cast(target_zc)->get();
NESTED_TEMPLATE_DISPATCH_FP(
const teca_variant_array_impl,
source_xc.get(),
- 2,
+ _SRC,
- const NT2 *p_source_xc = std::dynamic_pointer_cast(source_xc)->get();
- const NT2 *p_source_yc = std::dynamic_pointer_cast(source_yc)->get();
- const NT2 *p_source_zc = std::dynamic_pointer_cast(source_zc)->get();
+ const NT_SRC *p_source_xc = std::dynamic_pointer_cast(source_xc)->get();
+ const NT_SRC *p_source_yc = std::dynamic_pointer_cast(source_yc)->get();
+ const NT_SRC *p_source_zc = std::dynamic_pointer_cast(source_zc)->get();
size_t n_arrays = source_arrays.size();
for (size_t i = 0; i < n_arrays; ++i)
@@ -417,25 +559,37 @@ const_p_teca_dataset teca_cartesian_mesh_regrid::execute(
NESTED_TEMPLATE_DISPATCH(
teca_variant_array_impl,
target_a.get(),
- 3,
+ _DATA,
- const NT3 *p_source_a = std::static_pointer_cast(source_a)->get();
- NT3 *p_target_a = std::static_pointer_cast(target_a)->get();
+ const NT_DATA *p_source_a = std::static_pointer_cast(source_a)->get();
+ NT_DATA *p_target_a = std::static_pointer_cast(target_a)->get();
if (interpolate(this->interpolation_mode, target_nx, target_ny, target_nz,
p_target_xc, p_target_yc, p_target_zc, p_target_a, p_source_xc,
p_source_yc, p_source_zc, p_source_a, source_ihi, source_jhi,
- source_khi, source_nx, source_nxy))
+ source_khi, source_nx, source_ny, source_nz))
{
TECA_ERROR("Failed to move \"" << source_arrays[i] << "\"")
return nullptr;
}
-
- target_ac->set(source_arrays[i], target_a);
)
+ else
+ {
+ TECA_ERROR("Unsupported array type " << source_a->get_class_name())
}
- )
+
+ target_ac->set(source_arrays[i], target_a);
+ }
)
+ else
+ {
+ TECA_ERROR("Unupported coordinate type " << source_xc->get_class_name())
+ }
+ )
+ else
+ {
+ TECA_ERROR("Unupported coordinate type " << target_xc->get_class_name())
+ }
return target;
}
diff --git a/alg/teca_cartesian_mesh_regrid.h b/alg/teca_cartesian_mesh_regrid.h
index f9ae85da2..d67c5f7f0 100644
--- a/alg/teca_cartesian_mesh_regrid.h
+++ b/alg/teca_cartesian_mesh_regrid.h
@@ -4,26 +4,29 @@
#include "teca_shared_object.h"
#include "teca_algorithm.h"
#include "teca_metadata.h"
-#include "teca_variant_array_fwd.h"
+#include "teca_variant_array.h"
#include
#include
TECA_SHARED_OBJECT_FORWARD_DECL(teca_cartesian_mesh_regrid)
-/// transfer data between overlapping meshes of potentially different resolution
-/**
-an algorithm that transfers data between cartesian meshes defined in the same
-world coordinate system but potentially different resolutions. nearest or
-linear interpolation are supported.
-
-By default the first input is the target mesh. the second input is the source
-mesh. This can be changed by setting the target_input property.
-
-the arrays to move from source to target can be selected using add_array api or
-in the request key regrid_source_arrays. this is a spatial regriding operation
-for temporal regriding see teca_mesh_temporal_regrid.
-*/
+/** @brief
+ * Transfers data between spatially overlapping meshes of potentially different
+ * resolutions.
+ *
+ * @details
+ * an algorithm that transfers data between cartesian meshes defined in the
+ * same world coordinate system but potentially different resolutions. nearest
+ * or linear interpolation are supported.
+ *
+ * By default the first input is the target mesh. the second input is the
+ * source mesh. This can be changed by setting the target_input property.
+ *
+ * the arrays to move from source to target can be selected using add_array api
+ * or in the request key "arrays". this is a spatial regriding operation for
+ * temporal regriding see teca_mesh_temporal_regrid.
+ */
class teca_cartesian_mesh_regrid : public teca_algorithm
{
public:
@@ -37,23 +40,33 @@ class teca_cartesian_mesh_regrid : public teca_algorithm
TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION()
TECA_SET_ALGORITHM_PROPERTIES()
- // set the list of arrays to move from the source
- // to the target
+ /** @name array
+ * set the list of arrays to move from the source to the target
+ */
+ //@{
TECA_ALGORITHM_VECTOR_PROPERTY(std::string, array)
+ //@}
- // set the input connection from which metadata such as arrays
- // and time steps are taken from.
+ /** @name target_input
+ * set the input connection which provides the output geometry.
+ */
+ ///@{
TECA_ALGORITHM_PROPERTY(int, target_input)
-
- // set the interpolation mode used in transfering
- // data between meshes of differing resolution.
- // in nearest mode value at the nearest grid point
- // is used, in linear mode bi/tri linear interpolation
- // is used.
+ ///@}
+
+ /** @name interpolation_mode
+ * set the interpolation mode used in transfering data between meshes of
+ * differing resolution. in nearest mode value at the nearest grid point
+ * is used, in linear mode bi/tri linear interpolation is used.
+ */
+ //@{
enum {nearest=0, linear=1};
TECA_ALGORITHM_PROPERTY(int, interpolation_mode)
void set_interpolation_mode_nearest(){ interpolation_mode = nearest; }
void set_interpolation_mode_linear(){ interpolation_mode = linear; }
+ //@}
+
+
protected:
teca_cartesian_mesh_regrid();
diff --git a/alg/teca_cartesian_mesh_source.cxx b/alg/teca_cartesian_mesh_source.cxx
index 7946d188f..2139a9864 100644
--- a/alg/teca_cartesian_mesh_source.cxx
+++ b/alg/teca_cartesian_mesh_source.cxx
@@ -26,12 +26,19 @@ struct teca_cartesian_mesh_source::internals_t
// the world space [x0 x1 y0 y1 z0 z1 t0 t1] generate
// equally spaced coordinate axes x,y,z,t
static
- void initialize_axes(int type_code, unsigned long *extent,
- double *bounds, p_teca_variant_array &x_axis,
+ void initialize_axes(int type_code, const unsigned long *extent,
+ const double *bounds, p_teca_variant_array &x_axis,
p_teca_variant_array &y_axis, p_teca_variant_array &z_axis,
p_teca_variant_array &t_axis);
+ static
+ void initialize_axes(int type_code, const unsigned long *extent,
+ const double *bounds, p_teca_variant_array &x_axis,
+ p_teca_variant_array &y_axis, p_teca_variant_array &z_axis);
+
+ // cached metadata
teca_metadata metadata;
+ p_teca_variant_array t_axis;
};
@@ -45,21 +52,25 @@ void teca_cartesian_mesh_source::internals_t::initialize_axis(
unsigned long nx = i1 - i0 + 1;
x->resize(nx);
+ num_t *px = x->get();
+
// avoid divide by zero
if (nx < 2)
+ {
+ px[0] = x0;
return;
+ }
num_t dx = (x1 - x0)/(nx - 1l);
num_t xx = x0 + i0*dx;
- num_t *px = x->get();
for (unsigned long i = 0; i < nx; ++i)
px[i] = xx + dx*i;
}
// --------------------------------------------------------------------------
void teca_cartesian_mesh_source::internals_t::initialize_axes(int type_code,
- unsigned long *extent, double *bounds, p_teca_variant_array &x_axis,
+ const unsigned long *extent, const double *bounds, p_teca_variant_array &x_axis,
p_teca_variant_array &y_axis, p_teca_variant_array &z_axis,
p_teca_variant_array &t_axis)
{
@@ -86,6 +97,30 @@ void teca_cartesian_mesh_source::internals_t::initialize_axes(int type_code,
)
}
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::internals_t::initialize_axes(int type_code,
+ const unsigned long *extent, const double *bounds, p_teca_variant_array &x_axis,
+ p_teca_variant_array &y_axis, p_teca_variant_array &z_axis)
+{
+ // gernate equally spaced coordinate axes x,y,z,t
+ x_axis = teca_variant_array_factory::New(type_code);
+ y_axis = x_axis->new_instance();
+ z_axis = x_axis->new_instance();
+
+ TEMPLATE_DISPATCH(teca_variant_array_impl,
+ x_axis.get(),
+
+ internals_t::initialize_axis(std::static_pointer_cast(x_axis),
+ extent[0], extent[1], bounds[0], bounds[1]);
+
+ internals_t::initialize_axis(std::static_pointer_cast(y_axis),
+ extent[2], extent[3], bounds[2], bounds[3]);
+
+ internals_t::initialize_axis(std::static_pointer_cast(z_axis),
+ extent[4], extent[5], bounds[4], bounds[5]);
+ )
+}
+
// --------------------------------------------------------------------------
@@ -93,9 +128,7 @@ teca_cartesian_mesh_source::teca_cartesian_mesh_source() :
coordinate_type_code(teca_variant_array_code::get()),
field_type_code(teca_variant_array_code::get()),
x_axis_variable("lon"), y_axis_variable("lat"), z_axis_variable("plev"),
- t_axis_variable("time"), x_axis_units("degrees_east"),
- y_axis_units("degrees_north"), z_axis_units("pascals"),
- calendar("Gregorian"), time_units("seconds since 1970-01-01 00:00:00"),
+ t_axis_variable("time"),
whole_extents{0l, 359l, 0l, 179l, 0l, 0l, 0l, 0l},
bounds{0., 360, -90., 90., 0., 0., 0., 0.},
internals(new internals_t)
@@ -129,7 +162,6 @@ void teca_cartesian_mesh_source::set_properties(const std::string &prefix,
}
#endif
-
// --------------------------------------------------------------------------
void teca_cartesian_mesh_source::set_modified()
{
@@ -143,11 +175,311 @@ void teca_cartesian_mesh_source::set_modified()
void teca_cartesian_mesh_source::clear_cached_metadata()
{
this->internals->metadata.clear();
+ teca_algorithm::set_modified();
+}
+
+
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_x_axis_variable(const std::string &name)
+{
+ this->x_axis_variable = name;
+ this->x_axis_attributes.clear();
+ teca_algorithm::set_modified();
+}
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_x_axis_variable(const std::string &name,
+ const teca_metadata &atts)
+{
+ this->x_axis_variable = name;
+ this->x_axis_attributes = atts;
+ teca_algorithm::set_modified();
+}
+
+// --------------------------------------------------------------------------
+int teca_cartesian_mesh_source::set_x_axis_variable(const teca_metadata &md)
+{
+ teca_metadata coords;
+ if (md.get("coordinates", coords))
+ return -1;
+
+ if (coords.get("x_variable", this->x_axis_variable))
+ return -1;
+
+ teca_metadata atts;
+ if (md.get("attributes", atts))
+ return -1;
+
+ if (atts.get(this->x_axis_variable, this->x_axis_attributes))
+ return -1;
+
+ return 0;
+}
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_y_axis_variable(const std::string &name)
+{
+ this->y_axis_variable = name;
+ this->y_axis_attributes.clear();
+ teca_algorithm::set_modified();
+}
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_y_axis_variable(const std::string &name,
+ const teca_metadata &atts)
+{
+ this->y_axis_variable = name;
+ this->y_axis_attributes = atts;
+ teca_algorithm::set_modified();
+}
+
+// --------------------------------------------------------------------------
+int teca_cartesian_mesh_source::set_y_axis_variable(const teca_metadata &md)
+{
+ teca_metadata coords;
+ if (md.get("coordinates", coords))
+ return -1;
+
+ if (coords.get("y_variable", this->y_axis_variable))
+ return -1;
+
+ teca_metadata atts;
+ if (md.get("attributes", atts))
+ return -1;
+
+ if (atts.get(this->y_axis_variable, this->y_axis_attributes))
+ return -1;
+
+ return 0;
+}
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_z_axis_variable(const std::string &name)
+{
+ this->z_axis_variable = name;
+ this->z_axis_attributes.clear();
+ teca_algorithm::set_modified();
+}
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_z_axis_variable(const std::string &name,
+ const teca_metadata &atts)
+{
+ this->z_axis_variable = name;
+ this->z_axis_attributes = atts;
+ teca_algorithm::set_modified();
+}
+
+// --------------------------------------------------------------------------
+int teca_cartesian_mesh_source::set_z_axis_variable(const teca_metadata &md)
+{
+ // get coordinates and attributes, fail if either are missing
+ teca_metadata coords;
+ if (md.get("coordinates", coords))
+ return -1;
+
+ if (coords.get("x_variable", this->z_axis_variable))
+ return -1;
+
+ teca_metadata atts;
+ if (md.get("attributes", atts))
+ return -1;
+
+ if (atts.get(this->z_axis_variable, this->z_axis_attributes))
+ return -1;
+
+ return 0;
+}
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_t_axis_variable(const std::string &name)
+{
+ this->t_axis_variable = name;
+ this->t_axis_attributes.clear();
+ teca_algorithm::set_modified();
+}
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_calendar(
+ const std::string &calendar, const std::string &units)
+{
+ this->t_axis_attributes.clear();
+ this->t_axis_attributes.set("calendar", calendar);
+ this->t_axis_attributes.set("units", units);
+ teca_algorithm::set_modified();
+}
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_t_axis_variable(const std::string &name,
+ const teca_metadata &atts)
+{
+ this->t_axis_variable = name;
+ this->t_axis_attributes = atts;
+ teca_algorithm::set_modified();
+}
+
+// --------------------------------------------------------------------------
+int teca_cartesian_mesh_source::set_t_axis_variable(const teca_metadata &md)
+{
+ teca_metadata coords;
+ if (md.get("coordinates", coords))
+ return -1;
+
+ if (coords.get("t_variable", this->t_axis_variable))
+ return -1;
+
+ teca_metadata atts;
+ if (md.get("attributes", atts))
+ return -1;
+
+ if (atts.get(this->t_axis_variable, this->t_axis_attributes))
+ return -1;
+
+ return 0;
+}
+
+// --------------------------------------------------------------------------
+int teca_cartesian_mesh_source::set_t_axis(const teca_metadata &md)
+{
+ teca_metadata coords;
+ if (md.get("coordinates", coords))
+ return -1;
+
+ this->internals->t_axis = coords.get("t");
+
+ return 0;
+}
+
+// --------------------------------------------------------------------------
+void teca_cartesian_mesh_source::set_t_axis(const p_teca_variant_array &t)
+{
+ this->internals->t_axis = t;
+}
+
+// --------------------------------------------------------------------------
+int teca_cartesian_mesh_source::set_output_metadata(const teca_metadata &md)
+{
+ teca_metadata coords;
+ if (md.get("coordinates", coords))
+ return -1;
+
+ teca_metadata atts;
+ if (md.get("attributes", atts))
+ return -1;
+
+ // get the coordinate axes.
+ const_p_teca_variant_array x = coords.get("x");
+ const_p_teca_variant_array y = coords.get("y");
+ const_p_teca_variant_array z = coords.get("z");
+ const_p_teca_variant_array t = coords.get("t");
+
+ // because of assumptions made in execute, all must be provided
+ if (!x || !y || !z || !t)
+ return -1;
+
+ unsigned long nx = x->size();
+ unsigned long ny = y->size();
+ unsigned long nz = z->size();
+ unsigned long nxyz = nx*ny*nz;
+
+ // clear out any variables, and replace with those that we provide.
+ std::vector vars;
+ std::vector::iterator it = this->field_generators.begin();
+ std::vector::iterator end = this->field_generators.end();
+ for (; it != end; ++it)
+ {
+ vars.push_back(it->name);
+
+ // correct size
+ teca_metadata var_atts = it->attributes;
+ var_atts.set("size", nxyz);
+
+ atts.set(it->name, var_atts);
+ }
+
+ // copy the metadata
+ this->set_modified();
+
+ this->internals->metadata = md;
+ this->internals->metadata.set("variables", vars);
+
+ return 0;
+}
+
+// --------------------------------------------------------------------------
+int teca_cartesian_mesh_source::set_spatial_extents(const teca_metadata &md,
+ bool three_d)
+{
+ // get coordinates and attributes, fail if either are missing
+ teca_metadata coords;
+ if (md.get("coordinates", coords))
+ return -1;
+
+ teca_metadata attributes;
+ if (md.get("attributes", attributes))
+ return -1;
+
+ // get the coordinate axes
+ p_teca_variant_array x = coords.get("x");
+ p_teca_variant_array y = coords.get("y");
+ p_teca_variant_array z = coords.get("z");
+
+ // verify
+ if (!x || !y || (three_d && !z))
+ return -1;
+
+ // set the extents
+ this->whole_extents[0] = 0;
+ this->whole_extents[1] = x->size() - 1;
+ this->whole_extents[2] = 0;
+ this->whole_extents[3] = y->size() - 1;
+ this->whole_extents[4] = 0;
+ this->whole_extents[5] = three_d ? z->size() - 1 : 0;
+
+ return 0;
+}
+// --------------------------------------------------------------------------
+int teca_cartesian_mesh_source::set_spatial_bounds(const teca_metadata &md,
+ bool three_d)
+{
+ // get coordinates and attributes, fail if either are missing
+ teca_metadata coords;
+ if (md.get("coordinates", coords))
+ return -1;
+
+ teca_metadata attributes;
+ if (md.get("attributes", attributes))
+ return -1;
+
+ // get the coordinate axes
+ p_teca_variant_array x = coords.get("x");
+ p_teca_variant_array y = coords.get("y");
+ p_teca_variant_array z = coords.get("z");
+
+ // verify
+ if (!x || !y || (three_d && !z))
+ return -1;
+
+ // get the bounds
+ x->get(0lu, this->bounds[0]);
+ x->get(x->size() - 1lu, this->bounds[1]);
+ y->get(0lu, this->bounds[2]);
+ y->get(y->size() - 1lu, this->bounds[3]);
+ z->get(0lu, this->bounds[4]);
+
+ unsigned long khi = three_d ? z->size() - 1lu : 0lu;
+ z->get(khi, this->bounds[5]);
+
+ // set the coordinate type
+ this->set_coordinate_type_code(x->type_code());
+
+ return 0;
}
// --------------------------------------------------------------------------
void teca_cartesian_mesh_source::append_field_generator(
- const std::string &name, const teca_array_attributes &atts,
+ const std::string &name, const teca_metadata &atts,
field_generator_callback &callback)
{
this->append_field_generator({name, atts, callback});
@@ -181,59 +513,59 @@ teca_metadata teca_cartesian_mesh_source::get_output_metadata(
// generate cooridnate axes
p_teca_variant_array x_axis, y_axis, z_axis, t_axis;
- internals_t::initialize_axes(this->coordinate_type_code,
- this->whole_extents.data(), this->bounds.data(), x_axis,
- y_axis, z_axis, t_axis);
+ if (this->internals->t_axis)
+ {
+ // generate x,y,z axes but use cached time axis
+ internals_t::initialize_axes(this->coordinate_type_code,
+ this->whole_extents.data(), this->bounds.data(), x_axis,
+ y_axis, z_axis);
- size_t nx = this->whole_extents[1] - this->whole_extents[0] + 1;
- size_t ny = this->whole_extents[3] - this->whole_extents[2] + 1;
- size_t nz = this->whole_extents[5] - this->whole_extents[4] + 1;
- size_t nt = this->whole_extents[7] - this->whole_extents[6] + 1;
- size_t nxyz = nx*ny*nz;
+ t_axis = this->internals->t_axis;
+ }
+ else
+ {
+ // generate x,y,z and t axes
+ internals_t::initialize_axes(this->coordinate_type_code,
+ this->whole_extents.data(), this->bounds.data(), x_axis,
+ y_axis, z_axis, t_axis);
+ }
- std::string x_ax_var_name = (this->x_axis_variable.empty() ? "x" : this->x_axis_variable);
- std::string y_ax_var_name = (this->y_axis_variable.empty() ? "y" : this->y_axis_variable);
- std::string z_ax_var_name = (this->z_axis_variable.empty() ? "z" : this->z_axis_variable);
- std::string t_ax_var_name = (this->t_axis_variable.empty() ? "t" : this->t_axis_variable);
+ size_t nx = x_axis->size();
+ size_t ny = y_axis->size();
+ size_t nz = z_axis->size();
+ size_t nt = t_axis->size();
// construct attributes
- teca_metadata x_atts;
- x_atts.set("units", (this->x_axis_units.empty() ? "meters" : this->x_axis_units));
- x_atts.set("type_code", this->coordinate_type_code);
+ teca_metadata x_atts = this->x_axis_attributes;
+ x_atts.set("type_code", x_axis->type_code());
x_atts.set("size", nx);
- teca_metadata y_atts;
- y_atts.set("units", (this->y_axis_units.empty() ? "meters" : this->y_axis_units));
- y_atts.set("type_code", this->coordinate_type_code);
+ teca_metadata y_atts = this->y_axis_attributes;
+ y_atts.set("type_code", y_axis->type_code());
y_atts.set("size", ny);
- teca_metadata z_atts;
- z_atts.set("units", (this->z_axis_units.empty() ? "meters" : this->z_axis_units));
- z_atts.set("type_code", this->coordinate_type_code);
+ teca_metadata z_atts = this->z_axis_attributes;
+ z_atts.set("type_code", z_axis->type_code());
z_atts.set("size", nz);
- teca_metadata t_atts;
- t_atts.set("units", (this->time_units.empty() ?
- "seconds since 1970-01-01 00:00:00" : this->time_units));
-
- t_atts.set("calendar", (this->calendar.empty() ?
- "standard" : this->calendar));
-
- t_atts.set("type_code", this->coordinate_type_code);
+ teca_metadata t_atts = this->t_axis_attributes;
+ t_atts.set("type_code", t_axis->type_code());
t_atts.set("size", nt);
teca_metadata atts;
- atts.set(x_ax_var_name, x_atts);
- atts.set(y_ax_var_name, y_atts);
- atts.set(z_ax_var_name, z_atts);
- atts.set(t_ax_var_name, t_atts);
+ atts.set(this->x_axis_variable, x_atts);
+ atts.set(this->y_axis_variable, y_atts);
+ atts.set(this->z_axis_variable, z_atts);
+
+ if (!this->t_axis_variable.empty())
+ atts.set(this->t_axis_variable, t_atts);
// construct dataset metadata
teca_metadata coords;
- coords.set("x_variable", x_ax_var_name);
- coords.set("y_variable", y_ax_var_name);
- coords.set("z_variable", z_ax_var_name);
- coords.set("t_variable", t_ax_var_name);
+ coords.set("x_variable", this->x_axis_variable);
+ coords.set("y_variable", this->y_axis_variable);
+ coords.set("z_variable", this->z_axis_variable);
+ coords.set("t_variable", this->t_axis_variable);
coords.set("x", x_axis);
coords.set("y", y_axis);
@@ -243,6 +575,7 @@ teca_metadata teca_cartesian_mesh_source::get_output_metadata(
this->internals->metadata.set("whole_extent", this->whole_extents);
this->internals->metadata.set("coordinates", coords);
+ size_t nxyz = nx*ny*nz;
std::vector vars;
std::vector::iterator it = this->field_generators.begin();
std::vector::iterator end = this->field_generators.end();
@@ -251,18 +584,24 @@ teca_metadata teca_cartesian_mesh_source::get_output_metadata(
vars.push_back(it->name);
// correct size
- teca_array_attributes var_atts = it->attributes;
- var_atts.size = nxyz;
+ teca_metadata var_atts = it->attributes;
+ var_atts.set("size", nxyz);
- atts.set(it->name, teca_metadata(var_atts));
+ atts.set(it->name, var_atts);
}
this->internals->metadata.set("variables", vars);
this->internals->metadata.set("attributes", atts);
- this->internals->metadata.set("number_of_time_steps", t_axis->size());
- this->internals->metadata.set("index_initializer_key", std::string("number_of_time_steps"));
- this->internals->metadata.set("index_request_key", std::string("time_step"));
+ // setup the execution control keys
+ this->internals->metadata.set("number_of_time_steps",
+ t_axis->size());
+
+ this->internals->metadata.set("index_initializer_key",
+ std::string("number_of_time_steps"));
+
+ this->internals->metadata.set("index_request_key",
+ std::string("time_step"));
return this->internals->metadata;
}
@@ -320,25 +659,44 @@ const_p_teca_dataset teca_cartesian_mesh_source::execute(unsigned int port,
{
// bounds key was present, convert the bounds to an
// an extent that covers them.
- if (teca_coordinate_util::bounds_to_extent(
- req_bounds, in_x, in_y, in_z, req_extent))
+ if (teca_coordinate_util::bounds_to_extent(req_bounds,
+ in_x, in_y, in_z, req_extent) ||
+ teca_coordinate_util::validate_extent(in_x->size(),
+ in_y->size(), in_z->size(), req_extent, true))
{
TECA_ERROR("invalid bounds requested.")
return nullptr;
}
}
- // get the timestep
- unsigned long time_step = 0;
- if (request.get("time_step", time_step))
+ // get the timestep, no matter what the key is named we treat it as
+ // a time step. this is to support metadata provided by another source
+ // eg. a different reader.
+ std::string request_key;
+ if (request.get("index_request_key", request_key))
+ {
+ TECA_ERROR("Request is missing the \"index_request_key\"")
+ return nullptr;
+ }
+
+ unsigned long req_index = 0;
+ if (request.get(request_key, req_index))
+ {
+ TECA_ERROR("Request is missing \"" << request_key << "\"")
+ return nullptr;
+ }
+
+ // check that the we have a time value for the requested index.
+ if (req_index >= in_t->size())
{
- TECA_ERROR("Request is missing time_step")
+ TECA_ERROR("The requested index " << req_index
+ << " is out of bounds [0, " << in_t->size() << "]")
return nullptr;
}
// get the time
double t = 0.;
- in_t->get(time_step, t);
+ in_t->get(req_index, t);
// slice axes on the requested extent
p_teca_variant_array out_x = in_x->new_copy(req_extent[0], req_extent[1]);
@@ -351,22 +709,29 @@ const_p_teca_dataset teca_cartesian_mesh_source::execute(unsigned int port,
std::string x_variable = this->x_axis_variable.empty() ? "x" : this->x_axis_variable;
std::string y_variable = this->y_axis_variable.empty() ? "y" : this->y_axis_variable;
std::string z_variable = this->z_axis_variable.empty() ? "z" : this->z_axis_variable;
- std::string t_variable = this->t_axis_variable.empty() ? "t" : this->t_axis_variable;
mesh->set_x_coordinates(x_variable, out_x);
mesh->set_y_coordinates(y_variable, out_y);
mesh->set_z_coordinates(z_variable, out_z);
+ // get the calendar
+ std::string calendar;
+ std::string units;
+ teca_metadata atts;
+ this->internals->metadata.get("attributes", atts);
+ atts.get("calendar", calendar);
+ atts.get("units", units);
+
// set metadata
mesh->set_whole_extent(md_whole_extent);
mesh->set_extent(req_extent);
- mesh->set_time_step(time_step);
+ mesh->set_time_step(req_index);
mesh->set_time(t);
- mesh->set_calendar(this->calendar);
- mesh->set_time_units(this->time_units);
+ mesh->set_calendar(calendar);
+ mesh->set_time_units(units);
teca_metadata &mesh_md = mesh->get_metadata();
- mesh_md.set("index_request_key", std::string("time_step"));
+ mesh_md.set("index_request_key", request_key);
// generate fields over the requested subset
std::vector::iterator it = this->field_generators.begin();
@@ -378,8 +743,6 @@ const_p_teca_dataset teca_cartesian_mesh_source::execute(unsigned int port,
}
// pass the attributes
- teca_metadata atts;
- this->internals->metadata.get("attributes", atts);
mesh_md.set("attributes", atts);
return mesh;
diff --git a/alg/teca_cartesian_mesh_source.h b/alg/teca_cartesian_mesh_source.h
index be1e5039a..d357ad0ba 100644
--- a/alg/teca_cartesian_mesh_source.h
+++ b/alg/teca_cartesian_mesh_source.h
@@ -2,7 +2,7 @@
#define teca_cartesian_mesh_source_h
#include "teca_algorithm.h"
-#include "teca_array_attributes.h"
+#include "teca_metadata.h"
#include
#include