refactor: addLakesInDeepDepressions

src/scripts/generation/grid/grid.ts

@@ -0,0 +1,72 @@
+import {defineMapSize} from "modules/coordinates";
+import {generateGrid} from "scripts/generation/graph";
+import {markupGridFeatures} from "scripts/generation/markup";
+import {rn} from "utils/numberUtils";
+import {byId} from "utils/shorthands";
+import {generatePrecipitation} from "./precipitation";
+import {calculateTemperatures} from "./temperature";
+
+const {HeightmapGenerator} = window;
+
+export async function createGrid(globalGrid: IGrid, precreatedGrid?: IGrid): Promise<IGrid> {
+  const shouldRegenerate = shouldRegenerateGridPoints(globalGrid);
+  const {spacing, cellsDesired, boundary, points, cellsX, cellsY, cells, vertices} = shouldRegenerate
+    ? (precreatedGrid && undressPrecreatedGrid(precreatedGrid)) || generateGrid()
+    : undressPrecreatedGrid(globalGrid);
+
+  const heights: Uint8Array = await HeightmapGenerator.generate({
+    vertices,
+    points,
+    cells,
+    cellsDesired,
+    spacing,
+    cellsX,
+    cellsY
+  });
+  if (!heights) throw new Error("Heightmap generation failed");
+
+  const {featureIds, distanceField, features} = markupGridFeatures(cells.c, cells.b, heights);
+
+  const touchesEdges = features.some(feature => feature && feature.land && feature.border);
+  defineMapSize(touchesEdges);
+
+  const temp = calculateTemperatures(heights, cellsX, points);
+  const prec = generatePrecipitation(heights, temp, cellsX, cellsY);
+
+  return {
+    cellsDesired,
+    cellsX,
+    cellsY,
+    spacing,
+    boundary,
+    points,
+    vertices,
+    cells: {
+      ...cells,
+      h: heights,
+      f: featureIds,
+      t: distanceField,
+      prec,
+      temp
+    },
+    features
+  };
+}
+
+function undressPrecreatedGrid(extendedGrid: IGrid) {
+  const {spacing, cellsDesired, boundary, points, cellsX, cellsY, cells, vertices} = extendedGrid;
+  const {i, b, c, v} = cells;
+  return {spacing, cellsDesired, boundary, points, cellsX, cellsY, cells: {i, b, c, v}, vertices};
+}
+
+// check if new grid graph should be generated or we can use the existing one
+export function shouldRegenerateGridPoints(grid: IGrid) {
+  const cellsDesired = Number(byId("pointsInput")?.dataset.cells);
+  if (cellsDesired !== grid.cellsDesired) return true;
+
+  const newSpacing = rn(Math.sqrt((graphWidth * graphHeight) / cellsDesired), 2);
+  const newCellsX = Math.floor((graphWidth + 0.5 * newSpacing - 1e-10) / newSpacing);
+  const newCellsY = Math.floor((graphHeight + 0.5 * newSpacing - 1e-10) / newSpacing);
+
+  return grid.spacing !== newSpacing || grid.cellsX !== newCellsX || grid.cellsY !== newCellsY;
+}

src/scripts/generation/grid/lakes.ts

