hydrology.hpp

Go to the documentation of this file.

/* Copyright (c) 2023 Otto Link. Distributed under the terms of the GNU General
   Public License. The full license is in the file LICENSE, distributed with
   this software. */
 
#pragma once
#include <functional>
#include <limits>
 
#include "highmap/array.hpp"
#include "highmap/geometry/path.hpp"
 
namespace hmap
{
 
enum FlowDirectionMethod : int
{
  FDM_D8,
  FDM_PRIORITY_FLOOD
};
 
class DrainageBasinCellBased
{
public:
  std::vector<std::vector<glm::ivec2>> upstream_traversal;
  Mat<glm::ivec2>                      next = Mat<glm::ivec2>({0, 0});
  glm::ivec2                           null_cell = glm::ivec2(-1, -1);
 
  size_t                               get_basins_number() const;
  std::vector<std::vector<glm::ivec2>> get_main_channels();
  std::vector<glm::ivec2>              get_outlets() const;
  std::vector<glm::ivec2>              get_ridges();
  std::vector<std::vector<glm::ivec2>> get_ridges_neighbors();
 
  void generate_traversal(
      const Array                   &z,
      FlowDirectionMethod            fd_method = FlowDirectionMethod::FDM_D8,
      bool                           remove_lakes = true,
      const std::vector<glm::ivec2> &outlets = {});
 
  void generate_traversal_d8(const Array                   &z,
                             bool                           remove_lakes = true,
                             const std::vector<glm::ivec2> &outlets = {});
  void generate_traversal_priority_flood(
      const Array                   &z,
      const std::vector<glm::ivec2> &outlets = {});
 
  void accumulate(Array &acc) const;
  void traverse_downstream(std::function<void(int, int, int, int, int)> op);
  void traverse_downstream(std::function<void(int, int, int)> op);
  void traverse_upstream(std::function<void(int, int, int, int, int)> op);
  void traverse_upstream(std::function<void(int, int, int)> op);
 
private:
  void remove_lakes_d8(const Array &z,
                       float        dz_weight = 1.f,
                       float        dz_downstream_cost_ratio = 0.f);
  void update_traversal();
};
 
Array basin_id(const Array        &z,
               FlowDirectionMethod fd_method = FlowDirectionMethod::FDM_D8,
               bool                remove_lakes = true);
 
Array d8_compute_ndip(const Array &d8);
 
void find_flow_apex(const Array &z, std::vector<int> &is, std::vector<int> &js);
 
void find_flow_sinks(const Array      &z,
                     std::vector<int> &is,
                     std::vector<int> &js);
 
std::vector<glm::ivec2> find_flow_sinks(const Array &z);
 
Array flooding_uniform_level(const Array &z, float zref);
 
Array flooding_from_boundaries(const Array &z,
                               float        zref,
                               bool         from_east = true,
                               bool         from_west = true,
                               bool         from_north = true,
                               bool         from_south = true);
 
Array flooding_from_point(const Array &z,
                          int          i,
                          int          j,
                          float depth_min = std::numeric_limits<float>::max());
 
Array flooding_from_point(const Array            &z,
                          const std::vector<int> &i,
                          const std::vector<int> &j,
                          float depth_min = std::numeric_limits<float>::max());
 
Array flooding_lake_system(const Array &z, float surface_threshold = 0);
 
Array flow_accumulation_d8(const Array &z);
 
Array flow_accumulation_dinf(const Array &z, float talus_ref);
 
Array flow_accumulation_dinf_perturbed(const Array &z,
                                       float        talus_ref,
                                       int          nsamples,
                                       glm::vec2    kw,
                                       uint         seed,
                                       float        amp,
                                       const Array *p_perturb_scaling = nullptr,
                                       glm::vec4 bbox = {0.f, 1.f, 0.f, 1.f});
 
Array flow_direction_d8(const Array &z);
 
std::vector<Array> flow_direction_dinf(const Array &z, float talus_ref);
 
std::vector<float> flow_direction_dinf_flat(const Array &z, float talus_ref);
 
Array flow_direction_dinf_angle(const Array &z, float talus_ref);
 
Array flow_fixing(const Array &z,
                  float        riverbed_talus = 0.f,
                  int          iterations = 5,
                  int          prefilter_ir = 8,
                  bool         carve_riverbed = true,
                  bool         smooth_river_bottom = true,
                  float        talus_riverbank = 0.01f, // 4 / shape.X
                  uint         seed = 0,
                  float        riverbank_noise_ratio = 0.f,
                  float        merging_distance = 8.f, // pixels
                  const Array *p_noise_x = nullptr,
                  const Array *p_noise_y = nullptr);
 
Path flow_stream(const Array     &z,
                 const glm::ivec2 ij_start,
                 const float      elevation_ratio = 0.5f,
                 const float      distance_exponent = 2.f,
                 const float      upward_penalization = 100.f);
 
Array generate_riverbed(const Path &path,
                        glm::ivec2  shape,
                        glm::vec4   bbox = {0.f, 1.f, 0.f, 1.f},
                        bool        bezier_smoothing = false,
                        float       depth_start = 0.01f,
                        float       depth_end = 1.f,
                        float       slope_start = 64.f,
                        float       slope_end = 32.f,
                        float       shape_exponent_start = 1.f,
                        float       shape_exponent_end = 10.f,
                        float       k_smoothing = 0.5f,
                        int         post_filter_ir = 0,
                        Array      *p_noise_x = nullptr,
                        Array      *p_noise_y = nullptr,
                        Array      *p_noise_r = nullptr);
 
Array merge_water_depths(const Array &depth1,
                         const Array &depth2,
                         float        k_smooth = 0.f);
 
Array snow_melting_map(const Array &z,
                       float        melt_start_elevation = 0.f,
                       float        melt_end_elevation = 0.5f,
                       float        elevation_strength = 1.f,
                       float        elevation_exp = 1.f,
                       float        sun_azimuth = 0.f,
                       float        sun_zenith = 60.f,
                       float        aspect_strength = 0.f,
                       float        slope_exp = 1.f,
                       float        slope_strength = 0.f);
 
Array water_depth_from_mask(const Array &z,
                            const Array &mask,
                            float        mask_threshold = 0.f,
                            int          iterations_max = 10000,
                            float        tolerance = 1e-2f,
                            float        omega = 1.8f);
 
void water_depth_dry_out(Array       &water_depth,
                         float        dry_out_ratio = 0.5f,
                         const Array *p_mask = nullptr,
                         float depth_max = std::numeric_limits<float>::max());
 
Array water_depth_increase(const Array &water_depth,
                           const Array &z,
                           float        additional_depth);
 
Array water_mask(const Array &water_depth);
 
Array water_mask(const Array &water_depth,
                 const Array &z,
                 float        additional_depth);
 
} // namespace hmap
 
