Merge pull request #8 from CWFMF/wip-validation

improve spread algorithm

.gitignore

@@ -15,4 +15,6 @@ pyfire/venv
 **/logs
 **/build
 **/private
-
+**/output*
+tbd/validate
+tbd/dir_test

cmp_grid.sh

@@ -0,0 +1,11 @@
+#!/bin/bash
+
+fire="11N_54669"
+# file="probability_215_2023-08-03.tif"
+# file="007_000001_215_arrival.tif"
+file="020_000030_215_.tif"
+
+docker run -it -v /mnt/data/data/tmp/:/home/out \
+    ghcr.io/osgeo/gdal:ubuntu-full-latest gdalcompare.py -skip_binary \
+    /home/out/${fire}/${file} \
+    /home/out/${fire}.orig/${file}

pts_to_shp.py

@@ -0,0 +1,164 @@
+import math
+import os
+import re
+import shutil
+
+import geopandas as gpd
+import pandas as pd
+
+CRS_WGS84 = 4326
+
+DIR = "../data/test_output"
+DIR_SIMS = f"{DIR}/sims"
+DIR_OUT = f"{DIR}/gis"
+
+os.makedirs(DIR_OUT, exist_ok=True)
+
+
+def to_gdf(df, crs=CRS_WGS84):
+    geometry = df["geometry"] if "geometry" in df else gpd.points_from_xy(df["x"], df["y"], crs=crs)
+    return gpd.GeoDataFrame(df, crs=crs, geometry=geometry)
+
+
+STAGES = {"N": 0, "S": 1, "C": 2}
+
+
+def find_groups(s):
+    g = re.fullmatch("(\d*)(.)(\d)*", s).groups()
+    return [int(g[0]), STAGES[g[1]], int(g[2])]
+
+
+def find_id(s):
+    scenario, stage, step = find_groups(s)
+    return step * 10 + stage
+
+
+def to_shp(file, force=False, dir_to=DIR_OUT):
+    d = os.path.dirname(file)
+    name = os.path.basename(d)
+    save_to = f"{dir_to}/{name}.shp"
+    if not (force or not os.path.exists(save_to)):
+        return
+    df = pd.read_csv(file)
+    print(f"Converting {len(df)} points from {file}")
+    # df[["scenario", "stage", "step"]] = list(df["step_id"].apply(find_groups))
+    df["id"] = df["step_id"].apply(find_id)
+    gdf = to_gdf(df, crs=None)
+    # gdf["id"] = gdf.apply(lambda x: x["step"] * 10 + x["stage"], axis=1)
+    # del gdf["scenario"]
+    del gdf["step_id"]
+    del gdf["x"]
+    del gdf["y"]
+    gdf.to_file(
+        save_to,
+        # schema=schema,
+        layer=name,
+    )
+
+
+def to_gpkg(file, force=False, dir_to=DIR_OUT):
+    d = os.path.dirname(file)
+    name = os.path.basename(d)
+    save_to = f"{dir_to}/{name}.gpkg"
+    if not (force or not os.path.exists(save_to)):
+        return
+    df = pd.read_csv(file)
+    print(f"Converting {len(df)} points from {file}")
+    df[["scenario", "stage", "step"]] = list(df["step_id"].apply(find_groups))
+    gdf = to_gdf(df, crs=None)
+    del gdf["step_id"]
+    del gdf["x"]
+    del gdf["y"]
+    gdf.to_file(
+        save_to,
+        driver="GPKG",
+        layer=name,
+    )
+    shutil.make_archive(save_to, "zip", dir_to, os.path.basename(save_to))
+
+
+def to_parquet(file, force=False, dir_to=DIR_OUT):
+    d = os.path.dirname(file)
+    name = os.path.basename(d)
+    save_to = f"{dir_to}/{name}.parquet"
+    if not (force or not os.path.exists(save_to)):
+        df = gpd.read_parquet(save_to)
+        print(f"Already converted {len(df)} points from {name}")
+        return
+    df = pd.read_csv(file)
+    print(f"Converting {len(df)} points from {name}")
+    df[["scenario", "stage", "step"]] = list(df["step_id"].apply(find_groups))
+    # df = df.loc[(df["step"] < 2) | (df["step"] == df["step"].max())]
+    df = df.loc[(df["step"] < 1)]
+    gdf = to_gdf(df, crs=None)
+    del gdf["step_id"]
+    del gdf["x"]
+    del gdf["y"]
+    gdf.to_parquet(save_to, index=False)
+
+
+angles = []
+
+
+def add_offsets_calc_ros(theta):
+    d = math.degrees(theta)
+    print(d)
+    angles.append(d)
+    return True
+
+
+def real_angle(l_b, rads):
+    # return rads
+    rx = 1.0
+    ry = l_b
+    x = math.cos(rads)
+    y = math.sin(rads)
+    aa = rx * rx
+    bb = ry * ry
+    t = aa * bb / ((x * x * bb) + (y * y * aa))
+    x *= t
+    y *= t
+    y *= rx / ry
+    return math.atan2(y, x)
+
+
+# angles = []
+# # head_ros
+# l_b = 5.0
+# theta = 0
+# cur_x = 0
+# step_x = 0.1
+# added = True
+# while added and cur_x < 1:
+#     theta = abs(math.acos(cur_x) - math.radians(90))
+#     added = add_offsets_calc_ros(real_angle(l_b, theta))
+#     # added = add_offsets_calc_ros(theta)
+#     cur_x += step_x
+
+# if added:
+#     # added = add_offsets(math.radians(90), flank_ros * math.sqrt(a_sq_sub_c_sq) / a)
+#     added = add_offsets_calc_ros(math.radians(90))
+#     cur_x -= step_x
+
+# while added and cur_x > 0:
+#     # theta = math.acos(1.0 - cur_x)
+#     theta = math.radians(90) + math.acos(cur_x)
+#     added = add_offsets_calc_ros(real_angle(l_b, theta))
+#     # added = add_offsets_calc_ros(theta)
+#     cur_x -= step_x
+
+# if added:
+#     add_offsets_calc_ros(math.radians(180))
+
+# angles = [180 - x for x in angles]
+
+
+to_convert = []
+for root, dirs, files in os.walk(DIR_SIMS):
+    for f in files:
+        if "_points.txt" in f:
+            to_convert.append(os.path.join(root, f))
+
+to_convert = sorted(to_convert)
+for f in to_convert:
+    to_parquet(f)

