Environment.h

/* Copyright (c) Queen's Printer for Ontario, 2020. */
/* Copyright (c) His Majesty the King in Right of Canada as represented by the Minister of Natural Resources, 2021-2025. */
 
/* SPDX-License-Identifier: AGPL-3.0-or-later */
 
#pragma once
#include "stdafx.h"
#include "Cell.h"
#include "ConstantGrid.h"
#include "Event.h"
#include "FuelType.h"
#include "GridMap.h"
#include "IntensityMap.h"
#include "Point.h"
#include "Settings.h"
 
namespace fs::topo
{
using FuelGrid = data::ConstantGrid<const fuel::FuelType*, FuelSize>;
using ElevationGrid = data::ConstantGrid<ElevationSize>;
using CellGrid = data::ConstantGrid<Cell, Topo>;
class Environment
{
public:
  [[nodiscard]] static Environment loadEnvironment(
    const string dir_out,
    const string& path,
    const Point& point,
    const string& perimeter,
    int year);
  [[nodiscard]] static Environment load(const string dir_out,
                                        const Point& point,
                                        const string& in_fuel,
                                        const string& in_elevation);
  ~Environment();
  [[nodiscard]] unique_ptr<Coordinates> findCoordinates(
    const Point& point,
    bool flipped) const;
  Environment(Environment&& rhs) noexcept = default;
  Environment(const Environment& rhs) noexcept = default;
  Environment& operator=(Environment&& rhs) noexcept = default;
  Environment& operator=(const Environment& rhs) noexcept = default;
  [[nodiscard]] constexpr const string& proj4() const
  {
    return cells_->proj4();
  }
  [[nodiscard]] constexpr Idx rows() const
  {
    return cells_->rows();
  }
  [[nodiscard]] constexpr Idx columns() const
  {
    return cells_->columns();
  }
  [[nodiscard]] constexpr MathSize cellSize() const
  {
    return cells_->cellSize();
  }
  [[nodiscard]] constexpr ElevationSize elevation() const
  {
    return elevation_;
  }
  [[nodiscard]]
#ifdef NDEBUG
  constexpr
#endif
    Cell
    cell(const Idx row, const Idx column) const
  {
    return cells_->at(Location(row, column));
  }
  template <class P>
  [[nodiscard]] constexpr Cell cell(const Position<P>& position) const
  {
    return cells_->at(position);
  }
  //  [[nodiscard]] constexpr Cell cell(const HashSize hash_size) const
  //  {
  //    return cells_->at(hash_size);
  //  }
  [[nodiscard]]
#ifdef NDEBUG
  constexpr
#endif
    Cell
    offset(const sim::Event& event,
           const Idx row,
           const Idx column) const
  {
    const auto& p = event.cell();
    // return cell(p.hash() + column + static_cast<HashSize>(MAX_COLUMNS) * row);
    return cell(Location(p.row() + row, p.column() + column));
  }
  [[nodiscard]] sim::ProbabilityMap* makeProbabilityMap(DurationSize time,
                                                        DurationSize start_time,
                                                        int min_value,
                                                        int low_max,
                                                        int med_max,
                                                        int max_value) const;
  template <class Other>
  [[nodiscard]] unique_ptr<data::GridMap<Other>> makeMap(const Other nodata) const
  {
    return make_unique<data::GridMap<Other>>(*cells_, nodata);
  }
  [[nodiscard]] unique_ptr<sim::BurnedData> makeBurnedData() const
  {
    auto result = make_unique<sim::BurnedData>(*not_burnable_);
    return result;
  }
  void resetBurnedData(sim::BurnedData* data) const noexcept
  {
    // setting to {} probably makes a bitset of the same size on the stack?
    // *data = {};
    *data = *not_burnable_;
  }
protected:
  [[nodiscard]] static CellGrid* makeCells(
    const FuelGrid& fuel,
    const ElevationGrid& elevation)
  {
    logging::check_equal(fuel.yllcorner(), elevation.yllcorner(), "yllcorner");
    static Cell nodata{};
    auto values = vector<Cell>{fuel.data.size()};
    vector<HashSize> hashes{};
    //    for (HashSize h = 0; h < static_cast<size_t>(MAX_ROWS) * MAX_COLUMNS; ++h)
    //    {
    //      hashes.emplace_back(h);
    //    }
    for (Idx r = 0; r < fuel.rows(); ++r)
    {
      for (Idx c = 0; c < fuel.columns(); ++c)
      {
        const auto h = hashes.emplace_back(Location(r, c).hash());
        //        hashes.emplace_back(Location(r, c).hash());
        //      }
        //    }
        //    std::for_each(
        //      std::execution::par_unseq,
        //      hashes.begin(),
        //      hashes.end(),
        //      [&fuel, &values, &elevation](auto&& h)
        //      {
        const topo::Location loc{r, c, h};
        //        const topo::Location loc{static_cast<Idx>(h / MAX_COLUMNS),
        //                                 static_cast<Idx>(h % MAX_COLUMNS),
        //                                 h};
        //        const auto r = loc.row();
        //        const auto c = loc.column();
        //        const auto f = fuel::FuelType::safeCode(fuel.at(h));
        if (r >= 0 && r < fuel.rows() && c >= 0 && c < fuel.columns())
        {
          // NOTE: this needs to translate to internal codes?
          const auto f = fuel::FuelType::safeCode(fuel.at(loc));
          auto s = static_cast<SlopeSize>(INVALID_SLOPE);
          auto a = static_cast<AspectSize>(INVALID_ASPECT);
          // HACK: don't calculate for outside box of cells
          if (r > 0 && r < fuel.rows() - 1 && c > 0 && c < fuel.columns() - 1)
          {
            MathSize dem[9];
            bool valid = true;
            for (int i = -1; i < 2; ++i)
            {
              for (int j = -1; j < 2; ++j)
              {
                // grid is (0, 0) at bottom left, but want [0] in array to be NW corner
                auto actual_row = static_cast<Idx>(r - i);
                auto actual_column = static_cast<Idx>(c + j);
                auto cur_loc = Location{actual_row, actual_column};
                //                auto cur_h = cur_loc.hash();
                //              dem[3 * (i + 1) + (j + 1)] = 1.0 * elevation.at(cur_h);
                const auto v = elevation.at(cur_loc);
                // can't calculate slope & aspect if any surrounding cell is nodata
                if (elevation.nodataValue() == v)
                {
                  valid = false;
                  break;
                }
                dem[3 * (i + 1) + (j + 1)] = 1.0 * v;
              }
              if (!valid)
              {
                break;
              }
            }
            if (valid)
            {
              // Horn's algorithm
              const MathSize dx = ((dem[2] + dem[5] + dem[5] + dem[8])
                                   - (dem[0] + dem[3] + dem[3] + dem[6]))
                                / elevation.cellSize();
              const MathSize dy = ((dem[6] + dem[7] + dem[7] + dem[8])
                                   - (dem[0] + dem[1] + dem[1] + dem[2]))
                                / elevation.cellSize();
              const MathSize key = (dx * dx + dy * dy);
              auto slope_pct = static_cast<float>(100 * (sqrt(key) / 8.0));
              s = min(static_cast<SlopeSize>(MAX_SLOPE_FOR_DISTANCE), static_cast<SlopeSize>(round(slope_pct)));
              static_assert(std::numeric_limits<SlopeSize>::max() >= MAX_SLOPE_FOR_DISTANCE);
              MathSize aspect_azimuth = 0.0;
 
              if (s > 0 && (dx != 0 || dy != 0))
              {
                aspect_azimuth = atan2(dy, -dx) * M_RADIANS_TO_DEGREES;
                // NOTE: need to change this out of 'math' direction into 'real' direction (i.e. N is 0, not E)
                aspect_azimuth = (aspect_azimuth > 90.0) ? (450.0 - aspect_azimuth) : (90.0 - aspect_azimuth);
                if (aspect_azimuth == 360.0)
                {
                  aspect_azimuth = 0.0;
                }
              }
 
              a = static_cast<AspectSize>(round(aspect_azimuth));
            }
          }
          const auto cell = Cell{h, s, a, f};
          values.at(h) = cell;
#ifdef DEBUG_GRIDS
#ifndef VLD_RPTHOOK_INSTALL
          logging::check_equal(cell.row(), r, "Cell row");
          logging::check_equal(cell.column(), c, "Cell column");
          const auto v = values.at(h);
          logging::check_equal(v.row(), r, "Row");
          logging::check_equal(v.column(), c, "Column");
          if (!(INVALID_SLOPE == cell.slope() || INVALID_ASPECT == cell.aspect() || INVALID_FUEL_CODE == cell.fuelCode()))
          {
            logging::check_equal(cell.slope(), s, "Cell slope");
            logging::check_equal(v.slope(), s, "Slope");
            if (0 != s)
            {
              logging::check_equal(cell.aspect(), a, "Cell aspect");
              logging::check_equal(v.aspect(), a, "Aspect");
            }
            else
            {
              logging::check_equal(cell.aspect(), a, "Cell aspect when slope is 0");
              logging::check_equal(v.aspect(), 0, "Aspect when slope is 0");
            }
            logging::check_equal(v.fuelCode(), f, "Fuel");
            logging::check_equal(cell.fuelCode(), f, "Cell fuel");
          }
          else
          {
            logging::check_equal(cell.slope(), INVALID_SLOPE, "Invalid Cell slope");
            logging::check_equal(cell.aspect(), INVALID_ASPECT, "Invalid Cell aspect");
            logging::check_equal(cell.fuelCode(), INVALID_FUEL_CODE, "Invalid Cell fuel");
            logging::check_equal(v.slope(), INVALID_SLOPE, "Invalid slope");
            logging::check_equal(v.aspect(), INVALID_ASPECT, "Invalid aspect");
            logging::check_equal(v.fuelCode(), INVALID_FUEL_CODE, "Invalid fuel");
          }
#endif
#endif
        }
        //      });
      }
    }
    return new topo::CellGrid(
      fuel.cellSize(),
      fuel.rows(),
      fuel.columns(),
      nodata.fullHash(),
      nodata,
      fuel.xllcorner(),
      fuel.yllcorner(),
      fuel.xurcorner(),
      fuel.yurcorner(),
      string(fuel.proj4()),
      std::move(values));
  }
  shared_ptr<sim::BurnedData> initializeNotBurnable(const CellGrid& cells) const
  {
    // shared_ptr<sim::BurnedData> result{};
    //     std::fill(not_burnable_.begin(), not_burnable_.end(), false);
    //  make a template we can copy to reset things
    auto result = make_shared<sim::BurnedData>();
    for (Idx r = 0; r < rows(); ++r)
    {
      for (Idx c = 0; c < columns(); ++c)
      {
        const Location location(r, c);
        // HACK: just mark outside edge as unburnable so we never need to check
        bool is_outer = 0 == r
                     || 0 == c
                     || (rows() - 1) == r
                     || (columns() - 1) == c;
        // (*result)[location.hash()] = (nullptr == fuel::fuel_by_code(cells.at(location).fuelCode()));
        (*result)[location.hash()] = is_outer || fuel::is_null_fuel(cells.at(location));
        //        if (fuel::is_null_fuel(cell(location)))
        //        {
        //          not_burnable_[location.hash()] = true;
        //        }
      }
    }
    return result;
  }
  Environment(const string dir_out,
              CellGrid* cells,
              const ElevationSize elevation) noexcept
    : dir_out_(dir_out),
      cells_(cells),
      not_burnable_(initializeNotBurnable(*cells)),
      elevation_(elevation)
  {
    //    try
    //    {
    //      initializeNotBurnable();
    //    }
    //    catch (...)
    //    {
    //      std::terminate();
    //    }
  }
  Environment(
    const string dir_out,
    const FuelGrid& fuel,
    const ElevationGrid& elevation,
    const Point& point)
    : Environment(dir_out,
                  makeCells(fuel,
                            elevation),
                  elevation.at(Location(*elevation.findCoordinates(point, false).get())))
  {
    // take elevation at point so that if max grid size changes elevation doesn't
    logging::note("Start elevation is %d", elevation_);
    if (sim::Settings::saveSimulationArea())
    {
      logging::debug("Saving simulation area");
      const auto lookup = sim::Settings::fuelLookup();
      auto convert_to_slope =
        [&elevation](
          const Cell& v) -> SlopeSize {
        return v.slope();
      };
      auto convert_to_aspect =
        [&elevation](
          const Cell& v) -> AspectSize {
        return v.aspect();
      };
      auto convert_to_area =
        [&elevation](
          const ElevationSize v) -> ElevationSize {
        // need to still be nodata if it was
        return (v == elevation.nodataValue()) ? v : 3;
      };
      auto convert_to_fuelcode =
        [&lookup](
          const fuel::FuelType* const value) -> FuelSize {
        return lookup.fuelToCode(value);
      };
      fuel.saveToFile<FuelSize>(
        dir_out_,
        "fuel",
        convert_to_fuelcode);
      elevation.saveToFile(
        dir_out_,
        "dem");
      // save slope & aspect grids
      cells_->saveToFile<SlopeSize>(
        dir_out_,
        "slope",
        convert_to_slope,
        static_cast<SlopeSize>(INVALID_SLOPE));
      cells_->saveToFile<AspectSize>(
        dir_out_,
        "aspect",
        convert_to_aspect,
        static_cast<AspectSize>(INVALID_ASPECT));
      // HACK: make a grid with "3" as the value so if we merge max with it it'll cover up anything else
      elevation.saveToFile<ElevationSize>(
        dir_out_,
        "simulation_area",
        convert_to_area);
      logging::debug("Done saving fuel grid");
    }
  }
private:
  const string dir_out_;
  CellGrid* cells_;
  shared_ptr<const sim::BurnedData> not_burnable_;
  ElevationSize elevation_;
};
}