diff --git a/.travis.yml b/.travis.yml index 1a686bd99..9193d6af8 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=112 + - TECA_DATA_REVISION=115 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/alg/CMakeLists.txt b/alg/CMakeLists.txt index cbcd774fe..7bb4bf52f 100644 --- a/alg/CMakeLists.txt +++ b/alg/CMakeLists.txt @@ -21,11 +21,13 @@ 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 diff --git a/alg/teca_bayesian_ar_detect.cxx b/alg/teca_bayesian_ar_detect.cxx index 48069c09a..c14dd6f1f 100644 --- a/alg/teca_bayesian_ar_detect.cxx +++ b/alg/teca_bayesian_ar_detect.cxx @@ -608,7 +608,9 @@ 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), + 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); @@ -631,18 +633,21 @@ 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") + "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") + "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") + "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") + "Set the number of threads to parallelize execution over.") ; this->teca_algorithm::get_properties_description(prefix, opts); @@ -660,6 +665,7 @@ void teca_bayesian_ar_detect::set_properties(const std::string &prefix, 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) } @@ -809,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; @@ -820,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(); @@ -879,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"); @@ -1043,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 a8cf0989b..b8fe86ed8 100644 --- a/alg/teca_bayesian_ar_detect.h +++ b/alg/teca_bayesian_ar_detect.h @@ -62,6 +62,14 @@ class teca_bayesian_ar_detect : public teca_algorithm TECA_ALGORITHM_PROPERTY(std::string, min_component_area_variable) TECA_ALGORITHM_PROPERTY(std::string, hwhm_latitude_variable) + /** @anchor probability variable + * @name probability variable + * Set the name of the variable to store output probability as. + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, ar_probability_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. @@ -99,6 +107,7 @@ 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; struct internals_t; diff --git a/alg/teca_cartesian_mesh_source.cxx b/alg/teca_cartesian_mesh_source.cxx index 3266b59c5..ee4349ec7 100644 --- a/alg/teca_cartesian_mesh_source.cxx +++ b/alg/teca_cartesian_mesh_source.cxx @@ -549,7 +549,7 @@ teca_metadata teca_cartesian_mesh_source::get_output_metadata( z_atts.set("size", nz); teca_metadata t_atts = this->t_axis_attributes; - t_atts.set("type_code", z_axis->type_code()); + t_atts.set("type_code", t_axis->type_code()); t_atts.set("size", nt); teca_metadata atts; @@ -707,7 +707,6 @@ 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); diff --git a/alg/teca_elevation_mask.cxx b/alg/teca_elevation_mask.cxx new file mode 100644 index 000000000..bc662a313 --- /dev/null +++ b/alg/teca_elevation_mask.cxx @@ -0,0 +1,384 @@ +#include "teca_elevation_mask.h" + +#include "teca_cartesian_mesh.h" +#include "teca_array_collection.h" +#include "teca_variant_array.h" +#include "teca_metadata.h" +#include "teca_array_attributes.h" + +#include "teca_dataset_source.h" +#include "teca_dataset_capture.h" +#include "teca_cartesian_mesh_regrid.h" +#include "teca_index_executive.h" + +#include +#include +#include +#include +#include + +#if defined(TECA_HAS_BOOST) +#include +#endif + +//#define TECA_DEBUG + +struct teca_elevation_mask::internals_t +{ + // compute the valid value mask such that for each point the mask + // is 1 where the mesh point is above the surface of the Earth and + // 0 otherwise + template + static void mask_by_surface_elevation( + size_t nx, size_t ny, size_t nz, + mask_t * __restrict__ mask, + const elev_num_t * __restrict__ surface_elev, + const mesh_num_t * __restrict__ mesh_height) + { + size_t nxy = nx*ny; + for (size_t k = 0; k < nz; ++k) + { + const mesh_num_t * __restrict__ mesh_height_k = mesh_height + k*nxy; + mask_t * __restrict__ mask_k = mask + k*nxy; + for (size_t q = 0; q < nxy; ++q) + { + mask_k[q] = mesh_height_k[q] >= (mesh_num_t)surface_elev[q] ? mask_t(1) : mask_t(0); + } + } + } +}; + + +// -------------------------------------------------------------------------- +teca_elevation_mask::teca_elevation_mask() : + mesh_height_variable("zg"), surface_elevation_variable("z") +{ + this->set_number_of_input_connections(2); + this->set_number_of_output_ports(1); +} + +// -------------------------------------------------------------------------- +teca_elevation_mask::~teca_elevation_mask() +{ +} + +#if defined(TECA_HAS_BOOST) +// -------------------------------------------------------------------------- +void teca_elevation_mask::get_properties_description( + const std::string &prefix, options_description &global_opts) +{ + options_description opts("Options for " + + (prefix.empty()?"teca_elevation_mask":prefix)); + + opts.add_options() + + TECA_POPTS_GET(std::string, prefix, surface_elevation_variable, + "Set the name of the variable containing surface elevation" + " values in meters above mean sea level") + + TECA_POPTS_GET(std::string, prefix, mesh_height_variable, + "Set the name of the variable containing point wise mesh height" + " values in meters above mean sea level") + + TECA_POPTS_MULTI_GET(std::vector, prefix, mask_variables, + "Set the names of the variables to store the generated mask in." + " Each name is assigned a reference to the mask.") + ; + + this->teca_algorithm::get_properties_description(prefix, opts); + + global_opts.add(opts); +} + +// -------------------------------------------------------------------------- +void teca_elevation_mask::set_properties( + const std::string &prefix, variables_map &opts) +{ + this->teca_algorithm::set_properties(prefix, opts); + + TECA_POPTS_SET(opts, std::string, prefix, surface_elevation_variable) + TECA_POPTS_SET(opts, std::string, prefix, mesh_height_variable) + TECA_POPTS_SET(opts, std::vector, prefix, mask_variables) +} +#endif + +// -------------------------------------------------------------------------- +teca_metadata teca_elevation_mask::get_output_metadata( + unsigned int port, + const std::vector &input_md) +{ +#ifdef TECA_DEBUG + cerr << teca_parallel_id() + << "teca_elevation_mask::get_output_metadata" << endl; +#endif + (void)port; + + // validate runtime provided settings + unsigned int n_mask_vars = this->mask_variables.size(); + if (n_mask_vars == 0) + { + TECA_ERROR("The names of the mask_variables were not provided") + return teca_metadata(); + } + + // pass metadata from the input mesh through. + const teca_metadata &mesh_md = input_md[0]; + teca_metadata out_md(mesh_md); + + // add the mask arrays we will generate + for (unsigned int i = 0; i < n_mask_vars; ++i) + out_md.append("variables", this->mask_variables[i]); + + // insert attributes to enable this to be written by the CF writer + teca_metadata attributes; + out_md.get("attributes", attributes); + + teca_metadata mesh_height_atts; + if (attributes.get(this->mesh_height_variable, mesh_height_atts)) + { + TECA_WARNING("Failed to get mesh_height_variable \"" + << this->mesh_height_variable << "\" attrbibutes." + " Writing the result will not be possible") + } + else + { + // get the centering and size from the array + unsigned int centering = 0; + mesh_height_atts.get("centering", centering); + + unsigned long size = 0; + mesh_height_atts.get("size", size); + + // construct output attributes + teca_array_attributes mask_atts( + teca_variant_array_code::get(), + centering, size, "none", "", "elevation mask"); + + // add one for each output + for (unsigned int i = 0; i < n_mask_vars; ++i) + attributes.set(this->mask_variables[i], (teca_metadata)mask_atts); + + // update the attributes collection + out_md.set("attributes", attributes); + } + + return out_md; +} + +// -------------------------------------------------------------------------- +std::vector teca_elevation_mask::get_upstream_request( + unsigned int port, + const std::vector &input_md, + const teca_metadata &request) +{ + (void)port; + (void)input_md; + + std::vector up_reqs; + + // get the names of the arrays we need to request + if (this->mesh_height_variable.empty()) + { + TECA_ERROR("The mesh_height_variable was not specified") + return up_reqs; + } + + if (this->surface_elevation_variable.empty()) + { + TECA_ERROR("The surface_elevation_variable was not specified") + return up_reqs; + } + + // need to make the request for the surface elevation field using bounds + double req_bounds[6] = {0.0}; + if (request.get("bounds", req_bounds, 6)) + { + // bounds not specified, try to get an extent and convert to a bounds + unsigned long req_extent[6]; + if (request.get("extent", req_extent, 6)) + { + TECA_ERROR("Neither bounds nor extent were specified in the request") + return up_reqs; + } + + const teca_metadata &md = input_md[0]; + + teca_metadata coords; + p_teca_variant_array x,y; + + if (md.get("coordinates", coords) || + !(x = coords.get("x")) || !(y = coords.get("y"))) + { + TECA_ERROR("Failed to get mesh coordinates") + return up_reqs; + } + + x->get(req_extent[0], req_bounds[0]); + x->get(req_extent[1], req_bounds[1]); + y->get(req_extent[2], req_bounds[2]); + y->get(req_extent[3], req_bounds[3]); + } + + // input port 0 will source the mesh height field, and any other data + // requested by the down stream. copy the incoming request to preserve the + // downstream requirements and add the mesh height variable + teca_metadata req_0(request); + + std::set mesh_arrays; + if (req_0.has("arrays")) + req_0.get("arrays", mesh_arrays); + + mesh_arrays.insert(this->mesh_height_variable); + + // intercept request for our output + int n_mask_vars = this->mask_variables.size(); + for (int i = 0; i < n_mask_vars; ++i) + mesh_arrays.erase(this->mask_variables[i]); + + req_0.set("arrays", mesh_arrays); + + + // input port 1 provides the surface elevation field, request it + // preserve bounds etc + const teca_metadata &elev_md = input_md[1]; + + std::string req_key; + if (elev_md.get("index_request_key", req_key)) + { + TECA_ERROR("Metadata is missing \"index_request_key\"") + return up_reqs; + } + + // surface elevations don't change over the timescale of concern + // always request index 0 + teca_metadata req_1; + req_1.set(req_key, 0ul); + req_1.set("index_request_key", req_key); + + // request the surface elevation + std::vector elev_arrays(1, this->surface_elevation_variable); + req_1.set("arrays", elev_arrays); + + // at the bounds of interest + req_1.set("bounds", req_bounds, 6); + + // package the requests and send them up + up_reqs.push_back(req_0); + up_reqs.push_back(req_1); + + + return up_reqs; +} + +// -------------------------------------------------------------------------- +const_p_teca_dataset teca_elevation_mask::execute( + unsigned int port, + const std::vector &input_data, + const teca_metadata &request) +{ +#ifdef TECA_DEBUG + cerr << teca_parallel_id() << "teca_elevation_mask::execute" << endl; +#endif + (void)port; + (void)request; + + // check for an error upstream + if ((input_data.size() != 2) || !input_data[0] || !input_data[1]) + { + TECA_ERROR("Invalid inputs detected") + return nullptr; + } + + // get the input 3D mesh + const_p_teca_cartesian_mesh in_mesh + = std::dynamic_pointer_cast(input_data[0]); + + if (!in_mesh) + { + TECA_ERROR("Data to mask on input port 0 is not a" + " teca_cartesian_mesh. Got " << input_data[0]->get_class_name()) + return nullptr; + } + + // get the mesh dimensions + unsigned long extent[6]; + in_mesh->get_extent(extent); + + unsigned long nx = extent[1] - extent[0] + 1; + unsigned long ny = extent[3] - extent[2] + 1; + unsigned long nz = extent[5] - extent[4] + 1; + + // get the mesh height, this is a 3d field with the altitude for + // each mesh point + const_p_teca_variant_array mesh_height = + in_mesh->get_point_arrays()->get(this->mesh_height_variable); + + if (!mesh_height) + { + TECA_ERROR("Mesh to mask is missing the height field \"" + << this->mesh_height_variable << "\"") + return nullptr; + } + + // get the surface elevations + const_p_teca_cartesian_mesh in_elev + = std::dynamic_pointer_cast(input_data[1]); + + if (!in_elev) + { + TECA_ERROR("Data to mask on input port 0 is not a" + " teca_cartesian_mesh. Got " << input_data[0]->get_class_name()) + return nullptr; + } + + // get the surface elevation, this is a 2d field with surface altitude + // + // at each mesh point. regridding has been performed so that the horizontal + // coordinates are the same as the 3d mesh for which masks will be generated + const_p_teca_variant_array surface_elev = + in_elev->get_point_arrays()->get(this->surface_elevation_variable); + + if (!surface_elev) + { + TECA_ERROR("Surface elevation data has no array \"" + << this->surface_elevation_variable << "\"") + return nullptr; + } + + // compute the mask + p_teca_char_array mask = teca_char_array::New(mesh_height->size()); + char *p_mask = mask->get(); + + NESTED_TEMPLATE_DISPATCH(const teca_variant_array_impl, + surface_elev.get(), + _SURF, + + const NT_SURF *p_surface_elev = + static_cast(surface_elev.get())->get(); + + NESTED_TEMPLATE_DISPATCH(const teca_variant_array_impl, + mesh_height.get(), + _MESH, + + const NT_MESH *p_mesh_height = + static_cast(mesh_height.get())->get(); + + internals_t::mask_by_surface_elevation(nx, ny, nz, + p_mask, p_surface_elev, p_mesh_height); + + ) + ) + + // allocate the output mesh + p_teca_cartesian_mesh out_mesh = std::dynamic_pointer_cast + (std::const_pointer_cast(in_mesh)->new_shallow_copy()); + + // store the results under the requested names + int n_mask_vars = this->mask_variables.size(); + for (int i = 0; i < n_mask_vars; ++i) + { + out_mesh->get_point_arrays()->set(this->mask_variables[i], mask); + } + + return out_mesh; +} diff --git a/alg/teca_elevation_mask.h b/alg/teca_elevation_mask.h new file mode 100644 index 000000000..dc9c4b4c6 --- /dev/null +++ b/alg/teca_elevation_mask.h @@ -0,0 +1,115 @@ +#ifndef teca_elevation_mask_h +#define teca_elevation_mask_h + +#include "teca_shared_object.h" +#include "teca_algorithm.h" +#include "teca_metadata.h" + +#include +#include + +TECA_SHARED_OBJECT_FORWARD_DECL(teca_elevation_mask) + +/** @brief + * Generates a mask indicating where mesh points with a vertical pressure + * coordinate lie above the surface of the Earth. The mask is set to 1 where + * data is above the Earth's surface and 0 otherwise. + * + * @details + * Given a 3D height field containing the altitude of each point in meters + * above mean sea level, and a 2D height field corresponding to height in + * meters above mean sea level of the surface of the Earth, generate a mask + * that is 1 where the 3D point is on or above the surface of the Earth and 0 + * where it is below. + * + * The name of the 3D height field is specified by the @ref mesh_height_variable + * property. The name of the 2D height field conntaining elveation of the + * Earth's surface is specified by the @ref surface_elevation_variable property. + * + * The 3D mesh height field must be provided on input 0, and the 2D surface + * height field on input 1. Use the mask_names property to name the output + * mask. If more than one name is provided each name will reference a pointer + * to the mask. Consider using names of the form X_valid in which case the + * output is compatible with the teca_valid_value_mask and will be treated + * as missing values by down stream algorithms. + * + * If the simulation does not provide the 3D height field, for simulations + * where the acceleration due to the Earth's gravity is assumed constant, + * teca_geopotential_height can generate the 3D height field. + * + * The primary use case of this algorithm is when dealing with calculations on + * 3D meshes with a vertical pressure coordinate and there is a need to + * identify and treat specially the mesh points that are below the surface of + * the Earth. There are a number of alternatives available depending on the + * data. If your data has a _FillValue where data is below the surface then + * use teca_valid_value_mask instead of this algorithm. If your data has + * surface pressure field use teca_pressure_level_mask instead of this + * algorithm. If your dataset has surface temperature, and mean sea level + * pressure fields then use teca_surface_pressure to generate the surface + * pressure field and use teca_pressure_level_mask instead of this algorithm. + */ +class teca_elevation_mask : public teca_algorithm +{ +public: + TECA_ALGORITHM_STATIC_NEW(teca_elevation_mask) + TECA_ALGORITHM_DELETE_COPY_ASSIGN(teca_elevation_mask) + TECA_ALGORITHM_CLASS_NAME(teca_elevation_mask) + ~teca_elevation_mask(); + + // report/initialize to/from Boost program options + // objects. + TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION() + TECA_SET_ALGORITHM_PROPERTIES() + + /** @anchor mesh_height_variable + * @name algorithm property: mesh_height_variable + * set/get the name of the 3D height field + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, mesh_height_variable) + ///@} + + /** @anchor surface_elevation_variable + * @name algorithm property: surface_elevation_variable + * set/get the name of the variable containing the elevation of the Earth's + * surface. + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, surface_elevation_variable) + ///@} + + /** @anchor mask_variables + * @name algorithm property: mask_variable + * set the names of the variables to store the generated mask in + * each variable will contain a reference to the mask + */ + ///@{ + TECA_ALGORITHM_VECTOR_PROPERTY(std::string, mask_variable) + ///@} + +protected: + teca_elevation_mask(); + +private: + teca_metadata get_output_metadata( + unsigned int port, + const std::vector &input_md) override; + + std::vector get_upstream_request( + unsigned int port, + const std::vector &input_md, + const teca_metadata &request) override; + + const_p_teca_dataset execute( + unsigned int port, + const std::vector &input_data, + const teca_metadata &request) override; + +private: + std::string mesh_height_variable; + std::string surface_elevation_variable; + std::vector mask_variables; + struct internals_t; +}; + +#endif diff --git a/alg/teca_integrated_vapor_transport.h b/alg/teca_integrated_vapor_transport.h index 44bd896a4..00655130e 100644 --- a/alg/teca_integrated_vapor_transport.h +++ b/alg/teca_integrated_vapor_transport.h @@ -10,19 +10,19 @@ TECA_SHARED_OBJECT_FORWARD_DECL(teca_integrated_vapor_transport) -/// an algorithm that computes integrated vapor transport (IVT) +/// An algorithm that computes integrated vapor transport (IVT) /** -Compute integrated vaport transport (IVT) from wind vector and -specific humidity. - -IVT = - \frac{1}{g} \int_{p_0}^{p_1} \vec{v} q dp - -where q is the specific humidity, and \vec{v} = (u, v) are the -longitudinal and latitudinal components of wind. - -This calculation is an instance of a vertical reduction where -a 3D mesh is transformed into a 2D one. -*/ + * Compute integrated vapor transport (IVT) from wind vector and + * specific humidity. + * + * IVT = \frac{1}{g} \int_{p_{sfc}}^{p_{top}} \vec{v} q dp + * + * where q is the specific humidity, and \vec{v} = (u, v) are the + * longitudinal and latitudinal components of wind. + * + * This calculation is an instance of a vertical reduction where + * a 3D mesh is transformed into a 2D one. + */ class teca_integrated_vapor_transport : public teca_vertical_reduction { public: diff --git a/alg/teca_integrated_water_vapor.cxx b/alg/teca_integrated_water_vapor.cxx new file mode 100644 index 000000000..481dfc690 --- /dev/null +++ b/alg/teca_integrated_water_vapor.cxx @@ -0,0 +1,362 @@ +#include "teca_integrated_water_vapor.h" + +#include "teca_cartesian_mesh.h" +#include "teca_array_collection.h" +#include "teca_variant_array.h" +#include "teca_metadata.h" +#include "teca_coordinate_util.h" + +#include +#include +#include +#include + +#if defined(TECA_HAS_BOOST) +#include +#endif + +using std::string; +using std::vector; +using std::cerr; +using std::endl; +using std::cos; + +//#define TECA_DEBUG + +namespace { +template +void cartesian_iwv(unsigned long nx, unsigned long ny, unsigned long nz, + const coord_t *plev, const num_t *q, num_t *iwv) +{ + unsigned long nxy = nx*ny; + + // initialize the result + memset(iwv, 0, nxy*sizeof(num_t)); + + // work an x-y slice at a time + unsigned long nzm1 = nz - 1; + for (unsigned long k = 0; k < nzm1; ++k) + { + // dp over the slice + num_t h2 = num_t(0.5) * (plev[k+1] - plev[k]); + + // the current two x-y-planes of data + unsigned long knxy = k*nxy; + const num_t *q_k0 = q + knxy; + const num_t *q_k1 = q_k0 + nxy; + + // accumulate this plane of data using trapazoid rule + for (unsigned long i = 0; i < nxy; ++i) + { + iwv[i] += h2 * (q_k0[i] + q_k1[i]); + } + } + + // check the sign, in this way we can handle both increasing and decreasing + // pressure coordinates + num_t s = plev[1] - plev[0] < num_t(0) ? num_t(-1) : num_t(1); + + // scale by -1/g + num_t m1g = s/num_t(9.80665); + for (unsigned long i = 0; i < nxy; ++i) + iwv[i] *= m1g; +} + +template +void cartesian_iwv(unsigned long nx, unsigned long ny, unsigned long nz, + const coord_t *plev, const num_t *q, const char *q_valid, num_t *iwv) +{ + unsigned long nxy = nx*ny; + + // initialize the result + memset(iwv, 0, nxy*sizeof(num_t)); + + // work an x-y slice at a time + unsigned long nzm1 = nz - 1; + for (unsigned long k = 0; k < nzm1; ++k) + { + // dp over the slice + num_t h2 = num_t(0.5) * (plev[k+1] - plev[k]); + + // the current two x-y-planes of data + unsigned long knxy = k*nxy; + const num_t *q_k0 = q + knxy; + const num_t *q_k1 = q_k0 + nxy; + + const char *q_valid_k0 = q_valid + knxy; + const char *q_valid_k1 = q_valid_k0 + nxy; + + // accumulate this plane of data using trapazoid rule + for (unsigned long i = 0; i < nxy; ++i) + { + iwv[i] += ((q_valid_k0[i] && q_valid_k1[i]) ? + h2 * (q_k0[i] + q_k1[i]) : num_t(0)); + } + } + + // check the sign, in this way we can handle both increasing and decreasing + // pressure coordinates + num_t s = plev[1] - plev[0] < num_t(0) ? num_t(-1) : num_t(1); + + // scale by -1/g + num_t m1g = s/num_t(9.80665); + for (unsigned long i = 0; i < nxy; ++i) + iwv[i] *= m1g; +} +} + +// -------------------------------------------------------------------------- +teca_integrated_water_vapor::teca_integrated_water_vapor() : + specific_humidity_variable("Q"), iwv_variable("IWV"), + fill_value(1.0e20) +{ + this->set_number_of_input_connections(1); + this->set_number_of_output_ports(1); +} + +// -------------------------------------------------------------------------- +teca_integrated_water_vapor::~teca_integrated_water_vapor() +{} + +#if defined(TECA_HAS_BOOST) +// -------------------------------------------------------------------------- +void teca_integrated_water_vapor::get_properties_description( + const string &prefix, options_description &global_opts) +{ + options_description opts("Options for " + + (prefix.empty()?"teca_integrated_water_vapor":prefix)); + + opts.add_options() + TECA_POPTS_GET(std::string, prefix, specific_humidity_variable, + "name of the variable containg the specific humidity") + TECA_POPTS_GET(double, prefix, fill_value, + "the value of the NetCDF _FillValue attribute") + ; + + this->teca_algorithm::get_properties_description(prefix, opts); + + global_opts.add(opts); +} + +// -------------------------------------------------------------------------- +void teca_integrated_water_vapor::set_properties( + const string &prefix, variables_map &opts) +{ + this->teca_algorithm::set_properties(prefix, opts); + + TECA_POPTS_SET(opts, std::string, prefix, specific_humidity_variable) + TECA_POPTS_SET(opts, double, prefix, fill_value) +} +#endif + +// -------------------------------------------------------------------------- +teca_metadata teca_integrated_water_vapor::get_output_metadata( + unsigned int port, + const std::vector &input_md) +{ +#ifdef TECA_DEBUG + std::cerr << teca_parallel_id() + << "teca_integrated_water_vapor::get_output_metadata" << std::endl; +#endif + (void)port; + + // set things up in the first pass, and don't modify in subsequent passes + // due to threading concerns + + if (this->get_number_of_derived_variables() == 0) + { + // the base class will handle dealing with the transformation of + // mesh dimensions and reporting the array we produce, but we have + // to determine the data type and tell the name of the produced array. + const teca_metadata &md = input_md[0]; + + teca_metadata attributes; + if (md.get("attributes", attributes)) + { + TECA_ERROR("Failed to determine output data type " + "because attributes are misisng") + return teca_metadata(); + } + + teca_metadata hus_atts; + if (attributes.get(this->specific_humidity_variable, hus_atts)) + { + TECA_ERROR("Failed to determine output data type " + "because attributes for \"" << this->specific_humidity_variable + << "\" are misisng") + return teca_metadata(); + } + + int type_code = 0; + if (hus_atts.get("type_code", type_code)) + { + TECA_ERROR("Failed to determine output data type " + "because attributes for \"" << this->specific_humidity_variable + << "\" is misisng a \"type_code\"") + return teca_metadata(); + } + + teca_array_attributes iwv_atts( + type_code, teca_array_attributes::point_centering, + 0, "kg m^{-1} s^{-1}", "longitudinal integrated vapor transport", + "the longitudinal component of integrated vapor transport", + 1, this->fill_value); + + + // install name and attributes of the output variables in the base classs + this->append_derived_variable(this->iwv_variable); + this->append_derived_variable_attribute(iwv_atts); + } + + if (this->get_number_of_dependent_variables() == 0) + { + // install the names of the input variables in the base class + this->append_dependent_variable(this->specific_humidity_variable); + } + + // invoke the base class method, which does the work of transforming + // the mesh and reporting the variables and their attributes. + return teca_vertical_reduction::get_output_metadata(port, input_md); +} + +// -------------------------------------------------------------------------- +std::vector teca_integrated_water_vapor::get_upstream_request( + unsigned int port, + const std::vector &input_md, + const teca_metadata &request) +{ + // invoke the base class method + return teca_vertical_reduction::get_upstream_request(port, input_md, request); +} + +// -------------------------------------------------------------------------- +const_p_teca_dataset teca_integrated_water_vapor::execute( + unsigned int port, + const std::vector &input_data, + const teca_metadata &request) +{ +#ifdef TECA_DEBUG + std::cerr << teca_parallel_id() + << "teca_integrated_water_vapor::execute" << std::endl; +#endif + (void)port; + + // get the input mesh + const_p_teca_cartesian_mesh in_mesh + = std::dynamic_pointer_cast(input_data[0]); + + if (!in_mesh) + { + TECA_ERROR("Failed to compute IWV because a cartesian mesh is required.") + return nullptr; + } + + // get the input dimensions + unsigned long extent[6] = {0}; + if (in_mesh->get_extent(extent)) + { + TECA_ERROR("Failed to compute IWV because mesh extent is missing.") + return nullptr; + } + + unsigned long nx = extent[1] - extent[0] + 1; + unsigned long ny = extent[3] - extent[2] + 1; + unsigned long nz = extent[5] - extent[4] + 1; + + // get the pressure coordinates + const_p_teca_variant_array p = in_mesh->get_z_coordinates(); + if (!p) + { + TECA_ERROR("Failed to compute IWV because pressure coordinates are missing") + return nullptr; + } + + if (p->size() < 2) + { + TECA_ERROR("Failed to compute IWV because z dimensions " + << p->size() << " < 2 as required by the integration method") + return nullptr; + } + + // gather the input arrays + const_p_teca_variant_array q = + in_mesh->get_point_arrays()->get(this->specific_humidity_variable); + + if (!q) + { + TECA_ERROR("Failed to compute IWV because specific humidity \"" + << this->specific_humidity_variable << "\" is missing") + return nullptr; + } + + const_p_teca_variant_array q_valid = + in_mesh->get_point_arrays()->get(this->specific_humidity_variable + "_valid"); + + // the base class will construct the output mesh + p_teca_cartesian_mesh out_mesh + = std::dynamic_pointer_cast( + std::const_pointer_cast( + teca_vertical_reduction::execute(port, input_data, request))); + + if (!out_mesh) + { + TECA_ERROR("Failed to compute IWV because the output mesh was " + "not constructed") + return nullptr; + } + + // allocate the output arrays + unsigned long nxy = nx*ny; + p_teca_variant_array iwv = q->new_instance(nxy); + + // store the result + out_mesh->get_point_arrays()->set(this->iwv_variable, iwv); + + // calculate IWV + NESTED_TEMPLATE_DISPATCH_FP(const teca_variant_array_impl, + p.get(), _COORDS, + + const NT_COORDS *p_p = static_cast(p.get())->get(); + + NESTED_TEMPLATE_DISPATCH_FP(teca_variant_array_impl, + iwv.get(), _DATA, + + NT_DATA *p_iwv = static_cast(iwv.get())->get(); + + const NT_DATA *p_q = static_cast(q.get())->get(); + + const char *p_q_valid = nullptr; + if (q_valid) + { + using TT_MASK = teca_char_array; + + p_q_valid = dynamic_cast(q_valid.get())->get(); + + ::cartesian_iwv(nx, ny, nz, p_p, p_q, p_q_valid, p_iwv); + } + else + { + ::cartesian_iwv(nx, ny, nz, p_p, p_q, p_iwv); + } + ) + ) + + // pass 2D arrays through. + p_teca_array_collection in_arrays = + std::const_pointer_cast(in_mesh->get_point_arrays()); + + p_teca_array_collection out_arrays = out_mesh->get_point_arrays(); + + int n_arrays = in_arrays->size(); + for (int i = 0; i < n_arrays; ++i) + { + p_teca_variant_array array = in_arrays->get(i); + if (array->size() == nxy) + { + // pass the array. + out_arrays->append(in_arrays->get_name(i), array); + } + } + + return out_mesh; +} diff --git a/alg/teca_integrated_water_vapor.h b/alg/teca_integrated_water_vapor.h new file mode 100644 index 000000000..263ff02a2 --- /dev/null +++ b/alg/teca_integrated_water_vapor.h @@ -0,0 +1,86 @@ +#ifndef teca_integrated_water_vapor_h +#define teca_integrated_water_vapor_h + +#include "teca_shared_object.h" +#include "teca_vertical_reduction.h" +#include "teca_metadata.h" + +#include +#include + +TECA_SHARED_OBJECT_FORWARD_DECL(teca_integrated_water_vapor) + +/// An algorithm that computes integrated water vapor (IWV) +/** + * Compute integrated vaport transport (IWV) from the specific humidity. + * + * IWV = \frac{1}{g} \int_{p_{sfc}}^{p_{top}} q dp + * + * where q is the specific humidity. + * + * This calculation is an instance of a vertical reduction where + * a 3D mesh is transformed into a 2D one. + */ +class teca_integrated_water_vapor : public teca_vertical_reduction +{ +public: + TECA_ALGORITHM_STATIC_NEW(teca_integrated_water_vapor) + TECA_ALGORITHM_DELETE_COPY_ASSIGN(teca_integrated_water_vapor) + TECA_ALGORITHM_CLASS_NAME(teca_integrated_water_vapor) + ~teca_integrated_water_vapor(); + + // report/initialize to/from Boost program options + // objects. + TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION() + TECA_SET_ALGORITHM_PROPERTIES() + + /** @anchor specific_humidity_variable + * @name specific_humidity_variable + * set the name of the variable that contains the specific humidity ("hus") + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, specific_humidity_variable) + ///@} + + /** @anchor iwv_variable + * @name iwv_variable + * set the name of the varaiable that contains the integrated water vapor + * ("iwv"). + */ + ///@{ + TECA_ALGORITHM_PROPERTY(std::string, iwv_variable) + ///@} + + /** @anchor fill_value + * @name fill_value + * set the _fillValue attribute for the output data. default 1.0e20 + */ + ///@{ + TECA_ALGORITHM_PROPERTY(double, fill_value) + ///@} + +protected: + teca_integrated_water_vapor(); + +private: + teca_metadata get_output_metadata( + unsigned int port, + const std::vector &input_md) override; + + std::vector get_upstream_request( + unsigned int port, + const std::vector &input_md, + const teca_metadata &request) override; + + const_p_teca_dataset execute( + unsigned int port, + const std::vector &input_data, + const teca_metadata &request) override; + +private: + std::string specific_humidity_variable; + std::string iwv_variable; + double fill_value; +}; + +#endif diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt index 2d8a1be85..67b18de57 100644 --- a/apps/CMakeLists.txt +++ b/apps/CMakeLists.txt @@ -22,6 +22,9 @@ teca_add_app(teca_bayesian_ar_detect LIBS ${teca_app_link} teca_add_app(teca_integrated_vapor_transport LIBS ${teca_app_link} FEATURES ${TECA_HAS_BOOST} ${TECA_HAS_NETCDF} ${TECA_HAS_UDUNITS}) +teca_add_app(teca_integrated_water_vapor LIBS ${teca_app_link} + FEATURES ${TECA_HAS_BOOST} ${TECA_HAS_NETCDF} ${TECA_HAS_UDUNITS}) + teca_add_app(teca_tc_detect LIBS ${teca_app_link} FEATURES ${TECA_HAS_BOOST} ${TECA_HAS_NETCDF} ${TECA_HAS_UDUNITS}) diff --git a/apps/teca_bayesian_ar_detect.cpp b/apps/teca_bayesian_ar_detect.cpp index d72271f75..ea4d4c368 100644 --- a/apps/teca_bayesian_ar_detect.cpp +++ b/apps/teca_bayesian_ar_detect.cpp @@ -12,6 +12,10 @@ #include "teca_multi_cf_reader.h" #include "teca_integrated_vapor_transport.h" #include "teca_valid_value_mask.h" +#include "teca_cartesian_mesh_source.h" +#include "teca_cartesian_mesh_regrid.h" +#include "teca_indexed_dataset_cache.h" +#include "teca_elevation_mask.h" #include "teca_unpack_data.h" #include "teca_mpi_manager.h" #include "teca_coordinate_util.h" @@ -84,6 +88,17 @@ int main(int argc, char **argv) ("write_ivt", "\nwhen this flag is present IVT vector is written to disk with" " the result\n") + ("dem", value(), "\nA teca_cf_reader regex identifying the" + " file containing surface elevation field or DEM.\n") + ("dem_variable", value()->default_value("Z"), + "\nSets the name of the variable containing the surface elevation field\n") + ("mesh_height", value()->default_value("Zg"), + "\nSets the name of the variable containing the point wise vertical height" + " in meters above mean sea level\n") + + ("ar_probability", value()->default_value("ar_probability"), + "\nSets the name of the variable to store the computed AR probability mask in.\n") + ("ar_weighted_variables", value>()->multitoken(), "\nAn optional list of variables to weight with the computed AR probability." " Each such variable will be multiplied by the computed AR probability, and" @@ -99,8 +114,16 @@ int main(int argc, char **argv) ("periodic_in_x", value()->default_value(1), "\nFlags whether the x dimension (typically longitude) is periodic.\n") - ("binary_ar_threshold", value()->default_value(2.0/3.0,"0.667"), - "\nprobability threshold for segmenting ar_probability to produce ar_binary_tag\n") + ("segment_ar_probability", "\nA flag that enables a binary segmentation of AR" + " probability to be produced. --segment_threshold controls the segmentation." + " threshold and --segment_variable to set the name of the variable to store" + " the result in.\n") + ("segment_threshold", value()->default_value(2.0/3.0,"0.667"), + "\nSets the threshold value that is used when segmenting ar_probability." + " See also --segment_ar_probability\n") + ("segment_variable", value()->default_value("ar_binary_tag"), + "\nSet the name of the variable to store the result of a binary" + " segmentation of AR probabilty. See also --segment_ar_probability.") ("output_file", value()->default_value(std::string("TECA_BARD_%t%.nc")), "\nA path and file name pattern for the output NetCDF files. %t% is replaced with a" @@ -149,6 +172,38 @@ int main(int argc, char **argv) p_teca_multi_cf_reader mcf_reader = teca_multi_cf_reader::New(); mcf_reader->get_properties_description("mcf_reader", advanced_opt_defs); + p_teca_valid_value_mask vv_mask = teca_valid_value_mask::New(); + vv_mask->get_properties_description("vv_mask", advanced_opt_defs); + + p_teca_unpack_data unpack = teca_unpack_data::New(); + unpack->get_properties_description("unpack", advanced_opt_defs); + + p_teca_normalize_coordinates norm_coords = teca_normalize_coordinates::New(); + norm_coords->get_properties_description("norm_coords", advanced_opt_defs); + + p_teca_cf_reader elev_reader = teca_cf_reader::New(); + elev_reader->get_properties_description("elev_reader", advanced_opt_defs); + elev_reader->set_t_axis_variable(""); + + p_teca_normalize_coordinates elev_coords = teca_normalize_coordinates::New(); + elev_coords->get_properties_description("elev_coords", advanced_opt_defs); + elev_coords->set_enable_periodic_shift_x(1); + + p_teca_indexed_dataset_cache elev_cache = teca_indexed_dataset_cache::New(); + elev_cache->get_properties_description("elev_cache", advanced_opt_defs); + elev_cache->set_max_cache_size(1); + + p_teca_cartesian_mesh_source elev_mesh = teca_cartesian_mesh_source::New(); + elev_mesh->get_properties_description("elev_mesh", advanced_opt_defs); + + p_teca_cartesian_mesh_regrid elev_regrid = teca_cartesian_mesh_regrid::New(); + elev_regrid->get_properties_description("elev_regrid", advanced_opt_defs); + + p_teca_elevation_mask elev_mask = teca_elevation_mask::New(); + elev_mask->get_properties_description("elev_mask", advanced_opt_defs); + elev_mask->set_surface_elevation_variable("Z"); + elev_mask->set_mesh_height_variable("ZG"); + p_teca_l2_norm l2_norm = teca_l2_norm::New(); l2_norm->get_properties_description("ivt_magnitude", advanced_opt_defs); l2_norm->set_component_0_variable("IVT_U"); @@ -163,15 +218,6 @@ int main(int argc, char **argv) ivt_int->set_ivt_u_variable("IVT_U"); ivt_int->set_ivt_v_variable("IVT_V"); - p_teca_valid_value_mask vv_mask = teca_valid_value_mask::New(); - vv_mask->get_properties_description("vv_mask", advanced_opt_defs); - - p_teca_unpack_data unpack = teca_unpack_data::New(); - unpack->get_properties_description("unpack", advanced_opt_defs); - - p_teca_normalize_coordinates norm_coords = teca_normalize_coordinates::New(); - norm_coords->get_properties_description("norm_coords", advanced_opt_defs); - // parameter source p_teca_bayesian_ar_detect_parameters params = teca_bayesian_ar_detect_parameters::New(); @@ -182,6 +228,7 @@ int main(int argc, char **argv) p_teca_bayesian_ar_detect ar_detect = teca_bayesian_ar_detect::New(); ar_detect->get_properties_description("ar_detect", advanced_opt_defs); ar_detect->set_ivt_variable("IVT"); + ar_detect->set_ar_probability_variable("ar_probability"); // segment the ar probability field p_teca_binary_segmentation ar_tag = teca_binary_segmentation::New(); @@ -189,7 +236,7 @@ int main(int argc, char **argv) ar_tag->set_threshold_variable("ar_probability"); ar_tag->set_segmentation_variable("ar_binary_tag"); - // mask any requested variables by "ar_probability" + // mask any requested variables by the AR probability p_teca_apply_binary_mask ar_mask = teca_apply_binary_mask::New(); ar_mask->get_properties_description("ar_mask", advanced_opt_defs); ar_mask->set_mask_variable("ar_probability"); @@ -226,11 +273,17 @@ int main(int argc, char **argv) // options will override them cf_reader->set_properties("cf_reader", opt_vals); mcf_reader->set_properties("mcf_reader", opt_vals); - l2_norm->set_properties("ivt_magnitude", opt_vals); - ivt_int->set_properties("ivt_integral", opt_vals); vv_mask->set_properties("vv_mask", opt_vals); - unpack->set_properties("unpack", opt_vals); norm_coords->set_properties("norm_coords", opt_vals); + unpack->set_properties("unpack", opt_vals); + elev_reader->set_properties("elev_reader", opt_vals); + elev_coords->set_properties("elev_coords", opt_vals); + elev_mesh->set_properties("elev_mesh", opt_vals); + elev_cache->set_properties("elev_cache", opt_vals); + elev_regrid->set_properties("elev_regrid", opt_vals); + elev_mask->set_properties("elev_mask", opt_vals); + l2_norm->set_properties("ivt_magnitude", opt_vals); + ivt_int->set_properties("ivt_integral", opt_vals); params->set_properties("parameter_table", opt_vals); ar_detect->set_properties("ar_detect", opt_vals); ar_mask->set_properties("ar_mask", opt_vals); @@ -343,6 +396,7 @@ int main(int argc, char **argv) return -1; } + teca_metadata md; if (do_ivt) { std::string z_var = "plev"; @@ -352,6 +406,50 @@ int main(int argc, char **argv) cf_reader->set_z_axis_variable(z_var); mcf_reader->set_z_axis_variable(z_var); + // add the elevation mask stages + if (opt_vals.count("dem")) + { + if (mpi_man.get_comm_rank() == 0) + TECA_STATUS("Generating elevation mask") + + elev_reader->set_files_regex(opt_vals["dem"].as()); + + elev_coords->set_input_connection(elev_reader->get_output_port()); + + md = head->update_metadata(); + + elev_mesh->set_spatial_bounds(md, false); + elev_mesh->set_spatial_extents(md, false); + elev_mesh->set_x_axis_variable(md); + elev_mesh->set_y_axis_variable(md); + elev_mesh->set_z_axis_variable(md); + elev_mesh->set_t_axis_variable(md); + elev_mesh->set_t_axis(md); + + elev_regrid->set_input_connection(0, elev_mesh->get_output_port()); + elev_regrid->set_input_connection(1, elev_coords->get_output_port()); + + elev_cache->set_input_connection(elev_regrid->get_output_port()); + + elev_mask->set_input_connection(0, head->get_output_port()); + elev_mask->set_input_connection(1, elev_cache->get_output_port()); + + if (!opt_vals["dem_variable"].defaulted()) + elev_mask->set_surface_elevation_variable( + opt_vals["dem_variable"].as()); + + if (!opt_vals["mesh_height"].defaulted()) + elev_mask->set_mesh_height_variable( + opt_vals["mesh_height"].as()); + + elev_mask->set_mask_variables({ + ivt_int->get_specific_humidity_variable() + "_valid", + ivt_int->get_wind_u_variable() + "_valid", + ivt_int->get_wind_v_variable() + "_valid"}); + + head = elev_mask; + } + ivt_int->set_input_connection(head->get_output_port()); l2_norm->set_input_connection(ivt_int->get_output_port()); @@ -367,8 +465,44 @@ int main(int argc, char **argv) ar_detect->set_input_connection(0, params->get_output_port()); ar_detect->set_input_connection(1, head->get_output_port()); - ar_tag->set_input_connection(0, ar_detect->get_output_port()); - head = ar_tag; + if (!opt_vals["ar_probability"].defaulted()) + ar_detect->set_ar_probability_variable( + opt_vals["ar_probability"].as()); + + head = ar_detect; + + // configure binary segmentation + bool do_segment = opt_vals.count("segment_ar_probability"); + if (do_segment) + { + ar_tag->set_input_connection(0, ar_detect->get_output_port()); + + // set input and output variable names + ar_tag->set_threshold_variable( + ar_detect->get_ar_probability_variable()); + + if (!opt_vals["segement_variable"].defaulted()) + ar_tag->set_segmentation_variable( + opt_vals["segment_variable"].as()); + + // set threshold value + double ar_tag_threshold = opt_vals["segment_threshold"].as(); + ar_tag->set_low_threshold_value(ar_tag_threshold); + + // add I/O metadata + teca_metadata seg_atts; + seg_atts.set("long_name", std::string("binary indicator of atmospheric river")); + seg_atts.set("description", std::string("binary indicator of atmospheric river")); + seg_atts.set("scheme", std::string("TECA_BARD")); + seg_atts.set("version", std::string("1.0")); + seg_atts.set("note", + std::string("derived by thresholding ar_probability >= ") + + std::to_string(ar_tag_threshold)); + + ar_tag->set_segmentation_variable_attributes(seg_atts); + + head = ar_tag; + } // configure weight by AR probability std::vector weighted_vars; @@ -378,8 +512,9 @@ int main(int argc, char **argv) weighted_vars = opt_vals["ar_weighted_variables"].as>(); - ar_mask->set_input_connection(0, ar_tag->get_output_port()); + ar_mask->set_input_connection(0, head->get_output_port()); ar_mask->set_masked_variables(weighted_vars); + ar_mask->set_mask_variable(ar_detect->get_ar_probability_variable()); head = ar_mask; } @@ -388,7 +523,9 @@ int main(int argc, char **argv) cf_writer->set_input_connection(head->get_output_port()); // tell the writer to write ivt if needed - std::vector point_arrays({"ar_probability", "ar_binary_tag"}); + std::vector point_arrays( + {ar_detect->get_ar_probability_variable()}); + if ((do_ivt || do_ivt_magnitude) && opt_vals.count("write_ivt_magnitude")) { point_arrays.push_back(l2_norm->get_l2_norm_variable()); @@ -400,6 +537,11 @@ int main(int argc, char **argv) point_arrays.push_back(ivt_int->get_ivt_v_variable()); } + if (do_segment) + { + point_arrays.push_back(ar_tag->get_segmentation_variable()); + } + // tell the writer to write ar weighted variables if needed if (do_weighted) { @@ -450,7 +592,8 @@ int main(int argc, char **argv) if (parse_start_date || parse_end_date) { // run the reporting phase of the pipeline - teca_metadata md = reader->update_metadata(); + if (md.empty()) + md = reader->update_metadata(); teca_metadata atrs; if (md.get("attributes", atrs)) @@ -509,20 +652,6 @@ int main(int argc, char **argv) } } - // set the threshold for calculating ar_binary_tag - double ar_tag_threshold = opt_vals["binary_ar_threshold"].as(); - ar_tag->set_low_threshold_value(ar_tag_threshold); - - // add metadata for ar_binary_tag - teca_metadata seg_atts; - seg_atts.set("long_name", std::string("binary indicator of atmospheric river")); - seg_atts.set("description", std::string("binary indicator of atmospheric river")); - seg_atts.set("scheme", std::string("TECA_BARD")); - seg_atts.set("version", std::string("1.0")); - seg_atts.set("note", - std::string("derived by thresholding ar_probability >= ") + - std::to_string(ar_tag_threshold)); - ar_tag->set_segmentation_variable_attributes(seg_atts); // run the pipeline cf_writer->set_executive(exec); diff --git a/apps/teca_integrated_vapor_transport.cpp b/apps/teca_integrated_vapor_transport.cpp index 86a735e27..15e4d6a97 100644 --- a/apps/teca_integrated_vapor_transport.cpp +++ b/apps/teca_integrated_vapor_transport.cpp @@ -11,6 +11,10 @@ #include "teca_integrated_vapor_transport.h" #include "teca_valid_value_mask.h" #include "teca_unpack_data.h" +#include "teca_cartesian_mesh_source.h" +#include "teca_cartesian_mesh_regrid.h" +#include "teca_elevation_mask.h" +#include "teca_indexed_dataset_cache.h" #include "teca_mpi_manager.h" #include "teca_coordinate_util.h" #include "teca_table.h" @@ -23,6 +27,8 @@ #include #include +#include "teca_cartesian_mesh_writer.h" + using namespace std; using boost::program_options::value; @@ -96,6 +102,16 @@ int main(int argc, char **argv) ("z_axis_variable", value()->default_value("plev"), "\nname of z coordinate variable\n") + ("dem", value(), "\nA teca_cf_reader regex identifying the" + " file containing surface elevation field or DEM.\n") + + ("dem_variable", value()->default_value("Z"), + "\nSets the name of the variable containing the surface elevation field\n") + + ("mesh_height", value()->default_value("Zg"), + "\nSets the name of the variable containing the point wise vertical height" + " in meters above mean sea level\n") + ("first_step", value()->default_value(0), "\nfirst time step to process\n") ("last_step", value()->default_value(-1), "\nlast time step to process\n") @@ -150,6 +166,29 @@ int main(int argc, char **argv) p_teca_unpack_data unpack = teca_unpack_data::New(); unpack->get_properties_description("unpack", advanced_opt_defs); + p_teca_cf_reader elev_reader = teca_cf_reader::New(); + elev_reader->get_properties_description("elev_reader", advanced_opt_defs); + elev_reader->set_t_axis_variable(""); + + p_teca_normalize_coordinates elev_coords = teca_normalize_coordinates::New(); + elev_coords->get_properties_description("elev_coords", advanced_opt_defs); + elev_coords->set_enable_periodic_shift_x(1); + + p_teca_indexed_dataset_cache elev_cache = teca_indexed_dataset_cache::New(); + elev_cache->get_properties_description("elev_cache", advanced_opt_defs); + elev_cache->set_max_cache_size(1); + + p_teca_cartesian_mesh_source elev_mesh = teca_cartesian_mesh_source::New(); + elev_mesh->get_properties_description("elev_mesh", advanced_opt_defs); + + p_teca_cartesian_mesh_regrid elev_regrid = teca_cartesian_mesh_regrid::New(); + elev_regrid->get_properties_description("elev_regrid", advanced_opt_defs); + + p_teca_elevation_mask elev_mask = teca_elevation_mask::New(); + elev_mask->get_properties_description("elev_mask", advanced_opt_defs); + elev_mask->set_surface_elevation_variable("Z"); + elev_mask->set_mesh_height_variable("ZG"); + p_teca_integrated_vapor_transport ivt_int = teca_integrated_vapor_transport::New(); ivt_int->get_properties_description("ivt_integral", advanced_opt_defs); ivt_int->set_specific_humidity_variable("Q"); @@ -197,6 +236,12 @@ int main(int argc, char **argv) norm_coords->set_properties("norm_coords", opt_vals); vv_mask->set_properties("vv_mask", opt_vals); unpack->set_properties("unpack", opt_vals); + elev_reader->set_properties("elev_reader", opt_vals); + elev_coords->set_properties("elev_coords", opt_vals); + elev_mesh->set_properties("elev_mesh", opt_vals); + elev_cache->set_properties("elev_cache", opt_vals); + elev_regrid->set_properties("elev_regrid", opt_vals); + elev_mask->set_properties("elev_mask", opt_vals); ivt_int->set_properties("ivt_integral", opt_vals); l2_norm->set_properties("ivt_magnitude", opt_vals); cf_writer->set_properties("cf_writer", opt_vals); @@ -209,6 +254,16 @@ int main(int argc, char **argv) // configure the reader bool have_file = opt_vals.count("input_file"); bool have_regex = opt_vals.count("input_regex"); + if ((have_file && have_regex) || !(have_file || have_regex)) + { + if (mpi_man.get_comm_rank() == 0) + { + TECA_ERROR("Extacly one of --input_file or --input_regex can be specified. " + "Use --input_file to activate the multi_cf_reader (HighResMIP datasets) " + "and --input_regex to activate the cf_reader (CAM like datasets)") + } + return -1; + } if (opt_vals.count("input_file")) { @@ -234,6 +289,13 @@ int main(int argc, char **argv) mcf_reader->set_y_axis_variable(opt_vals["y_axis_variable"].as()); } + std::string z_var = "plev"; + if (!opt_vals["z_axis_variable"].defaulted()) + z_var = opt_vals["z_axis_variable"].as(); + + cf_reader->set_z_axis_variable(z_var); + mcf_reader->set_z_axis_variable(z_var); + // set the inputs to the integrator if (!opt_vals["wind_u"].defaulted()) { @@ -278,18 +340,80 @@ int main(int argc, char **argv) // add the ivt caluation stages if needed bool do_ivt = opt_vals["write_ivt"].as(); bool do_ivt_magnitude = opt_vals["write_ivt_magnitude"].as(); + if (!(do_ivt || do_ivt_magnitude)) - std::string z_var = "plev"; - if (!opt_vals["z_axis_variable"].defaulted()) - z_var = opt_vals["z_axis_variable"].as(); + { + if (mpi_man.get_comm_rank() == 0) + { + TECA_ERROR("At least one of --write_ivt or --write_ivt_magnitude " + " must be set.") + } + return -1; + } - cf_reader->set_z_axis_variable(z_var); - mcf_reader->set_z_axis_variable(z_var); + // add the elevation mask stages + teca_metadata md; + if (opt_vals.count("dem")) + { + if (mpi_man.get_comm_rank() == 0) + TECA_STATUS("Generating elevation mask") + + elev_reader->set_files_regex(opt_vals["dem"].as()); + + elev_coords->set_input_connection(elev_reader->get_output_port()); + + md = head->update_metadata(); + + elev_mesh->set_spatial_bounds(md, false); + elev_mesh->set_spatial_extents(md, false); + elev_mesh->set_x_axis_variable(md); + elev_mesh->set_y_axis_variable(md); + elev_mesh->set_z_axis_variable(md); + elev_mesh->set_t_axis_variable(md); + elev_mesh->set_t_axis(md); + + elev_regrid->set_input_connection(0, elev_mesh->get_output_port()); + elev_regrid->set_input_connection(1, elev_coords->get_output_port()); + + elev_cache->set_input_connection(elev_regrid->get_output_port()); + + /*p_teca_cartesian_mesh_writer rdw = teca_cartesian_mesh_writer::New(); + rdw->set_input_connection(elev_cache->get_output_port()); + rdw->set_file_name("regrid_dem_%t%.vtk");*/ + + elev_mask->set_input_connection(0, head->get_output_port()); + elev_mask->set_input_connection(1, elev_cache->get_output_port()); + //elev_mask->set_input_connection(1, rdw->get_output_port()); + + if (!opt_vals["dem_variable"].defaulted()) + elev_mask->set_surface_elevation_variable( + opt_vals["dem_variable"].as()); + + if (!opt_vals["mesh_height"].defaulted()) + elev_mask->set_mesh_height_variable( + opt_vals["mesh_height"].as()); + + elev_mask->set_mask_variables({ + ivt_int->get_specific_humidity_variable() + "_valid", + ivt_int->get_wind_u_variable() + "_valid", + ivt_int->get_wind_v_variable() + "_valid"}); + + /*p_teca_cartesian_mesh_writer emw = teca_cartesian_mesh_writer::New(); + emw->set_input_connection(elev_mask->get_output_port()); + emw->set_file_name("elev_mask_%t%.vtk"); + emw->set_binary(1); + head = emw;*/ + + head = elev_mask; + } ivt_int->set_input_connection(head->get_output_port()); if (do_ivt_magnitude) { + if (mpi_man.get_comm_rank() == 0) + TECA_STATUS("Computing IVT magnitude") + l2_norm->set_input_connection(ivt_int->get_output_port()); head = l2_norm; } @@ -327,27 +451,6 @@ int main(int argc, char **argv) cf_writer->set_thread_pool_size(opt_vals["n_threads"].as()); // some minimal check for missing options - if ((have_file && have_regex) || !(have_file || have_regex)) - { - if (mpi_man.get_comm_rank() == 0) - { - TECA_ERROR("Extacly one of --input_file or --input_regex can be specified. " - "Use --input_file to activate the multi_cf_reader (HighResMIP datasets) " - "and --input_regex to activate the cf_reader (CAM like datasets)") - } - return -1; - } - - if (!(do_ivt || do_ivt_magnitude)) - { - if (mpi_man.get_comm_rank() == 0) - { - TECA_ERROR("AT least one of --write_ivt or --write_ivt_magnitude " - " must be set.") - } - return -1; - } - if (cf_writer->get_file_name().empty()) { if (mpi_man.get_comm_rank() == 0) @@ -367,7 +470,8 @@ int main(int argc, char **argv) if (parse_start_date || parse_end_date) { // run the reporting phase of the pipeline - teca_metadata md = reader->update_metadata(); + if (md.empty()) + md = reader->update_metadata(); teca_metadata atrs; if (md.get("attributes", atrs)) diff --git a/apps/teca_integrated_water_vapor.cpp b/apps/teca_integrated_water_vapor.cpp new file mode 100644 index 000000000..bd90145fb --- /dev/null +++ b/apps/teca_integrated_water_vapor.cpp @@ -0,0 +1,440 @@ +#include "teca_config.h" +#include "teca_cf_reader.h" +#include "teca_cf_writer.h" +#include "teca_index_executive.h" +#include "teca_normalize_coordinates.h" +#include "teca_metadata.h" +#include "teca_integrated_water_vapor.h" +#include "teca_binary_segmentation.h" +#include "teca_l2_norm.h" +#include "teca_multi_cf_reader.h" +#include "teca_integrated_water_vapor.h" +#include "teca_valid_value_mask.h" +#include "teca_unpack_data.h" +#include "teca_cartesian_mesh_source.h" +#include "teca_cartesian_mesh_regrid.h" +#include "teca_elevation_mask.h" +#include "teca_indexed_dataset_cache.h" +#include "teca_mpi_manager.h" +#include "teca_coordinate_util.h" +#include "teca_table.h" +#include "teca_dataset_source.h" +#include "teca_app_util.h" +#include "calcalcs.h" + +#include +#include +#include +#include + +#include "teca_cartesian_mesh_writer.h" + +using namespace std; + +using boost::program_options::value; + +// -------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + // initialize mpi + teca_mpi_manager mpi_man(argc, argv); + + // initialize command line options description + // set up some common options to simplify use for most + // common scenarios + int help_width = 100; + options_description basic_opt_defs( + "Basic usage:\n\n" + "The following options are the most commonly used. Information\n" + "on advanced options can be displayed using --advanced_help\n\n" + "Basic command line options", help_width, help_width - 4 + ); + basic_opt_defs.add_options() + ("input_file", value(), "\na teca_multi_cf_reader configuration file" + " identifying the set of NetCDF CF2 files to process. When present data is" + " read using the teca_multi_cf_reader. Use one of either --input_file or" + " --input_regex.\n") + + ("input_regex", value(), "\na teca_cf_reader regex identifying the" + " set of NetCDF CF2 files to process. When present data is read using the" + " teca_cf_reader. Use one of either --input_file or --input_regex.\n") + + ("specific_humidity", value()->default_value("Q"), + "\nname of variable with the 3D specific humidity field.\n") + + ("iwv", value()->default_value("IWV"), + "\nname to use for the longitudinal component of the integrated vapor transport vector.\n") + + ("output_file", value()->default_value("IWV_%t%.nc"), + "\nA path and file name pattern for the output NetCDF files. %t% is replaced with a" + " human readable date and time corresponding to the time of the first time step in" + " the file. Use --cf_writer::date_format to change the formatting\n") + + ("steps_per_file", value()->default_value(128), + "\nnumber of time steps per output file\n") + + ("x_axis_variable", value()->default_value("lon"), + "\nname of x coordinate variable\n") + ("y_axis_variable", value()->default_value("lat"), + "\nname of y coordinate variable\n") + ("z_axis_variable", value()->default_value("plev"), + "\nname of z coordinate variable\n") + + ("dem", value(), "\nA teca_cf_reader regex identifying the" + " file containing surface elevation field or DEM.\n") + + ("dem_variable", value()->default_value("Z"), + "\nSets the name of the variable containing the surface elevation field\n") + + ("mesh_height", value()->default_value("Zg"), + "\nSets the name of the variable containing the point wise vertical height" + " in meters above mean sea level\n") + + ("first_step", value()->default_value(0), "\nfirst time step to process\n") + ("last_step", value()->default_value(-1), "\nlast time step to process\n") + + ("start_date", value(), "\nThe first time to process in 'Y-M-D h:m:s'" + " format. Note: There must be a space between the date and time specification\n") + ("end_date", value(), "\nThe last time to process in 'Y-M-D h:m:s' format\n") + + ("n_threads", value()->default_value(-1), "\nSets the thread pool size on each MPI" + " rank. When the default value of -1 is used TECA will coordinate the thread pools" + " across ranks such each thread is bound to a unique physical core.\n") + + ("verbose", "\nenable extra terminal output\n") + + ("help", "\ndisplays documentation for application specific command line options\n") + ("advanced_help", "\ndisplays documentation for algorithm specific command line options\n") + ("full_help", "\ndisplays both basic and advanced documentation together\n") + ; + + // add all options from each pipeline stage for more advanced use + options_description advanced_opt_defs( + "Advanced usage:\n\n" + "The following list contains the full set options giving one full\n" + "control over all runtime modifiable parameters. The basic options\n" + "(see" "--help) map to these, and will override them if both are\n" + "specified.\n\n" + "integrated vapor transport pipeline:\n\n" + " (cf / mcf_reader)\n" + " \\\n" + " (iwv_integral)--(iwv_magnitude)\n" + " \\\n" + " (cf_writer)\n\n" + "Advanced command line options", help_width, help_width - 4 + ); + + // create the pipeline stages here, they contain the + // documentation and parse command line. + // objects report all of their properties directly + // set default options here so that command line options override + // them. while we are at it connect the pipeline + p_teca_cf_reader cf_reader = teca_cf_reader::New(); + cf_reader->get_properties_description("cf_reader", advanced_opt_defs); + + p_teca_multi_cf_reader mcf_reader = teca_multi_cf_reader::New(); + mcf_reader->get_properties_description("mcf_reader", advanced_opt_defs); + + p_teca_normalize_coordinates norm_coords = teca_normalize_coordinates::New(); + norm_coords->get_properties_description("norm_coords", advanced_opt_defs); + + p_teca_valid_value_mask vv_mask = teca_valid_value_mask::New(); + vv_mask->get_properties_description("vv_mask", advanced_opt_defs); + + p_teca_unpack_data unpack = teca_unpack_data::New(); + unpack->get_properties_description("unpack", advanced_opt_defs); + + p_teca_cf_reader elev_reader = teca_cf_reader::New(); + elev_reader->get_properties_description("elev_reader", advanced_opt_defs); + elev_reader->set_t_axis_variable(""); + + p_teca_normalize_coordinates elev_coords = teca_normalize_coordinates::New(); + elev_coords->get_properties_description("elev_coords", advanced_opt_defs); + elev_coords->set_enable_periodic_shift_x(1); + + p_teca_indexed_dataset_cache elev_cache = teca_indexed_dataset_cache::New(); + elev_cache->get_properties_description("elev_cache", advanced_opt_defs); + elev_cache->set_max_cache_size(1); + + p_teca_cartesian_mesh_source elev_mesh = teca_cartesian_mesh_source::New(); + elev_mesh->get_properties_description("elev_mesh", advanced_opt_defs); + + p_teca_cartesian_mesh_regrid elev_regrid = teca_cartesian_mesh_regrid::New(); + elev_regrid->get_properties_description("elev_regrid", advanced_opt_defs); + + p_teca_elevation_mask elev_mask = teca_elevation_mask::New(); + elev_mask->get_properties_description("elev_mask", advanced_opt_defs); + elev_mask->set_surface_elevation_variable("Z"); + elev_mask->set_mesh_height_variable("ZG"); + + p_teca_integrated_water_vapor iwv_int = teca_integrated_water_vapor::New(); + iwv_int->get_properties_description("iwv_integral", advanced_opt_defs); + iwv_int->set_specific_humidity_variable("Q"); + iwv_int->set_iwv_variable("IWV"); + + // Add an executive for the writer + p_teca_index_executive exec = teca_index_executive::New(); + + // Add the writer + p_teca_cf_writer cf_writer = teca_cf_writer::New(); + cf_writer->get_properties_description("cf_writer", advanced_opt_defs); + cf_writer->set_verbose(0); + cf_writer->set_steps_per_file(128); + + // package basic and advanced options for display + options_description all_opt_defs(help_width, help_width - 4); + all_opt_defs.add(basic_opt_defs).add(advanced_opt_defs); + + // parse the command line + int ierr = 0; + variables_map opt_vals; + if ((ierr = teca_app_util::process_command_line_help( + mpi_man.get_comm_rank(), argc, argv, basic_opt_defs, + advanced_opt_defs, all_opt_defs, opt_vals))) + { + if (ierr == 1) + return 0; + return -1; + } + + // pass command line arguments into the pipeline objects + // advanced options are processed first, so that the basic + // options will override them + cf_reader->set_properties("cf_reader", opt_vals); + mcf_reader->set_properties("mcf_reader", opt_vals); + norm_coords->set_properties("norm_coords", opt_vals); + vv_mask->set_properties("vv_mask", opt_vals); + unpack->set_properties("unpack", opt_vals); + elev_reader->set_properties("elev_reader", opt_vals); + elev_coords->set_properties("elev_coords", opt_vals); + elev_mesh->set_properties("elev_mesh", opt_vals); + elev_cache->set_properties("elev_cache", opt_vals); + elev_regrid->set_properties("elev_regrid", opt_vals); + elev_mask->set_properties("elev_mask", opt_vals); + iwv_int->set_properties("iwv_integral", opt_vals); + cf_writer->set_properties("cf_writer", opt_vals); + + // now pass in the basic options, these are processed + // last so that they will take precedence + // configure the pipeline from the command line options. + p_teca_algorithm reader; + + // configure the reader + bool have_file = opt_vals.count("input_file"); + bool have_regex = opt_vals.count("input_regex"); + if ((have_file && have_regex) || !(have_file || have_regex)) + { + if (mpi_man.get_comm_rank() == 0) + { + TECA_ERROR("Extacly one of --input_file or --input_regex can be specified. " + "Use --input_file to activate the multi_cf_reader (HighResMIP datasets) " + "and --input_regex to activate the cf_reader (CAM like datasets)") + } + return -1; + } + + if (opt_vals.count("input_file")) + { + mcf_reader->set_input_file(opt_vals["input_file"].as()); + reader = mcf_reader; + } + else if (opt_vals.count("input_regex")) + { + cf_reader->set_files_regex(opt_vals["input_regex"].as()); + reader = cf_reader; + } + p_teca_algorithm head = reader; + + if (!opt_vals["x_axis_variable"].defaulted()) + { + cf_reader->set_x_axis_variable(opt_vals["x_axis_variable"].as()); + mcf_reader->set_x_axis_variable(opt_vals["x_axis_variable"].as()); + } + + if (!opt_vals["y_axis_variable"].defaulted()) + { + cf_reader->set_y_axis_variable(opt_vals["y_axis_variable"].as()); + mcf_reader->set_y_axis_variable(opt_vals["y_axis_variable"].as()); + } + + std::string z_var = "plev"; + if (!opt_vals["z_axis_variable"].defaulted()) + z_var = opt_vals["z_axis_variable"].as(); + + cf_reader->set_z_axis_variable(z_var); + mcf_reader->set_z_axis_variable(z_var); + + // set the inputs to the integrator + if (!opt_vals["specific_humidity"].defaulted()) + { + iwv_int->set_specific_humidity_variable( + opt_vals["specific_humidity"].as()); + } + + // set all that use or produce iwv + if (!opt_vals["iwv"].defaulted()) + iwv_int->set_iwv_variable(opt_vals["iwv"].as()); + + // add the valid value mask stage + norm_coords->set_input_connection(head->get_output_port()); + vv_mask->set_input_connection(norm_coords->get_output_port()); + unpack->set_input_connection(vv_mask->get_output_port()); + head = unpack; + + // add the elevation mask stages + teca_metadata md; + if (opt_vals.count("dem")) + { + if (mpi_man.get_comm_rank() == 0) + TECA_STATUS("Generating elevation mask") + + elev_reader->set_files_regex(opt_vals["dem"].as()); + + elev_coords->set_input_connection(elev_reader->get_output_port()); + + md = head->update_metadata(); + + elev_mesh->set_spatial_bounds(md, false); + elev_mesh->set_spatial_extents(md, false); + elev_mesh->set_x_axis_variable(md); + elev_mesh->set_y_axis_variable(md); + elev_mesh->set_z_axis_variable(md); + elev_mesh->set_t_axis_variable(md); + elev_mesh->set_t_axis(md); + + elev_regrid->set_input_connection(0, elev_mesh->get_output_port()); + elev_regrid->set_input_connection(1, elev_coords->get_output_port()); + + elev_cache->set_input_connection(elev_regrid->get_output_port()); + + elev_mask->set_input_connection(0, head->get_output_port()); + elev_mask->set_input_connection(1, elev_cache->get_output_port()); + + if (!opt_vals["dem_variable"].defaulted()) + elev_mask->set_surface_elevation_variable( + opt_vals["dem_variable"].as()); + + if (!opt_vals["mesh_height"].defaulted()) + elev_mask->set_mesh_height_variable( + opt_vals["mesh_height"].as()); + + elev_mask->set_mask_variables({ + iwv_int->get_specific_humidity_variable() + "_valid"}); + + head = elev_mask; + } + + iwv_int->set_input_connection(head->get_output_port()); + head = iwv_int; + + // tell the writer to write iwv if needed + std::vector point_arrays; + point_arrays.push_back(iwv_int->get_iwv_variable()); + + cf_writer->set_point_arrays(point_arrays); + + cf_writer->set_file_name(opt_vals["output_file"].as()); + + if (!opt_vals["steps_per_file"].defaulted()) + cf_writer->set_steps_per_file(opt_vals["steps_per_file"].as()); + + if (!opt_vals["first_step"].defaulted()) + cf_writer->set_first_step(opt_vals["first_step"].as()); + + if (!opt_vals["last_step"].defaulted()) + cf_writer->set_last_step(opt_vals["last_step"].as()); + + if (opt_vals.count("verbose")) + { + cf_writer->set_verbose(1); + exec->set_verbose(1); + } + + cf_writer->set_thread_pool_size(opt_vals["n_threads"].as()); + + // some minimal check for missing options + if (cf_writer->get_file_name().empty()) + { + if (mpi_man.get_comm_rank() == 0) + { + TECA_ERROR("missing file name pattern for netcdf writer. " + "See --help for a list of command line options.") + } + return -1; + } + + // connect the fixed stages of the pipeline + cf_writer->set_input_connection(head->get_output_port()); + + // look for requested time step range, start + bool parse_start_date = opt_vals.count("start_date"); + bool parse_end_date = opt_vals.count("end_date"); + if (parse_start_date || parse_end_date) + { + // run the reporting phase of the pipeline + if (md.empty()) + md = reader->update_metadata(); + + teca_metadata atrs; + if (md.get("attributes", atrs)) + { + TECA_ERROR("metadata missing attributes") + return -1; + } + + teca_metadata time_atts; + std::string calendar; + std::string units; + if (atrs.get("time", time_atts) + || time_atts.get("calendar", calendar) + || time_atts.get("units", units)) + { + TECA_ERROR("failed to determine the calendaring parameters") + return -1; + } + + teca_metadata coords; + p_teca_variant_array time; + if (md.get("coordinates", coords) || !(time = coords.get("t"))) + { + TECA_ERROR("failed to determine time coordinate") + return -1; + } + + // convert date string to step, start date + if (parse_start_date) + { + unsigned long first_step = 0; + std::string start_date = opt_vals["start_date"].as(); + if (teca_coordinate_util::time_step_of(time, true, true, calendar, + units, start_date, first_step)) + { + TECA_ERROR("Failed to locate time step for start date \"" + << start_date << "\"") + return -1; + } + cf_writer->set_first_step(first_step); + } + + // and end date + if (parse_end_date) + { + unsigned long last_step = 0; + std::string end_date = opt_vals["end_date"].as(); + if (teca_coordinate_util::time_step_of(time, false, true, calendar, + units, end_date, last_step)) + { + TECA_ERROR("Failed to locate time step for end date \"" + << end_date << "\"") + return -1; + } + cf_writer->set_last_step(last_step); + } + } + + // run the pipeline + cf_writer->set_executive(exec); + cf_writer->update(); + + return 0; +} diff --git a/apps/teca_tc_wind_radii.cpp b/apps/teca_tc_wind_radii.cpp index 0dcbaae9b..8654d748a 100644 --- a/apps/teca_tc_wind_radii.cpp +++ b/apps/teca_tc_wind_radii.cpp @@ -42,7 +42,7 @@ int main(int argc, char **argv) "Basic command line options", help_width, help_width - 4 ); basic_opt_defs.add_options() - ("track_file", value(), "\na file containing cyclone tracks (tracks.bin)\n") + ("track_file", value(), "\na file containing cyclone tracks\n") ("input_file", value(), "\na teca_multi_cf_reader configuration file" " identifying the set of NetCDF CF2 files to process. When present data is" @@ -160,12 +160,33 @@ int main(int argc, char **argv) // now pass in the basic options, these are processed // last so that they will take precedence - if (!opt_vals["track_file"].defaulted()) - track_reader->set_file_name(opt_vals["track_file"].as()); + if (!opt_vals.count("track_file")) + { + if (mpi_man.get_comm_rank() == 0) + { + TECA_ERROR("A file with previously calculated storm tracks must be " + "specified with --track_file") + } + return -1; + } + + track_reader->set_file_name(opt_vals["track_file"].as()); bool have_file = opt_vals.count("input_file"); bool have_wind_files = opt_vals.count("wind_files"); bool have_regex = opt_vals.count("input_regex"); + + if ((have_file && have_regex) || !(have_file || have_regex)) + { + if (mpi_man.get_comm_rank() == 0) + { + TECA_ERROR("Extacly one of --input_file or --input_regex can be specified. " + "Use --input_file to activate the multi_cf_reader (HighResMIP datasets) " + "and --input_regex to activate the cf_reader (CAM like datasets)") + } + return -1; + } + p_teca_algorithm wind_reader; if (have_file) { @@ -240,18 +261,6 @@ int main(int argc, char **argv) else map_reduce->set_thread_pool_size(-1); - // some minimal check for missing options - if ((have_file && have_regex) || !(have_file || have_regex)) - { - if (mpi_man.get_comm_rank() == 0) - { - TECA_ERROR("Extacly one of --input_file or --input_regex can be specified. " - "Use --input_file to activate the multi_cf_reader (HighResMIP datasets) " - "and --input_regex to activate the cf_reader (CAM like datasets)") - } - return -1; - } - // connect the pipeline wind_coords->set_input_connection(wind_reader->get_output_port()); wind_radii->set_input_connection(0, track_input->get_output_port()); diff --git a/doc/rtd/applications.rst b/doc/rtd/applications.rst index e01f6a515..5955188f7 100644 --- a/doc/rtd/applications.rst +++ b/doc/rtd/applications.rst @@ -21,6 +21,8 @@ batch script rather than in your shell. +----------------------------------------+--------------------------------------------------+ | :ref:`teca_integrated_vapor_transport` | Computes IVT (integrated vapor transport) | +----------------------------------------+--------------------------------------------------+ +| :ref:`teca_integrated_water_vapor` | Computes IWV (integrated water vapor) | ++----------------------------------------+--------------------------------------------------+ | :ref:`teca_bayesian_ar_detect` | AR detection with uncertainty quantification | +----------------------------------------+--------------------------------------------------+ | :ref:`teca_deeplab_ar_detect` | A machine learning based AR detector | @@ -1002,6 +1004,19 @@ Command Line Arguments --write_ivt when this flag is present IVT vector is written to disk with the result +--dem arg + A teca_cf_reader regex identifying the file containing surface elevation field or DEM. + +--dem_variable arg (=Z) + Sets the name of the variable containing the surface elevation field + +--mesh_height arg (=Zg) + Sets the name of the variable containing the point wise vertical height in meters above mean + sea level + +--ar_probability arg (=ar_probability) + Sets the name of the variable to store the computed AR probability mask in. + --ar_weighted_variables arg An optional list of variables to weight with the computed AR probability. Each such variable will be multiplied by the computed AR probability, and written to disk as "NAME_ar_wgtd". @@ -1018,10 +1033,20 @@ Command Line Arguments --periodic_in_x arg (=1) Flags whether the x dimension (typically longitude) is periodic. ---binary_ar_threshold arg (=0.667) - probability threshold for segmenting ar_probability to produce ar_binary_tag +--segment_ar_probability + A flag that enables a binary segmentation of AR probability to be produced. `--segment_threshold` + controls the segmentation. threshold and `--segment_variable` to set the name of the variable to + store the result in. + +--segment_threshold arg (=0.667) + Sets the threshold value that is used when segmenting ar_probability. See also + `--segment_ar_probability` + +--segment_variable arg (=ar_binary_tag) + Set the name of the variable to store the result of a binary segmentation of AR probabilty. See + also `--segment_ar_probability`. ---output_file arg (=CASCADE_BARD_%t%.nc) +--output_file arg (=TECA_BARD_%t%.nc) A path and file name pattern for the output NetCDF files. %t% is replaced with a human readable date and time corresponding to the time of the first time step in the file. Use `--cf_writer::date_format` to change the formatting @@ -1029,6 +1054,8 @@ Command Line Arguments --steps_per_file arg (=128) number of time steps per output file +--first_step arg (=0) + first time step to process --last_step arg (=-1) last time step to process @@ -1194,8 +1221,14 @@ taken into account. See the :ref:`teca_metadata_probe` ARTMIP :ref:`example` was used to configure the readers. + +.. _teca_integrated_water_vapor: + +teca_integrated_water_vapor +------------------------------- +The integrated water vapor(IWV) command line application computes: + +.. math:: + + IWV = \frac{1}{g} \int_{p_{sfc}}^{p_{top}} q dp + +where g is the acceleration due to Earth's gravity, p is atmospheric pressure, +and q is specific humidity. + + +Inputs +~~~~~~ +A 3D time dependent mesh in NetCDF CF2 format with: + +1. specific humidity + +Outputs +~~~~~~~ +A 2D mesh with: + +1. IWV + + +Command Line Arguments +~~~~~~~~~~~~~~~~~~~~~~ +--input_file arg + a teca_multi_cf_reader configuration file identifying the set of NetCDF CF2 files to process. + When present data is read using the teca_multi_cf_reader. Use one of either `--input_file` or + `--input_regex`. + +--input_regex arg + a teca_cf_reader regex identifying the set of NetCDF CF2 files to process. When present data is + read using the teca_cf_reader. Use one of either `--input_file` or `--input_regex`. + +--specific_humidity arg (=Q) + name of variable with the 3D specific humidity field. + +--iwv arg (=IWV) + name to use for the longitudinal component of the integrated vapor transport vector. + +--output_file arg (=IWV_%t%.nc) + A path and file name pattern for the output NetCDF files. %t% is replaced with a human readable + date and time corresponding to the time of the first time step in the file. Use + `--cf_writer::date_format` to change the formatting + +--steps_per_file arg (=128) + number of time steps per output file + +--x_axis_variable arg (=lon) + name of x coordinate variable + +--y_axis_variable arg (=lat) + name of y coordinate variable + +--z_axis_variable arg (=plev) + name of z coordinate variable + +--dem arg + A teca_cf_reader regex identifying the file containing surface elevation field or DEM. + +--dem_variable arg (=Z) + Sets the name of the variable containing the surface elevation field + +--mesh_height arg (=Zg) + Sets the name of the variable containing the point wise vertical height in meters above mean + sea level + +--first_step arg (=0) + first time step to process + +--last_step arg (=-1) + last time step to process + +--start_date arg + The first time to process in 'Y-M-D h:m:s' format. Note: There must be a space between the date + and time specification + +--end_date arg + The last time to process in 'Y-M-D h:m:s' format + +--n_threads arg (=-1) + Sets the thread pool size on each MPI rank. When the default value of -1 is used TECA will + coordinate the thread pools across ranks such each thread is bound to a unique physical core. + +--verbose + enable extra terminal output + +--help + displays documentation for application specific command line options + +--advanced_help + displays documentation for algorithm specific command line options + +--full_help + displays both basic and advanced documentation together + + .. _teca_tc_detect: teca_tc_detect diff --git a/io/teca_cf_layout_manager.cxx b/io/teca_cf_layout_manager.cxx index 2f99c5a87..d049f462e 100644 --- a/io/teca_cf_layout_manager.cxx +++ b/io/teca_cf_layout_manager.cxx @@ -40,7 +40,10 @@ int teca_cf_layout_manager::create(const std::string &file_name, } coords.get("t_variable", this->t_variable); - p_teca_variant_array t = coords.get("t"); + + p_teca_variant_array t; + if (!this->t_variable.empty()) + t = coords.get("t"); // construct the file name if (!date_format.empty()) @@ -593,7 +596,11 @@ int teca_cf_layout_manager::write(long index, { size_t starts[4] = {0, 0, 0, 0}; size_t counts[4] = {1, 0, 0, 0}; - for (int i = 1; i < this->n_dims; ++i) + + // make space for the time dimension + int i0 = this->t ? 1 : 0; + + for (int i = i0; i < this->n_dims; ++i) counts[i] = this->dims[i]; // get this data's position in the file diff --git a/python/teca_py_alg.i b/python/teca_py_alg.i index 97bc81bda..80367d2af 100644 --- a/python/teca_py_alg.i +++ b/python/teca_py_alg.i @@ -19,6 +19,7 @@ #include "teca_derived_quantity_numerics.h" #include "teca_descriptive_statistics.h" #include "teca_evaluate_expression.h" +#include "teca_elevation_mask.h" #include "teca_integrated_vapor_transport.h" #include "teca_l2_norm.h" #include "teca_laplacian.h" @@ -420,3 +421,11 @@ struct teca_tc_saffir_simpson %shared_ptr(teca_unpack_data) %ignore teca_unpack_data::operator=; %include "teca_unpack_data.h" + +/*************************************************************************** + elevation_mask + ***************************************************************************/ +%ignore teca_elevation_mask::shared_from_this; +%shared_ptr(teca_elevation_mask) +%ignore teca_elevation_mask::operator=; +%include "teca_elevation_mask.h" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index cac540072..be18a6705 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -646,10 +646,28 @@ teca_add_test(test_cf_writer_bad_type REQ_TECA_DATA WILL_FAIL) -teca_add_test(test_integrated_vapor_transport +teca_add_test(test_integrated_water_vapor_vv_mask + SOURCES test_integrated_water_vapor.cpp + EXEC_NAME test_integrated_water_vapor + LIBS teca_core teca_data teca_alg teca_io ${teca_test_link} + COMMAND test_integrated_water_vapor 0 0 0 + FEATURES ${TECA_HAS_NETCDF}) + +teca_add_test(test_integrated_water_vapor_elev_mask + LIBS teca_core teca_data teca_alg teca_io ${teca_test_link} + COMMAND test_integrated_water_vapor 1 0 0 + FEATURES ${TECA_HAS_NETCDF}) + +teca_add_test(test_integrated_vapor_transport_vv_mask SOURCES test_integrated_vapor_transport.cpp + EXEC_NAME test_integrated_vapor_transport LIBS teca_core teca_data teca_alg teca_io ${teca_test_link} - COMMAND test_integrated_vapor_transport + COMMAND test_integrated_vapor_transport 0 0 0 + FEATURES ${TECA_HAS_NETCDF}) + +teca_add_test(test_integrated_vapor_transport_elev_mask + LIBS teca_core teca_data teca_alg teca_io ${teca_test_link} + COMMAND test_integrated_vapor_transport 1 0 0 FEATURES ${TECA_HAS_NETCDF}) teca_add_test(test_cf_time_axis_reader @@ -689,3 +707,13 @@ teca_add_test(test_bounds_to_extent SOURCES test_bounds_to_extent.cpp LIBS teca_core teca_data ${teca_test_link} COMMAND test_bounds_to_extent) + +teca_add_test(test_elevation_mask + SOURCES test_elevation_mask.cpp + LIBS teca_core teca_data teca_io teca_alg ${teca_test_link} + COMMAND test_elevation_mask cf + "${TECA_DATA_ROOT}/ECMWF-IFS-HR-SST-present_himalaya\\.nc" + "${TECA_DATA_ROOT}/GTOPO_DEM_025deg\\.nc" + "${TECA_DATA_ROOT}/test_elevation_mask" + FEATURES ${TECA_HAS_NETCDF} + REQ_TECA_DATA) diff --git a/test/apps/CMakeLists.txt b/test/apps/CMakeLists.txt index 8abb36140..607c7aa85 100644 --- a/test/apps/CMakeLists.txt +++ b/test/apps/CMakeLists.txt @@ -109,6 +109,28 @@ teca_add_app_test(test_bayesian_ar_detect_app_packed_data_mpi FEATURES ${TECA_HAS_NETCDF_MPI} ${TECA_HAS_MPI} ${TEST_MPI_THREADS} REQ_TECA_DATA) +teca_add_app_test(test_integrated_water_vapor_app_threads + COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_integrated_water_vapor_app.sh + ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT} -1 + FEATURES ${TECA_HAS_NETCDF} + REQ_TECA_DATA) + +teca_add_app_test(test_integrated_water_vapor_app_mpi + teca_integrated_water_vapor + COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_integrated_water_vapor_app.sh + ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT} -1 + ${MPIEXEC} ${TEST_CORES} + FEATURES ${TECA_HAS_NETCDF_MPI} ${TECA_HAS_MPI} + REQ_TECA_DATA) + +teca_add_app_test(test_integrated_water_vapor_app_mpi_threads + teca_integrated_water_vapor + COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_integrated_water_vapor_app.sh + ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT} -1 + ${MPIEXEC} ${HALF_TEST_CORES} + FEATURES ${TECA_HAS_NETCDF_MPI} ${TECA_HAS_MPI} ${TEST_MPI_THREADS} + REQ_TECA_DATA) + teca_add_app_test(test_integrated_vapor_transport_app_threads COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/test_integrated_vapor_transport_app.sh ${CMAKE_BINARY_DIR}/${BIN_PREFIX} ${TECA_DATA_ROOT} -1 @@ -473,5 +495,8 @@ foreach (app_name ${app_names}) teca_add_app_test(${app_name}_advanced_help ${app_name} COMMAND ${CMAKE_BINARY_DIR}/${BIN_PREFIX}/${app_name} --advanced_help) + teca_add_app_test(${app_name}_no_args ${app_name} + COMMAND ${CMAKE_BINARY_DIR}/${BIN_PREFIX}/${app_name} + WILL_FAIL) endif() endforeach() diff --git a/test/apps/test_bayesian_ar_detect_app.sh b/test/apps/test_bayesian_ar_detect_app.sh index de93bce27..a46d1e881 100755 --- a/test/apps/test_bayesian_ar_detect_app.sh +++ b/test/apps/test_bayesian_ar_detect_app.sh @@ -23,7 +23,7 @@ set -x # run the app ${launcher} ${app_prefix}/teca_bayesian_ar_detect \ --input_regex "${data_root}/ARTMIP_MERRA_2D_2017-05.*\.nc$" \ - --ar_weighted_variables IVT \ + --ar_weighted_variables IVT --segment_ar_probability \ --output_file test_bayesian_ar_detect_app_output_%t%.nc \ --steps_per_file 365 --n_threads ${n_threads} --verbose diff --git a/test/apps/test_bayesian_ar_detect_app_mcf.sh b/test/apps/test_bayesian_ar_detect_app_mcf.sh index c273c363b..8cb80b068 100755 --- a/test/apps/test_bayesian_ar_detect_app_mcf.sh +++ b/test/apps/test_bayesian_ar_detect_app_mcf.sh @@ -24,7 +24,7 @@ set -x ${launcher} ${app_prefix}/teca_bayesian_ar_detect \ --input_file "${app_prefix}/../test/ECMWF-IFS-HR-SST-present.mcf" \ --compute_ivt --wind_u ua --wind_v va --specific_humidity hus \ - --write_ivt --write_ivt_magnitude \ + --segment_ar_probability --write_ivt --write_ivt_magnitude \ --output_file test_bayesian_ar_detect_app_mcf_output_%t%.nc \ --steps_per_file 365 --first_step 8 --last_step 23 \ --n_threads ${n_threads} --verbose diff --git a/test/apps/test_bayesian_ar_detect_app_packed_data.sh b/test/apps/test_bayesian_ar_detect_app_packed_data.sh index 3d1268618..40d131a5e 100755 --- a/test/apps/test_bayesian_ar_detect_app_packed_data.sh +++ b/test/apps/test_bayesian_ar_detect_app_packed_data.sh @@ -24,8 +24,9 @@ set -x ${launcher} ${app_prefix}/teca_bayesian_ar_detect \ --input_regex "${data_root}/ERAinterim_1979-01-0.*\.nc$" \ --x_axis_variable longitude --y_axis_variable latitude --z_axis_variable level \ - --wind_u uwnd --wind_v vwnd --specific_humidity shum --compute_ivt --write_ivt \ - --write_ivt_magnitude --steps_per_file 256 --n_threads ${n_threads} --verbose \ + --wind_u uwnd --wind_v vwnd --specific_humidity shum --segment_ar_probability \ + --compute_ivt --write_ivt --write_ivt_magnitude --steps_per_file 256 \ + --n_threads ${n_threads} --verbose \ --output_file test_bayesian_ar_detect_app_packed_data_output_%t%.nc do_test=1 diff --git a/test/apps/test_integrated_water_vapor_app.sh b/test/apps/test_integrated_water_vapor_app.sh new file mode 100755 index 000000000..06edffd8e --- /dev/null +++ b/test/apps/test_integrated_water_vapor_app.sh @@ -0,0 +1,47 @@ +#!/bin/bash + +if [[ $# < 3 ]] +then + echo "usage: test_integrated_water_vapor_app.sh [app prefix] " \ + "[data root] [num threads] [mpiexec] [num ranks]" + exit -1 +fi + +app_prefix=${1} +data_root=${2} +n_threads=${3} + +if [[ $# -eq 5 ]] +then + mpi_exec=${4} + test_cores=${5} + launcher="${mpi_exec} -n ${test_cores}" +fi + +set -x +set -e + +# run the app +${launcher} ${app_prefix}/teca_integrated_water_vapor \ + --input_file "${app_prefix}/../test/ECMWF-IFS-HR-SST-present.mcf" \ + --specific_humidity hus --n_threads ${n_threads} \ + --output_file test_integrated_water_vapor_app_mcf_output_%t%.nc \ + --verbose --steps_per_file 365 + +do_test=1 +if [[ ${do_test} == 1 ]] +then + # run the diff + ${app_prefix}/teca_cartesian_mesh_diff \ + --reference_dataset ${data_root}/test_integrated_water_vapor_app_mcf_ref'.*\.nc' \ + --test_dataset test_integrated_water_vapor_app_mcf_output'.*\.nc' --arrays IWV \ + --verbose +else + for f in `ls test_integrated_water_vapor_app_mcf_output_*.nc` + do + cp -vd $f ${data_root}/${f/output/ref} + done +fi + +# clean up +rm test_integrated_water_vapor_app_mcf_output*.nc diff --git a/test/test_elevation_mask.cpp b/test/test_elevation_mask.cpp new file mode 100644 index 000000000..ee7b71e43 --- /dev/null +++ b/test/test_elevation_mask.cpp @@ -0,0 +1,171 @@ +#include "teca_config.h" +#include "teca_cf_reader.h" +#include "teca_cartesian_mesh_regrid.h" +#include "teca_cartesian_mesh_source.h" +#include "teca_normalize_coordinates.h" +#include "teca_elevation_mask.h" +#include "teca_indexed_dataset_cache.h" +#include "teca_cf_writer.h" +#include "teca_cf_reader.h" +#include "teca_multi_cf_reader.h" +#include "teca_dataset_diff.h" +#include "teca_file_util.h" +#include "teca_system_util.h" + +#include "teca_index_executive.h" +#include "teca_mpi_manager.h" +#include "teca_system_interface.h" +#include "teca_mpi.h" + +#include +#include +#include +using namespace std; + +int main(int argc, char **argv) +{ + teca_mpi_manager mpi_man(argc, argv); + int rank = mpi_man.get_comm_rank(); + + teca_system_interface::set_stack_trace_on_error(); + teca_system_interface::set_stack_trace_on_mpi_error(); + + if (argc < 5) + { + if (rank == 0) + { + std::cerr << std::endl << "Usage error:" << std::endl + << "test_elevation_mask [reader type (cf,mcf)] [mesh files]" + " [surface elev file] [out file] [first step] [last step]" + << std::endl << std::endl; + } + return -1; + } + std::string reader_type = argv[1]; + std::string mesh_regex = argv[2]; + std::string elev_regex = argv[3]; + std::string baseline = argv[4]; + int first_step = (argc > 5 ? atoi(argv[5]) : 0); + int last_step = (argc > 6 ? atoi(argv[6]) : -1); + + // mesh + p_teca_algorithm mesh_reader; + if (reader_type == "cf") + { + p_teca_cf_reader cfr = teca_cf_reader::New(); + cfr->set_z_axis_variable("plev"); + cfr->set_files_regex(mesh_regex); + mesh_reader = cfr; + } + else if (reader_type == "mcf") + { + p_teca_multi_cf_reader cfr = teca_multi_cf_reader::New(); + cfr->set_z_axis_variable("plev"); + cfr->set_input_file(mesh_regex); + mesh_reader = cfr; + } + + p_teca_normalize_coordinates mesh_coords = teca_normalize_coordinates::New(); + mesh_coords->set_input_connection(mesh_reader->get_output_port()); + + teca_metadata md = mesh_coords->update_metadata(); + + // surface elevations + p_teca_cf_reader elev_cfr = teca_cf_reader::New(); + elev_cfr->set_files_regex(elev_regex); + elev_cfr->set_t_axis_variable(""); + + p_teca_normalize_coordinates elev_coords = teca_normalize_coordinates::New(); + elev_coords->set_input_connection(elev_cfr->get_output_port()); + elev_coords->set_enable_periodic_shift_x(1); + + // regrid the surface elevations onto the mesh coordinates + p_teca_cartesian_mesh_source elev_tgt = teca_cartesian_mesh_source::New(); + elev_tgt->set_spatial_bounds(md, false); + elev_tgt->set_spatial_extents(md, false); + elev_tgt->set_x_axis_variable(md); + elev_tgt->set_y_axis_variable(md); + elev_tgt->set_z_axis_variable(md); + elev_tgt->set_t_axis_variable(md); + elev_tgt->set_t_axis(md); + + p_teca_cartesian_mesh_regrid regrid = teca_cartesian_mesh_regrid::New(); + regrid->set_input_connection(0, elev_tgt->get_output_port()); + regrid->set_input_connection(1, elev_coords->get_output_port()); + + p_teca_algorithm head = regrid; +#ifdef TECA_DEBUG + p_teca_cartesian_mesh_writer vtkwr = teca_cartesian_mesh_writer::New(); + vtkwr->set_input_connection(regrid->get_output_port()); + vtkwr->set_file_name("out_%t%.vtk"); + head = vtkwr; +#endif + + p_teca_indexed_dataset_cache elev_cache = teca_indexed_dataset_cache::New(); + elev_cache->set_input_connection(head->get_output_port()); + elev_cache->set_max_cache_size(1); + + // mask + p_teca_elevation_mask mask = teca_elevation_mask::New(); + mask->set_input_connection(0, mesh_coords->get_output_port()); + mask->set_input_connection(1, elev_cache->get_output_port()); + mask->set_mesh_height_variable("zg"); + mask->set_surface_elevation_variable("z"); + mask->set_mask_variables({"hus_valid","ua_valid","va_valid"}); + + p_teca_index_executive rex = teca_index_executive::New(); + rex->set_verbose(1); + + std::vector arrays({"zg", "hus", "hus_valid", + "ua", "ua_valid", "va", "va_valid"}); + + + // run the test or generate the baseline + bool do_test = true; + teca_system_util::get_environment_variable("TECA_DO_TEST", do_test); + + // find files named by a regex. + std::string tmp_baseline = baseline + ".*\\.nc$"; + std::vector test_files; + std::string regex = teca_file_util::filename(tmp_baseline); + std::string tmp_path = teca_file_util::path(tmp_baseline); + if (do_test && !teca_file_util::locate_files(tmp_path, regex, test_files)) + { + std::cerr << "running the test ..." << std::endl; + + p_teca_cf_reader cfr = teca_cf_reader::New(); + cfr->set_files_regex(baseline); + cfr->set_z_axis_variable("plev"); + + rex->set_arrays(arrays); + + p_teca_dataset_diff diff = teca_dataset_diff::New(); + diff->set_input_connection(0, cfr->get_output_port()); + diff->set_input_connection(1, mask->get_output_port()); + diff->set_executive(rex); + diff->set_verbose(1); + + diff->update(); + } + else + { + std::cerr << "writing the baseline ..." << std::endl; + tmp_baseline = baseline + ".nc"; + + p_teca_cf_writer cmw = teca_cf_writer::New(); + cmw->set_input_connection(mask->get_output_port()); + cmw->set_thread_pool_size(1); + cmw->set_file_name(tmp_baseline); + cmw->set_steps_per_file(10000); + cmw->set_point_arrays(arrays); + + cmw->set_first_step(first_step); + cmw->set_last_step(last_step); + cmw->set_verbose(1); + + cmw->set_executive(rex); + cmw->update(); + } + + return 0; +} diff --git a/test/test_integrated_vapor_transport.cpp b/test/test_integrated_vapor_transport.cpp index b8b63613b..45a22952b 100644 --- a/test/test_integrated_vapor_transport.cpp +++ b/test/test_integrated_vapor_transport.cpp @@ -2,7 +2,9 @@ #include "teca_cartesian_mesh_writer.h" #include "teca_cf_writer.h" #include "teca_valid_value_mask.h" +#include "teca_elevation_mask.h" #include "teca_integrated_vapor_transport.h" +#include "teca_index_executive.h" #include "teca_coordinate_util.h" #include "teca_dataset_capture.h" #include "teca_system_interface.h" @@ -89,8 +91,9 @@ int main(int argc, char **argv) double p_sfc = 92500e-4; double p_top = 5000e-4; double fill_value = 1.0e14; - int write_input = 0; - int write_output = 0; + int vv_mask = (argc > 0 ? atoi(argv[1]) : 1); + int write_input = (argc > 1 ? atoi(argv[2]) : 0); + int write_output = (argc > 2 ? atoi(argv[3]) : 0); // double the z axis, but hit all of the original points. if we set the // integrand to the fill_value where p > p_sfc and apply the valid value @@ -138,25 +141,81 @@ int main(int argc, char **argv) 1, fill_value), v}); - // generate the valid value mask - p_teca_valid_value_mask mask = teca_valid_value_mask::New(); - mask->set_input_connection(mesh->get_output_port()); - mask->set_verbose(0); + p_teca_algorithm head; + if (vv_mask) + { + // generate the valid value mask + std::cerr << "Testing with the valid_value_mask" << std::endl; + + p_teca_valid_value_mask mask = teca_valid_value_mask::New(); + mask->set_input_connection(mesh->get_output_port()); + mask->set_verbose(0); + + head = mask; + } + else + { + // generate the elevation mask + std::cerr << "Testing with the elevation_mask" << std::endl; + + // add generator for mesh height + // let zg = -p + function_of_z zg([](double z) -> double { return -z; }, + 1e6, fill_value); + + mesh->append_field_generator({"zg", + teca_array_attributes(teca_variant_array_code::get(), + teca_array_attributes::point_centering, 0, "m", + "mesh height", "test data where zg = p", + 1, fill_value), + zg}); + + // generate surface elevation + p_teca_cartesian_mesh_source elev = teca_cartesian_mesh_source::New(); + elev->set_whole_extents({0, 2, 0, 2, 0, 0, 0, 0}); + elev->set_bounds({-1.0, 1.0, -1.0, 1.0, p_sfc, p_sfc, 0.0, 0.0}); + elev->set_t_axis_variable(""); + + elev->append_field_generator({"z", + teca_array_attributes(teca_variant_array_code::get(), + teca_array_attributes::point_centering, 0, "m", + "surface elevation", "test data where z = p", + 1, fill_value), + zg}); + + p_teca_elevation_mask mask = teca_elevation_mask::New(); + mask->set_input_connection(0, mesh->get_output_port()); + mask->set_input_connection(1, elev->get_output_port()); + mask->set_mask_variables({"q_valid", "u_valid", "v_valid"}); + mask->set_surface_elevation_variable("z"); + mask->set_mesh_height_variable("zg"); + + head = mask; + } // write the test input dataset if (write_input) { + std::string fn = std::string("test_integrated_vapor_transport_input_") + + std::string(vv_mask ? "vv_mask" : "elev_mask") + std::string("_%t%.nc"); + + p_teca_index_executive exec = teca_index_executive::New(); + p_teca_cf_writer w = teca_cf_writer::New(); - w->set_input_connection(mask->get_output_port()); - w->set_file_name("test_integrated_vapor_transport_input_%t%.nc"); + w->set_input_connection(head->get_output_port()); + w->set_file_name(fn); w->set_thread_pool_size(1); w->set_point_arrays({"q", "u", "v", "q_valid", "u_valid", "v_valid"}); + if (!vv_mask) + w->append_point_array("zg"); + w->set_executive(exec); + w->update(); } // compute IVT p_teca_integrated_vapor_transport ivt = teca_integrated_vapor_transport::New(); - ivt->set_input_connection(mask->get_output_port()); + ivt->set_input_connection(head->get_output_port()); ivt->set_wind_u_variable("u"); ivt->set_wind_v_variable("v"); ivt->set_specific_humidity_variable("q"); @@ -164,9 +223,12 @@ int main(int argc, char **argv) // write the result if (write_output) { + std::string fn = std::string("test_integrated_vapor_transport_output_") + + std::string(vv_mask ? "vv_mask" : "elev_mask") + std::string("_%t%.nc"); + p_teca_cf_writer w = teca_cf_writer::New(); w->set_input_connection(ivt->get_output_port()); - w->set_file_name("test_integrated_vapor_transport_output_%t%.nc"); + w->set_file_name(fn); w->set_thread_pool_size(1); w->set_point_arrays({"ivt_u", "ivt_v"}); w->update(); diff --git a/test/test_integrated_water_vapor.cpp b/test/test_integrated_water_vapor.cpp new file mode 100644 index 000000000..53aa3e145 --- /dev/null +++ b/test/test_integrated_water_vapor.cpp @@ -0,0 +1,239 @@ +#include "teca_cartesian_mesh_source.h" +#include "teca_cartesian_mesh_writer.h" +#include "teca_cf_writer.h" +#include "teca_valid_value_mask.h" +#include "teca_elevation_mask.h" +#include "teca_integrated_water_vapor.h" +#include "teca_index_executive.h" +#include "teca_coordinate_util.h" +#include "teca_dataset_capture.h" +#include "teca_system_interface.h" + + +#include +#include + + +// iwv = - \frac{1}{g} \int_{p_{sfc}}^{p_{top}} q dp +// +// strategy: define integrand q = exp(p) and evaluate numerically using our +// implementation and analytically and compare the results as a check +// +// let q = exp(p) +// +// we can then validate our implementation against the analyitic solutions: +// +// iwv = exp(p) |_{p_{sfc}}^{p_{top}} +// + +template +struct function_of_z +{ + using f_type = std::function; + + function_of_z(const f_type &a_f, num_t max_z, num_t fill_value) : + m_max_z(max_z), m_fill_value(fill_value), m_f(a_f) {} + + p_teca_variant_array operator()(const const_p_teca_variant_array &x, + const const_p_teca_variant_array &y, const const_p_teca_variant_array &z, + double t) + { + (void)t; + + size_t nx = x->size(); + size_t ny = y->size(); + size_t nz = z->size(); + size_t nxy = nx*ny; + + p_teca_variant_array_impl fz = teca_variant_array_impl::New(nx*ny*nz); + num_t *pfz = fz->get(); + + TEMPLATE_DISPATCH(teca_variant_array_impl, + fz.get(), + const NT *pz = dynamic_cast(z.get())->get(); + for (size_t k = 0; k < nz; ++k) + { + for (size_t j = 0; j < ny; ++j) + { + for (size_t i = 0; i < nx; ++i) + { + num_t z = pz[k]; + pfz[k*nxy + j*nx + i] = z > m_max_z ? m_fill_value : this->m_f(z); + } + } + } + ) + + return fz; + } + + num_t m_max_z; + num_t m_fill_value; + f_type m_f; +}; + +// an actual sequence of pressure levels from the HighResMIP data +// plev = 92500, 85000, 70000, 60000, 50000, 25000, 5000 ; + +int main(int argc, char **argv) +{ + (void)argc; + (void)argv; + + teca_system_interface::set_stack_trace_on_error(); + + unsigned long i1 = 1024; + double p_sfc = 92500e-4; + double p_top = 5000e-4; + double fill_value = 1.0e14; + int vv_mask = (argc > 0 ? atoi(argv[1]) : 1); + int write_input = (argc > 1 ? atoi(argv[2]) : 0); + int write_output = (argc > 2 ? atoi(argv[3]) : 0); + + // double the z axis, but hit all of the original points. if we set the + // integrand to the fill_value where p > p_sfc and apply the valid value + // mask then the integral should have the value as if integrated from p_sfc + // to p_top. This lets us verify that the integrator works correctly in + // the presence of missing values + double p_sfc_2 = 2.0*p_sfc - p_top; + unsigned long j1 = 2*i1; + + p_teca_cartesian_mesh_source mesh = teca_cartesian_mesh_source::New(); + mesh->set_whole_extents({0, 2, 0, 2, 0, j1, 0, 0}); + mesh->set_bounds({-1.0, 1.0, -1.0, 1.0, p_sfc_2, p_top, 0.0, 0.0}); + mesh->set_calendar("standard", "days since 2020-09-30 00:00:00"); + + // let q = exp(p) + function_of_z q([](double p) -> double { return exp(p); }, + p_sfc, fill_value); + + mesh->append_field_generator({"q", + teca_array_attributes(teca_variant_array_code::get(), + teca_array_attributes::point_centering, 0, "g kg^-1", + "specific humidty", "test data where q = sin(p)", + 1, fill_value), + q}); + + p_teca_algorithm head; + if (vv_mask) + { + // generate the valid value mask + std::cerr << "Testing with the valid_value_mask" << std::endl; + + p_teca_valid_value_mask mask = teca_valid_value_mask::New(); + mask->set_input_connection(mesh->get_output_port()); + mask->set_verbose(0); + + head = mask; + } + else + { + // generate the elevation mask + std::cerr << "Testing with the elevation_mask" << std::endl; + + // add generator for mesh height + // let zg = -p + function_of_z zg([](double z) -> double { return -z; }, + 1e6, fill_value); + + mesh->append_field_generator({"zg", + teca_array_attributes(teca_variant_array_code::get(), + teca_array_attributes::point_centering, 0, "m", + "mesh height", "test data where zg = p", + 1, fill_value), + zg}); + + // generate surface elevation + p_teca_cartesian_mesh_source elev = teca_cartesian_mesh_source::New(); + elev->set_whole_extents({0, 2, 0, 2, 0, 0, 0, 0}); + elev->set_bounds({-1.0, 1.0, -1.0, 1.0, p_sfc, p_sfc, 0.0, 0.0}); + elev->set_t_axis_variable(""); + + elev->append_field_generator({"z", + teca_array_attributes(teca_variant_array_code::get(), + teca_array_attributes::point_centering, 0, "m", + "surface elevation", "test data where z = p", + 1, fill_value), + zg}); + + p_teca_elevation_mask mask = teca_elevation_mask::New(); + mask->set_input_connection(0, mesh->get_output_port()); + mask->set_input_connection(1, elev->get_output_port()); + mask->set_mask_variables({"q_valid"}); + mask->set_surface_elevation_variable("z"); + mask->set_mesh_height_variable("zg"); + + head = mask; + } + + // write the test input dataset + if (write_input) + { + std::string fn = std::string("test_integrated_water_vapor_input_") + + std::string(vv_mask ? "vv_mask" : "elev_mask") + std::string("_%t%.nc"); + + p_teca_index_executive exec = teca_index_executive::New(); + + p_teca_cf_writer w = teca_cf_writer::New(); + w->set_input_connection(head->get_output_port()); + w->set_file_name(fn); + w->set_thread_pool_size(1); + w->set_point_arrays({"q", "q_valid"}); + if (!vv_mask) + w->append_point_array("zg"); + w->set_executive(exec); + + w->update(); + } + + // compute IWV + p_teca_integrated_water_vapor iwv = teca_integrated_water_vapor::New(); + iwv->set_input_connection(head->get_output_port()); + iwv->set_specific_humidity_variable("q"); + iwv->set_iwv_variable("iwv"); + + // write the result + if (write_output) + { + std::string fn = std::string("test_integrated_water_vapor_output_") + + std::string(vv_mask ? "vv_mask" : "elev_mask") + std::string("_%t%.nc"); + + p_teca_cf_writer w = teca_cf_writer::New(); + w->set_input_connection(iwv->get_output_port()); + w->set_file_name(fn); + w->set_thread_pool_size(1); + w->set_point_arrays({"iwv"}); + w->update(); + } + + // capture the result + p_teca_dataset_capture dsc = teca_dataset_capture::New(); + dsc->set_input_connection(iwv->get_output_port()); + dsc->update(); + + const_p_teca_cartesian_mesh m = + std::dynamic_pointer_cast(dsc->get_dataset()); + + const_p_teca_double_array iwva = + std::dynamic_pointer_cast + (m->get_point_arrays()->get("iwv")); + + double test_iwv = iwva->get(0); + + // calculate the analytic solution + double g = 9.80665; + double base_iwv = -(1./g) * (exp(p_top) - exp(p_sfc)); + + // display + std::cerr << "base_iwv = " << base_iwv << std::endl; + std::cerr << "test_iwv = " << test_iwv << std::endl << std::endl; + + // compare against the analytic solution + if (!teca_coordinate_util::equal(base_iwv, test_iwv, 1e-4)) + { + TECA_ERROR("base_iwv = " << base_iwv << " test_iwv = " << test_iwv) + return -1; + } + + return 0; +}