@@ -0,0 +1,117 @@
+import {TIME} from "config/logging";
+import {getInputNumber, getInputValue} from "utils/nodeUtils";
+import {DISTANCE_FIELD, MAX_HEIGHT, MIN_LAND_HEIGHT} from "config/generation";
+import {drawPolygon} from "utils/debugUtils";
+
+const {LAND_COAST, WATER_COAST} = DISTANCE_FIELD;
+
+// near sea lakes usually get a lot of water inflow
+// most of them would brake threshold and flow out to sea (see Ancylus Lake)
+// connect these type of lakes to the main water body to improve the heightmap
+export function openNearSeaLakes(grid: IGraph & Partial<IGrid>) {
+  if (getInputValue("templateInput") === "Atoll") return; // no need for Atolls
+
+  const {cells, features} = grid;
+  if (!features?.find(f => f && f.type === "lake")) return; // no lakes
+
+  TIME && console.time("openNearSeaLakes");
+  const LIMIT = 22; // max height that can be breached by water
+
+  const isLake = (featureId: number) => featureId && (features[featureId] as IGridFeature).type === "lake";
+  const isOcean = (featureId: number) => featureId && (features[featureId] as IGridFeature).type === "ocean";
+
+  for (const cellId of cells.i) {
+    const featureId = cells.f[cellId];
+    if (!isLake(featureId)) continue; // not a lake cell
+
+    check_neighbours: for (const neibCellId of cells.c[cellId]) {
+      // water cannot brake the barrier
+      if (cells.t[neibCellId] !== WATER_COAST || cells.h[neibCellId] > LIMIT) continue;
+
+      for (const neibOfNeibCellId of cells.c[neibCellId]) {
+        const neibOfNeibFeatureId = cells.f[neibOfNeibCellId];
+        if (!isOcean(neibOfNeibFeatureId)) continue; // not an ocean
+        removeLake(neibCellId, featureId, neibOfNeibFeatureId);
+        break check_neighbours;
+      }
+    }
+  }
+
+  function removeLake(barrierCellId: number, lakeFeatureId: number, oceanFeatureId: number) {
+    cells.h[barrierCellId] = MIN_LAND_HEIGHT - 1;
+    cells.t[barrierCellId] = WATER_COAST;
+    cells.f[barrierCellId] = oceanFeatureId;
+
+    for (const neibCellId of cells.c[barrierCellId]) {
+      if (cells.h[neibCellId] >= MIN_LAND_HEIGHT) cells.t[neibCellId] = LAND_COAST;
+    }
+
+    if (features && lakeFeatureId) {
+      // mark former lake as ocean
+      (features[lakeFeatureId] as IGridFeature).type = "ocean";
+    }
+  }
+
+  TIME && console.timeEnd("openNearSeaLakes");
+}
+
+// some deeply depressed areas may not be resolved on river generation
+// this areas tend to collect precipitation, so we can add a lake there to help the resolver
+export function addLakesInDeepDepressions(
+  heights: Uint8Array,
+  neighbours: number[][],
+  cellVertices: number[][],
+  vertices: IGraphVertices,
+  indexes: UintArray
+) {
+  const ELEVATION_LIMIT = getInputNumber("lakeElevationLimitOutput");
+  if (ELEVATION_LIMIT === MAX_HEIGHT - MIN_LAND_HEIGHT) return heights; // any depression can be resolved
+
+  TIME && console.time("addLakesInDeepDepressions");
+
+  const landCells = indexes.filter(i => heights[i] >= MIN_LAND_HEIGHT);
+  landCells.sort((a, b) => heights[a] - heights[b]); // lower elevation first
+
+  const currentHeights = new Uint8Array(heights);
+  const checkedCells: Dict<true> = {[landCells[0]]: true};
+
+  for (const cellId of landCells) {
+    if (checkedCells[cellId]) continue;
+
+    const THESHOLD_HEIGHT = currentHeights[cellId] + ELEVATION_LIMIT;
+
+    let inDeepDepression = true;
+
+    const queue = [cellId];
+    const checkedPaths: Dict<true> = {[cellId]: true};
+
+    while (queue.length) {
+      const nextCellId = queue.pop()!;
+
+      if (currentHeights[nextCellId] < MIN_LAND_HEIGHT) {
+        inDeepDepression = false;
+        break;
+      }
+
+      for (const neibCellId of neighbours[nextCellId]) {
+        if (checkedPaths[neibCellId]) continue;
+
+        checkedPaths[neibCellId] = true;
+        checkedCells[neibCellId] = true;
+
+        if (currentHeights[neibCellId] < THESHOLD_HEIGHT) queue.push(neibCellId);
+      }
+    }
+
+    if (inDeepDepression) {
+      currentHeights[cellId] = MIN_LAND_HEIGHT - 1;
+      console.log(`ⓘ Added lake at deep depression. Cell: ${cellId}`);
+
+      const polygon = cellVertices[cellId].map(vertex => vertices.p[vertex]);
+      drawPolygon(polygon, {stroke: "red", strokeWidth: 1, fill: "none"});
+    }
+  }
+
+  TIME && console.timeEnd("addLakesInDeepDepressions");
+  return currentHeights;
+}

src/scripts/generation/grid/precipitation.ts