source.r

@@ -0,0 +1,140 @@
+setwd("//wsl$/Ubuntu-20.04/mnt/data/FireSTARR")
+library(data.table)
+# library(raster)
+library(terra)
+library(lattice)
+library(rasterVis)
+
+DIR <- "../data/test_condense16"
+# DIR_OUT <- "."
+DIR_OUT <- Sys.getenv("HOME")
+FILE_OUT <- file.path(DIR_OUT, paste0(basename(DIR), ".pdf"))
+# colours <- list(
+# "1"=list(r=255, g=0, b=0),
+# "2"=list(r=255, g=255, b=255),
+# "4"=list(r=247, g=250, b=21),
+# "8"=list(r=0, g=0, b=0),
+# #"16"=list(r=185, g=122, b=87),
+# #"32"=list(r=34, g=177, b=76),
+# #"64"=list(r=163, g=73, b=164),
+# #"128"=list(r=0, g=162, b=232)
+# "16"=list(r=0, g=255, b=0),
+# "32"=list(r=0, g=0, b=255),
+# "64"=list(r=200, g=200, b=0),
+# "128"=list(r=200, g=0, b=200)
+# )
+# colours <- list(
+#   "1"=list(r=255, g=0, b=0),
+#   "2"=list(r=0, g=255, b=0),
+#   "4"=list(r=0, g=0, b=255),
+#   "8"=list(r=255, g=255, b=255),
+#   "16"=list(r=255, g=255, b=0),
+#   "32"=list(r=0, g=255, b=255),
+#   "64"=list(r=255, g=0, b=255),
+#   "128"=list(r=0, g=0, b=0)
+# )
+# colours <- list(
+#   "1"=list(r=255, g=0, b=0),
+#   "2"=list(r=0, g=0, b=255),
+#   "4"=list(r=0, g=0, b=0),
+#   "8"=list(r=255, g=255, b=255),
+#   "16"=list(r=0, g=255, b=0),
+#   "32"=list(r=255, g=0, b=255),
+#   "64"=list(r=255, g=0, b=255),
+#   "128"=list(r=0, g=255, b=255)
+# )
+colours <- list(
+  "1" = list(r = 98, g = 45, b = 209),
+  "2" = list(r = 255, g = 255, b = 0),
+  "4" = list(r = 113, g = 143, b = 200),
+  "8" = list(r = 255, g = 0, b = 0),
+  "16" = list(r = 0, g = 0, b = 255),
+  "32" = list(r = 255, g = 144, b = 0),
+  "64" = list(r = 144, g = 0, b = 144),
+  "128" = list(r = 144, g = 255, b = 0)
+)
+results <- rbind(NULL, c(Value = 0, Red = 0, Green = 0, Blue = 0))
+for (i in 1:255) {
+  if (!(i %in% names(colours))) {
+    cur <- i
+    count <- 0
+    r <- 0
+    g <- 0
+    b <- 0
+    print(i)
+    for (j in 1:(i - 1)) {
+      if (bitwAnd(cur, j) > 0) {
+        print(paste0("Matches ", j))
+        co <- colours[[as.character(j)]]
+        print(co)
+        r <- r + co$r
+        g <- g + co$g
+        b <- b + co$b
+        count <- count + 1
+        cur <- bitwAnd(cur, bitwNot(j))
+      }
+    }
+    r <- as.integer(r / count)
+    g <- as.integer(g / count)
+    b <- as.integer(b / count)
+    colours[[as.character(i)]] <- list(r = r, g = g, b = b)
+  }
+  v <- colours[[as.character(i)]]
+  results <- rbind(results, c(Value = i, Red = v$r, Green = v$g, Blue = v$b))
+}
+
+to_hex <- function(r, g, b) {
+  return(rgb(r, g, b, maxColorValue = 255))
+}
+hex <- list()
+for (i in 1:(length(colours))) {
+  v <- colours[[as.character(i)]]
+  hex[as.character(i)] <- to_hex(v$r, v$g, v$b)
+}
+# hex <- sapply(colours,
+#               function(v) {
+#                 print(v)
+#                 return(to_hex(v$r, v$g, v$b))
+#               })
+hex <- unlist(hex)
+hex <- c(to_hex(0, 0, 0), hex)
+names(hex) <- 0:255
+results <- as.data.table(results)
+write.table(results, file = "source.clr", col.names = FALSE, row.names = FALSE, sep = " ")
+
+imgs <- list.files(DIR, pattern = "^source.tif$", full.names = TRUE, recursive = TRUE)
+imgs <- imgs[grep("C2", imgs)]
+rasters <- lapply(imgs, function(img) {
+  return(trim(rast(img)))
+})
+e <- ext(rasters[[1]])
+for (r in rasters) {
+  e <- ext(extend(r, e))
+}
+pdf(FILE_OUT, width = 11, height = 8.5)
+par(mfrow = c(3, 4), mar = c(2, 2, 2, 2))
+for (i in 1:length(rasters)) {
+  img <- imgs[[i]]
+  title <- basename(dirname(img))
+  r <- rasters[[i]]
+  # set largest extent between all rasters
+  r <- extend(r, e)
+  # HACK: have to filter to values in the raster so it picks right colours
+  # col <- hex[1 + unique(r)]
+  # col <- hex[unique(r)]
+  # col <- hex[1:max(unique(r))]
+  # col <- hex[as.character(unique(r))]
+  col <- hex[as.integer(unlist(unique(r)))]
+  # plot(r, col=col, main=title)
+  # plot(ratify(r), col=hex, main=title)
+  # colortable(r) <- hex
+  # plot(r, main=title)
+  # print(levelplot(ratify(r), col.regions=col, att='ID'))
+  r <- as.factor(r)
+  # lvls <- levels(r)[[1]]
+  # lvls$ID <- unique(r)
+  # levels(r) <- lvls
+  # print(levelplot(r, col.regions = col, att = "ID", main = title))
+  plot(r, col=col, main=title)
+}
+dev.off()