namespace hmap::gpu
{
 
Array flow_accumulation_from_velocity_field(const Array &u,
                                            const Array &v,
                                            int          iterations);
 
Array flow_direction_d8(const Array &z);
 
Array flow_simulation(const Array &z,
                      float        water_depth,
                      const Array &depth_map,
                      int          iterations,
                      float        dt = 0.5f,
                      bool         flux_diffusion = true,
                      float        flux_diffusion_strength = 0.01f,
                      float        dry_out_ratio = 0.f);
 
Array flow_simulation_viscous(const Array &z,
                              float        water_depth,
                              const Array &depth_map,
                              int          iterations,
                              float        dt = 0.5f,
                              float        dry_out_ratio = 0.f,
                              float        viscosity = 1.f,
                              float        power = 2.5f);
 
Array generate_riverbed(const Path &path,
                        glm::ivec2  shape,
                        glm::vec4   bbox = {0.f, 1.f, 0.f, 1.f},
                        bool        bezier_smoothing = false,
                        float       depth_start = 0.01f,
                        float       depth_end = 1.f,
                        float       slope_start = 64.f,
                        float       slope_end = 32.f,
                        float       shape_exponent_start = 1.f,
                        float       shape_exponent_end = 10.f,
                        float       k_smoothing = 0.5f,
                        int         post_filter_ir = 0,
                        Array      *p_noise_x = nullptr,
                        Array      *p_noise_y = nullptr,
                        Array      *p_noise_r = nullptr);
 
Array snow_simulation(const Array &z,
                      float        snow_depth,
                      const Array &fall_map,
                      const Array &melting_map,
                      const Array &talus,
                      int          iterations,
                      float        dt = 0.5f,
                      float        fall_iterations_ratio = 1.f,
                      float        k_snow = 0.5f,
                      float        k_visc = 0.1f,
                      float        k_melt_factor = 0.8f,
                      float        k_depth_ratio = 1.f,
                      float        k_depth_slope_ratio = 1.f,
                      bool         post_filter = true,
                      float        thermal_talus_ratio = 0.2f);
 
void water_depth_filter(Array &depth, const Array &z, int ir);
 
} // namespace hmap::gpu