@@ -0,0 +1,170 @@
+// @ts-nocheck
+import * as d3 from "d3";
+
+import {TIME} from "config/logging";
+import {minmax} from "utils/numberUtils";
+import {rand} from "utils/probabilityUtils";
+import {getInputNumber, getInputValue} from "utils/nodeUtils";
+import {byId} from "utils/shorthands";
+
+// simplest precipitation model
+export function generatePrecipitation(heights: Uint8Array, temperatures: Int8Array, cellsX: number, cellsY: number) {
+  TIME && console.time("generatePrecipitation");
+  prec.selectAll("*").remove();
+
+  const precipitation = new Uint8Array(heights.length); // precipitation array
+
+  const cellsNumberModifier = (byId("pointsInput").dataset.cells / 10000) ** 0.25;
+  const precInputModifier = getInputNumber("precInput") / 100;
+  const modifier = cellsNumberModifier * precInputModifier;
+
+  const westerly = [];
+  const easterly = [];
+  let southerly = 0;
+  let northerly = 0;
+
+  // precipitation modifier per latitude band
+  // x4 = 0-5 latitude: wet through the year (rising zone)
+  // x2 = 5-20 latitude: wet summer (rising zone), dry winter (sinking zone)
+  // x1 = 20-30 latitude: dry all year (sinking zone)
+  // x2 = 30-50 latitude: wet winter (rising zone), dry summer (sinking zone)
+  // x3 = 50-60 latitude: wet all year (rising zone)
+  // x2 = 60-70 latitude: wet summer (rising zone), dry winter (sinking zone)
+  // x1 = 70-85 latitude: dry all year (sinking zone)
+  // x0.5 = 85-90 latitude: dry all year (sinking zone)
+  const latitudeModifier = [4, 2, 2, 2, 1, 1, 2, 2, 2, 2, 3, 3, 2, 2, 1, 1, 1, 0.5];
+  const MAX_PASSABLE_ELEVATION = 85;
+
+  // define wind directions based on cells latitude and prevailing winds there
+  d3.range(0, heights.length, cellsX).forEach(function (c, i) {
+    const lat = mapCoordinates.latN - (i / cellsY) * mapCoordinates.latT;
+    const latBand = ((Math.abs(lat) - 1) / 5) | 0;
+    const latMod = latitudeModifier[latBand];
+    const windTier = (Math.abs(lat - 89) / 30) | 0; // 30d tiers from 0 to 5 from N to S
+    const {isWest, isEast, isNorth, isSouth} = getWindDirections(windTier);
+
+    if (isWest) westerly.push([c, latMod, windTier]);
+    if (isEast) easterly.push([c + cellsX - 1, latMod, windTier]);
+    if (isNorth) northerly++;
+    if (isSouth) southerly++;
+  });
+
+  // distribute winds by direction
+  if (westerly.length) passWind(westerly, 120 * modifier, 1, cellsX);
+  if (easterly.length) passWind(easterly, 120 * modifier, -1, cellsX);
+
+  const vertT = southerly + northerly;
+  if (northerly) {
+    const bandN = ((Math.abs(mapCoordinates.latN) - 1) / 5) | 0;
+    const latModN = mapCoordinates.latT > 60 ? d3.mean(latitudeModifier) : latitudeModifier[bandN];
+    const maxPrecN = (northerly / vertT) * 60 * modifier * latModN;
+    passWind(d3.range(0, cellsX, 1), maxPrecN, cellsX, cellsY);
+  }
+
+  if (southerly) {
+    const bandS = ((Math.abs(mapCoordinates.latS) - 1) / 5) | 0;
+    const latModS = mapCoordinates.latT > 60 ? d3.mean(latitudeModifier) : latitudeModifier[bandS];
+    const maxPrecS = (southerly / vertT) * 60 * modifier * latModS;
+    passWind(d3.range(heights.length - cellsX, heights.length, 1), maxPrecS, -cellsX, cellsY);
+  }
+
+  function getWindDirections(tier) {
+    const angle = options.winds[tier];
+
+    const isWest = angle > 40 && angle < 140;
+    const isEast = angle > 220 && angle < 320;
+    const isNorth = angle > 100 && angle < 260;
+    const isSouth = angle > 280 || angle < 80;
+
+    return {isWest, isEast, isNorth, isSouth};
+  }
+
+  function passWind(source, maxPrec, next, steps) {
+    const maxPrecInit = maxPrec;
+
+    for (let first of source) {
+      if (first[0]) {
+        maxPrec = Math.min(maxPrecInit * first[1], 255);
+        first = first[0];
+      }
+
+      let humidity = maxPrec - heights[first]; // initial water amount
+      if (humidity <= 0) continue; // if first cell in row is too elevated consider wind dry
+
+      for (let s = 0, current = first; s < steps; s++, current += next) {
+        if (temperatures[current] < -5) continue; // no flux in permafrost
+
+        if (heights[current] < 20) {
+          // water cell
+          if (heights[current + next] >= 20) {
+            precipitation[current + next] += Math.max(humidity / rand(10, 20), 1); // coastal precipitation
+          } else {
+            humidity = Math.min(humidity + 5 * modifier, maxPrec); // wind gets more humidity passing water cell
+            precipitation[current] += 5 * modifier; // water cells precipitation (need to correctly pour water through lakes)
+          }
+          continue;
+        }
+
+        // land cell
+        const isPassable = heights[current + next] <= MAX_PASSABLE_ELEVATION;
+        const cellPrec = isPassable ? getPrecipitation(humidity, current, next) : humidity;
+        precipitation[current] += cellPrec;
+        const evaporation = cellPrec > 1.5 ? 1 : 0; // some humidity evaporates back to the atmosphere
+        humidity = isPassable ? minmax(humidity - cellPrec + evaporation, 0, maxPrec) : 0;
+      }
+    }
+  }
+
+  function getPrecipitation(humidity, i, n) {
+    const normalLoss = Math.max(humidity / (10 * modifier), 1); // precipitation in normal conditions
+    const diff = Math.max(heights[i + n] - heights[i], 0); // difference in height
+    const mod = (heights[i + n] / 70) ** 2; // 50 stands for hills, 70 for mountains
+    return minmax(normalLoss + diff * mod, 1, humidity);
+  }
+
+  TIME && console.timeEnd("generatePrecipitation");
+  return precipitation;
+}
+
+// TODO: move to renderer
+function drawWindDirection() {
+  const wind = prec.append("g").attr("id", "wind");
+
+  d3.range(0, 6).forEach(function (t) {
+    if (westerly.length > 1) {
+      const west = westerly.filter(w => w[2] === t);
+      if (west && west.length > 3) {
+        const from = west[0][0];
+        const to = west[west.length - 1][0];
+        const y = (grid.points[from][1] + grid.points[to][1]) / 2;
+        wind.append("text").attr("x", 20).attr("y", y).text("\u21C9");
+      }
+    }
+    if (easterly.length > 1) {
+      const east = easterly.filter(w => w[2] === t);
+      if (east && east.length > 3) {
+        const from = east[0][0];
+        const to = east[east.length - 1][0];
+        const y = (grid.points[from][1] + grid.points[to][1]) / 2;
+        wind
+          .append("text")
+          .attr("x", graphWidth - 52)
+          .attr("y", y)
+          .text("\u21C7");
+      }
+    }
+  });
+
+  if (northerly)
+    wind
+      .append("text")
+      .attr("x", graphWidth / 2)
+      .attr("y", 42)
+      .text("\u21CA");
+  if (southerly)
+    wind
+      .append("text")
+      .attr("x", graphWidth / 2)
+      .attr("y", graphHeight - 20)
+      .text("\u21C8");
+}