tbd/CMakeLists.txt

@@ -11,6 +11,7 @@ message("CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}")
 
 # set(CMAKE_CXX_FLAGS "-Wall -Wextra -std=c++23")
 set(CMAKE_CXX_FLAGS "-Wall -Wextra")
+# set(CMAKE_CXX_FLAGS "-Wall -Wextra -Wfatal-errors")
 add_compile_options(
         "$<$<CONFIG:RELEASE>:-O3;>"
         "$<$<CONFIG:DEBUG>:-O0;-fno-omit-frame-pointer;-g;>"

tbd/mk_test.sh

@@ -0,0 +1,3 @@
+#!/bin/bash
+scripts/mk_clean.sh Release
+./tbd test ../data/test_output all

tbd/prep.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+set -x
+
+d="collected"
+g="11N_47610"
+r="m3_202308021910"
+
+scripts/mk_clean.sh
+
+rm -rf /appl/data/${d}
+
+dir_to="/appl/data/${d}/${g}/${r}"
+pushd /appl/data/sims/${r}/sims/${g}
+mkdir -p ${dir_to}
+cp * ${dir_to}/
+popd
+
+./get_grids.sh 00_pc_or_pdf_is_nonfuel
+./get_grids.sh 00_pc_or_pdf_is_d1
+./get_grids.sh 00_pc_or_pdf_is_d2
+./get_grids.sh normal
+
+pushd /appl/data/${d}
+FILE_7Z="${g}.7z"
+rm -rf ${FILE_7Z}
+7za a ${FILE_7Z} ${g}/*
+popd

tbd/scripts/mk.sh

@@ -0,0 +1,10 @@
+#!/bin/bash
+# didn't figure out how to do this with cmake yet but this works for now
+DIR_BUILD=/appl/tbd/build
+VARIANT="$*"
+if [ -z "${VARIANT}" ]; then
+  VARIANT=Release
+fi
+echo Set VARIANT=${VARIANT}
+/usr/bin/cmake --no-warn-unused-cli -DCMAKE_EXPORT_COMPILE_COMMANDS:BOOL=TRUE -DCMAKE_BUILD_TYPE:STRING=${VARIANT} -S/appl/tbd -B${DIR_BUILD} -G "Unix Makefiles" \
+  && /usr/bin/cmake --build ${DIR_BUILD} --config ${VARIANT} --target all -j 50 --

tbd/src/cpp/Cell.h

@@ -1,17 +1,6 @@
-// Copyright (c) 2020-2021, Queen's Printer for Ontario.
-//
-// This program is free software: you can redistribute it and/or modify
-// it under the terms of the GNU Affero General Public License as
-// published by the Free Software Foundation, either version 3 of the
-// License, or (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU Affero General Public License for more details.
-//
-// You should have received a copy of the GNU Affero General Public License
-// along with this program.  If not, see <https://www.gnu.org/licenses/>.
+/* Copyright (c) 2020,  Queen's Printer for Ontario */
+
+/* SPDX-License-Identifier: AGPL-3.0-or-later */
 
 #pragma once
 #include <limits>