src/scripts/generation/grid/temperature.ts

@@ -0,0 +1,44 @@
+import * as d3 from "d3";
+
+import {TIME} from "config/logging";
+import {minmax} from "utils/numberUtils";
+import {getInputNumber} from "utils/nodeUtils";
+import {MIN_LAND_HEIGHT} from "config/generation";
+
+const interpolate = d3.easePolyInOut.exponent(0.5); // interpolation function
+
+export function calculateTemperatures(heights: Uint8Array, cellsX: number, points: TPoints) {
+  TIME && console.time("calculateTemperatures");
+
+  const temperatures = new Int8Array(heights.length); // temperature array
+
+  // temperature decreases by 6.5 Celsius per kilometer
+  const heightExponent = getInputNumber("heightExponentInput");
+  function decreaseTempFromElevation(height: number) {
+    if (height < MIN_LAND_HEIGHT) return 0;
+
+    const realHeight = Math.pow(height - 18, heightExponent);
+    return (realHeight / 1000) * 6.5;
+  }
+
+  const tEq = getInputNumber("temperatureEquatorInput");
+  const tPole = getInputNumber("temperaturePoleInput");
+  const tDelta = tEq - tPole;
+
+  const {latN, latT} = window.mapCoordinates;
+
+  d3.range(0, heights.length, cellsX).forEach(rowStart => {
+    const y = points[rowStart][1];
+    const lat = Math.abs(latN - (y / graphHeight) * latT); // [0; 90]
+
+    const initTemp = tEq - interpolate(lat / 90) * tDelta;
+    for (let i = rowStart; i < rowStart + cellsX; i++) {
+      const elevationDecrease = decreaseTempFromElevation(heights[i]);
+      temperatures[i] = minmax(initTemp - elevationDecrease, -128, 127);
+    }
+  });
+
+  TIME && console.timeEnd("calculateTemperatures");
+
+  return temperatures;
+}