earth weather

src/mapgen/earth/png_holder.h

@@ -0,0 +1,181 @@
+#pragma once
+
+#include <png.h>
+#include <cstdio>
+#include <string>
+#include <vector>
+#include <stdexcept>
+#include <optional>
+#include <cstdint>
+//#include <iostream>
+#include "SColor.h"
+
+//struct RGBA {uint8_t r, g, b, a;};
+using RGBA = video::SColor;
+
+class PngImage {
+public:
+    // open and read PNG, convert to 8-bit RGBA
+    explicit PngImage(const std::string& filename) {
+        fp_ = std::fopen(filename.c_str(), "rb");
+        if (!fp_) throw std::runtime_error("Failed to open file: " + filename);
+
+        png_ptr_ = png_create_read_struct(PNG_LIBPNG_VER_STRING, nullptr, nullptr, nullptr);
+        if (!png_ptr_) {
+            std::fclose(fp_);
+            fp_ = nullptr;
+            throw std::runtime_error("png_create_read_struct failed");
+        }
+
+        info_ptr_ = png_create_info_struct(png_ptr_);
+        if (!info_ptr_) {
+            png_destroy_read_struct(&png_ptr_, nullptr, nullptr);
+            std::fclose(fp_);
+            fp_ = nullptr;
+            throw std::runtime_error("png_create_info_struct failed");
+        }
+
+        if (setjmp(png_jmpbuf(png_ptr_))) {
+            // libpng error
+            cleanup();
+            throw std::runtime_error("libpng error during init/read");
+        }
+
+        png_init_io(png_ptr_, fp_);
+        png_read_info(png_ptr_, info_ptr_);
+
+        png_uint_32 w, h;
+        int bit_depth, color_type, interlace_type, compression_type, filter_method;
+        png_get_IHDR(png_ptr_, info_ptr_, &w, &h, &bit_depth, &color_type,
+                     &interlace_type, &compression_type, &filter_method);
+
+        width_ = static_cast<int>(w);
+        height_ = static_cast<int>(h);
+
+        // Transformations to 8-bit RGBA:
+        // Expand palette, grayscale < 8, and tRNS to RGB(A)
+        png_set_expand(png_ptr_);
+
+        // Convert 16-bit to 8-bit
+        if (bit_depth == 16) png_set_strip_16(png_ptr_);
+
+        // Convert gray/gray+alpha to RGB(A)
+        if (color_type == PNG_COLOR_TYPE_GRAY || color_type == PNG_COLOR_TYPE_GRAY_ALPHA)
+            png_set_gray_to_rgb(png_ptr_);
+
+        // Ensure we have alpha: if image has no alpha, add an opaque alpha channel
+        if (!(color_type & PNG_COLOR_MASK_ALPHA)) {
+            // Put alpha after RGB (i.e., RGBA)
+            png_set_filler(png_ptr_, 0xFF, PNG_FILLER_AFTER);
+        }
+
+        // Handle interlace
+        png_set_interlace_handling(png_ptr_);
+
+        // Update info after transformations
+        png_read_update_info(png_ptr_, info_ptr_);
+
+        // Get bytes per row (after transforms)
+        png_uint_32 rowbytes = png_get_rowbytes(png_ptr_, info_ptr_);
+        if (rowbytes == 0) {
+            cleanup();
+            throw std::runtime_error("Invalid rowbytes");
+        }
+
+        rowbytes_ = static_cast<size_t>(rowbytes);
+        // allocate one contiguous buffer for all image data
+        image_data_.resize(rowbytes_ * static_cast<size_t>(height_));
+        row_pointers_.resize(height_);
+
+        for (int i = 0; i < height_; ++i) {
+            row_pointers_[i] = image_data_.data() + (size_t)i * rowbytes_;
+        }
+
+        // Read the image into our buffer
+        png_read_image(png_ptr_, row_pointers_.data());
+
+        // After read, we expect 4 bytes per pixel (RGBA)
+        // But verify: rowbytes must be width * 4 (or at least enough)
+        if (rowbytes_ < static_cast<size_t>(width_) * 4) {
+            cleanup();
+            throw std::runtime_error("Unexpected row size after transformation");
+        }
+    }
+
+    // Disable copy
+    PngImage(const PngImage&) = delete;
+    PngImage& operator=(const PngImage&) = delete;
+
+    // Enable move
+    PngImage(PngImage&& other) noexcept {
+        steal_from(std::move(other));
+    }
+    PngImage& operator=(PngImage&& other) noexcept {
+        if (this != &other) {
+            cleanup();
+            steal_from(std::move(other));
+        }
+        return *this;
+    }
+
+    ~PngImage() {
+        cleanup();
+    }
+
+    int width() const noexcept { return width_; }
+    int height() const noexcept { return height_; }
+
+    // Return RGBA for given coords, or std::nullopt if out of bounds.
+    std::optional<RGBA> get_pixel(int x, int y) const noexcept {
+        if (x < 0 || y < 0 || x >= width_ || y >= height_) return std::nullopt;
+        const uint8_t* row = row_pointers_[y];
+        const uint8_t* px = row + static_cast<size_t>(x) * 4; // 4 bytes per pixel (R,G,B,A)
+        RGBA c{ px[3], px[0], px[1], px[2] };
+        return c;
+    }
+
+private:
+    png_structp png_ptr_ = nullptr;
+    png_infop info_ptr_ = nullptr;
+    FILE* fp_ = nullptr;
+
+    int width_ = 0;
+    int height_ = 0;
+    size_t rowbytes_ = 0;
+
+    std::vector<uint8_t> image_data_;
+    std::vector<png_bytep> row_pointers_;
+
+    void cleanup() noexcept {
+        if (png_ptr_ || info_ptr_) {
+            // png_destroy_read_struct frees info_ptr_ as well
+            png_destroy_read_struct(&png_ptr_, &info_ptr_, nullptr);
+            png_ptr_ = nullptr;
+            info_ptr_ = nullptr;
+        }
+        if (fp_) {
+            std::fclose(fp_);
+            fp_ = nullptr;
+        }
+        image_data_.clear();
+        row_pointers_.clear();
+    }
+
+    void steal_from(PngImage&& other) noexcept {
+        png_ptr_ = other.png_ptr_;
+        info_ptr_ = other.info_ptr_;
+        fp_ = other.fp_;
+        width_ = other.width_;
+        height_ = other.height_;
+        rowbytes_ = other.rowbytes_;
+        image_data_ = std::move(other.image_data_);
+        row_pointers_ = std::move(other.row_pointers_);
+
+        other.png_ptr_ = nullptr;
+        other.info_ptr_ = nullptr;
+        other.fp_ = nullptr;
+        other.width_ = other.height_ = 0;
+        other.rowbytes_ = 0;
+    }
+};
+

src/mapgen/earth/rgb_temp.h

@@ -0,0 +1,116 @@
+// write function on c++ : params color in rgb encoding temperature from blue to red and returns temperature in celsius
+
+/*
+I'm considering how to solve for inverse mapping by finding the parameter \( t \) that minimizes error. A brute-force search over the range [0, 1] with small steps or Newton's method could work, but implementing a numerical search sounds more robust. I can sample 1,000 points and refine via binary search for better matching. I’ll provide both a straightforward linear method and a more intricate jet inverse approach using C++17 with relevant examples. I’ll format everything in markdown!**Coding RGB to Temperature Conversion**
+
+I’m writing code that includes the `cmath` and `algorithm` libraries. First, I’ll define a function, `rgbToTemperatureLinear`, that converts RGB values to a temperature range based on given minimum and maximum temperatures. I’ll also implement `rgbToTemperatureJet`, which uses formulas to calculate values for red, green, and blue based on a temperature parameter \( t \). Next, I’ll execute a coarse search with 2,000 samples to find the best \( t \), and refine that with smaller steps to minimize the squared error.**Refining the Temperature Mapping**
+
+Next, I plan to refine the search around the best \( t \) value found earlier by using either the golden-section method or a simpler approach with a loop from 0 to 5. I’ll search within the range defined by \( t_{\text{best}} \) plus and minus a small delta, using a finer step of \( K = 200 \) to compute any necessary updates. As I decrease the delta for convergence, I'll ensure the final result is clear: returning the calculated temperature. 
+
+I'll also include usage examples and note important considerations like potential ambiguities in mapping and the presence of noise.</div>Below are two C++ functions you can use to convert an RGB color (assumed to encode a temperature on a blue→red scale) back to a temperature in Celsius.
+
+- The first function (rgbToCelsiusLinear) assumes a simple linear gradient from pure blue (0,0,255) at the cold end to pure red (255,0,0) at the hot end. It is cheap and analytic.
+- The second function (rgbToCelsiusJet) handles a common “rainbow/jet-like” thermal map (blue→cyan→green→yellow→red) by numerically inverting the forward jet colormap. This is a more realistic thermal palette but requires sampling/search so it’s a bit heavier.
+
+Both functions take r,g,b in 0..255 and minTemp/maxTemp in Celsius (defaults given) and return a double temperature in Celsius.
+
+Code (C++17):
+*/
+#pragma once
+#include <cmath>
+#include <algorithm>
+#include <array>
+#include <limits>
+
+// clamp helper
+static inline double clamp01(double x) {
+    return std::max(0.0, std::min(1.0, x));
+}
+
+// 1) Simple linear blue->red gradient inversion
+// Maps: blue (0,0,255) -> t=0  ; red (255,0,0) -> t=1
+// r,g,b are 0..255, returns temperature in Celsius between minTemp and maxTemp.
+double rgbToCelsiusLinear(int r, int g, int b, double minTemp = -30.0, double maxTemp = 50.0) {
+    double rn = clamp01(r / 255.0);
+    double bn = clamp01(b / 255.0);
+
+    // projection onto the R<-B gradient: t = (r_norm - b_norm + 1) / 2
+    double t = (rn - bn + 1.0) / 2.0;
+    t = clamp01(t);
+
+    return minTemp + t * (maxTemp - minTemp);
+}
+
+// 2) Invert a typical "jet" / rainbow colormap numerically.
+// Forward jet mapping (t in [0,1]) used here:
+// r = clamp(1.5 - |4t - 3|)
+// g = clamp(1.5 - |4t - 2|)
+// b = clamp(1.5 - |4t - 1|)
+static inline std::array<double,3> jetColor(double t) {
+    t = clamp01(t);
+    double r = clamp01(1.5 - std::fabs(4.0 * t - 3.0));
+    double g = clamp01(1.5 - std::fabs(4.0 * t - 2.0));
+    double b = clamp01(1.5 - std::fabs(4.0 * t - 1.0));
+    return {r, g, b};
+}
+
+// invert jet colormap by coarse-to-fine search
+double rgbToCelsiusJet(int r, int g, int b, double minTemp = -30.0, double maxTemp = 50.0) {
+    const double rn = clamp01(r / 255.0);
+    const double gn = clamp01(g / 255.0);
+    const double bn = clamp01(b / 255.0);
+
+    // coarse sampling to get initial guess
+    int coarseN = 800; // coarse samples (adjust for speed/precision)
+    double bestT = 0.0;
+    double bestErr = std::numeric_limits<double>::infinity();
+
+    for (int i = 0; i < coarseN; ++i) {
+        double t = double(i) / double(coarseN - 1);
+        auto c = jetColor(t);
+        double err = (c[0] - rn)*(c[0] - rn) + (c[1] - gn)*(c[1] - gn) + (c[2] - bn)*(c[2] - bn);
+        if (err < bestErr) {
+            bestErr = err;
+            bestT = t;
+        }
+    }
+
+    // refine around bestT with a few finer searches
+    double delta = 1.0 / double(coarseN - 1);
+    for (int pass = 0; pass < 4; ++pass) {
+        double start = clamp01(bestT - delta);
+        double end   = clamp01(bestT + delta);
+        int fineN = 200;
+        double localBestT = bestT;
+        double localBestErr = bestErr;
+        for (int i = 0; i < fineN; ++i) {
+            double t = start + (end - start) * double(i) / double(fineN - 1);
+            auto c = jetColor(t);
+            double err = (c[0] - rn)*(c[0] - rn) + (c[1] - gn)*(c[1] - gn) + (c[2] - bn)*(c[2] - bn);
+            if (err < localBestErr) {
+                localBestErr = err;
+                localBestT = t;
+            }
+        }
+        bestT = localBestT;
+        bestErr = localBestErr;
+        delta *= 0.25; // shrink search radius
+    }
+
+    double temp = minTemp + bestT * (maxTemp - minTemp);
+    return temp;
+}
+
+/*
+Usage examples:
+- If your color map is exactly linear blue->red use rgbToCelsiusLinear.
+- If your color map is a rainbow/thermal (blue→cyan→green→yellow→red) use rgbToCelsiusJet.
+
+Notes and caveats:
+- The mapping from color to temperature is not unique without knowing the exact colormap and the numeric temperature range. You must supply the min and max temperatures that correspond to the cold and hot ends of the colormap.
+- Real images may contain noise, gamma correction or encoding differences; you may need to pre-correct for gamma or color-space conversions if colors appear non-linear.
+- If your colormap is different, provide its forward mapping and invert it similarly (analytic or numeric search).
+- If you expect many conversions, precompute a lookup table (LUT) from colors to temperatures for speed.
+
+If you tell me which exact colormap and numeric temperature range you use (and whether input is linear RGB or sRGB), I can tailor the function exactly to that mapping.
+*/

src/mapgen/mapgen_earth.cpp

@@ -32,6 +32,9 @@ along with Freeminer.  If not, see <http://www.gnu.org/licenses/>.
 #include "constants.h"
 #include "debug/dump.h"
 #include "emerge.h"
+#include "mapgen/earth/png_holder.h"
+#include "mapgen/earth/rgb_temp.h"
+#include "server.h"
 #include "filesys.h"
 #include "irr_v2d.h"
 #include "irr_v3d.h"
@@ -206,6 +209,14 @@ MapgenEarth::MapgenEarth(MapgenEarthParams *params_, EmergeParams *emerge) :
 	hgt_reader.debug = 0;
 */
 
+	{
+		const auto heat_img = m_emerge->env->m_server->getGameSpec()->path + "/" +
+							  "mods/earth/data/earth_heat.png";
+		if (std::filesystem::exists(heat_img)) {
+			heat_image = std::make_unique<PngImage>(heat_img);
+		}
+	}
+
 	if (!maps_holder) {
 		maps_holder = std::make_unique<maps_holder_t>();
 	}
@@ -507,16 +518,42 @@ void MapgenEarth::generateBuildings()
 #endif
 }
 
+/*
+Heat data:
+https://search.earthdata.nasa.gov/search/granules?p=C1276812859-GES_DISC&pg[0][v]=f&pg[0][gsk]=-start_date&g=G3503129918-GES_DISC&ff=Map%20Imagery&tl=1038183927.842!5!!&fsm0=Clouds&fst0=Atmosphere&long=2.3696682464454994&zoom=2.8910241901494977
+https://gibs-a.earthdata.nasa.gov/wmts/epsg4326/best/wmts.cgi?TIME=2025-03-01&layer=MERRA2_2m_Air_Temperature_Monthly&style=default&tilematrixset=2km&Service=WMTS&Request=GetTile&Version=1.0.0&Format=image%2Fpng&TileMatrix=0&TileCol=0&TileRow=0
+
+Snow data:
+https://cmr.earthdata.nasa.gov/search/concepts/C3050353608-NSIDC_CPRD.html
+https://gibs-a.earthdata.nasa.gov/wmts/epsg4326/best/wmts.cgi?TIME=2025-03-01&layer=MODIS_Terra_L3_Snow_Cover_Monthly_Average_Pct&style=default&tilematrixset=2km&Service=WMTS&Request=GetTile&Version=1.0.0&Format=image%2Fpng&TileMatrix=0&TileCol=0&TileRow=0
+*/
+
 weather::heat_t MapgenEarth::calcBlockHeat(const v3pos_t &p, uint64_t seed,
 		float timeofday, float totaltime, bool use_weather)
 {
 #if USE_OSMIUM
 	const auto ll = pos_to_ll(p);
-	const auto tile = osmium::geom::Tile(3, osmium::Location(ll.lon, ll.lat));
+	if (heat_image) {
+		const auto x = heat_image->width() * ((ll.lon + 180) / 360);
+		const auto y =  (std::min<int>( heat_image->height(), heat_image->width()/1.6 )) * ((ll.lat + 90) / 180);
+		const auto pixel = heat_image->get_pixel(x, y);
+		auto heat = rgbToCelsiusJet(pixel->getRed(), pixel->getGreen(), pixel->getBlue());
+		heat += m_emerge->biomemgr->weather_heat_daily *
+				(sin(cycle_shift(timeofday, -0.25) * M_PI) - 0.5); //-64..0..34
+		heat += p.Y / m_emerge->biomemgr->weather_heat_height;
+		if (m_emerge->biomemgr->weather_hot_core &&
+				p.Y < -(WEATHER_LIMIT - m_emerge->biomemgr->weather_hot_core))
+			heat += 6000 *
+					(1.0 - ((float)(p.Y - -WEATHER_LIMIT) /
+								   m_emerge->biomemgr
+										   ->weather_hot_core)); //hot core, later via realms
+		return heat;
+	}
 #endif
 	return m_emerge->biomemgr->calcBlockHeat(p, seed, timeofday, totaltime, use_weather);
 }
 
+// TODO: use cloud data
 weather::humidity_t MapgenEarth::calcBlockHumidity(const v3pos_t &p, uint64_t seed,
 		float timeofday, float totaltime, bool use_weather)
 {

src/mapgen/mapgen_earth.h

@@ -34,6 +34,8 @@ along with Freeminer.  If not, see <http://www.gnu.org/licenses/>.
 #include "threading/concurrent_vector.h"
 #include "util/lrucache.hpp"
 
+class PngImage;
+
 typedef core::vector2d<double> v2d;
 
 //using ll_t = float;
@@ -117,6 +119,7 @@ public:
 			float totaltime, bool use_weather) override;
 	weather::humidity_t calcBlockHumidity(const v3pos_t &p, uint64_t seed,
 			float timeofday, float totaltime, bool use_weather) override;
+	std::unique_ptr<PngImage> heat_image;
 
 	struct Stat
 	{

src/mapgen/mg_biome.cpp

@@ -307,36 +307,45 @@ Biome *BiomeGenOriginal::calcBiomeFromNoise(float heat, float humidity, v3s16 po
 }
 
 // Freeminer Weather
-weather::heat_t BiomeManager::calcBlockHeat(const v3pos_t &p, uint64_t seed, float timeofday, float totaltime, bool use_weather) {
+weather::heat_t BiomeManager::calcBlockHeat(const v3pos_t &p, uint64_t seed,
+		float timeofday, float totaltime, bool use_weather)
+{
 	//variant 1: full random
 	//f32 heat = NoisePerlin3D(np_heat, p.X, env->getGameTime()/100, p.Z, seed);
 
 	//variant 2: season change based on default heat map
-	auto heat = NoiseFractal2D(&(mapgen_params->bparams->np_heat), p.X, p.Z, seed); // -30..20..70
+	auto heat = NoiseFractal2D(
+			&(mapgen_params->bparams->np_heat), p.X, p.Z, seed); // -30..20..70
 	// auto heat =calcHeatAtPoint(p);
 
 	if (use_weather) {
 		f32 seasonv = totaltime;
 		seasonv /= 86400 * year_days; // season change speed
-		seasonv += (f32)p.X / weather_heat_width; // you can walk to area with other season
+		seasonv +=
+				(f32)p.X / weather_heat_width; // you can walk to area with other season
 		seasonv = sin(seasonv * M_PI);
 		//heat += (weather_heat_season * (heat < offset ? 2 : 0.5)) * seasonv; // -60..0..30
-		heat += (weather_heat_season) * seasonv; // -60..0..30
+		heat += (weather_heat_season)*seasonv; // -60..0..30
 
 		// daily change, hotter at sun +4, colder at night -4
-		heat += weather_heat_daily * (sin(cycle_shift(timeofday, -0.25) * M_PI) - 0.5); //-64..0..34
+		heat += weather_heat_daily *
+				(sin(cycle_shift(timeofday, -0.25) * M_PI) - 0.5); //-64..0..34
 	}
 	heat += p.Y / weather_heat_height; // upper=colder, lower=hotter, 3c per 1000
 
-	if (weather_hot_core && p.Y < -(WEATHER_LIMIT-weather_hot_core))
-		heat += 6000 * (1.0-((float)(p.Y - -WEATHER_LIMIT)/weather_hot_core)); //hot core, later via realms
+	if (weather_hot_core && p.Y < -(WEATHER_LIMIT - weather_hot_core))
+		heat += 6000 * (1.0 - ((float)(p.Y - -WEATHER_LIMIT) /
+									  weather_hot_core)); //hot core, later via realms
 
 	return heat;
 }
 
-weather::humidity_t BiomeManager::calcBlockHumidity(const v3pos_t &p, uint64_t seed, float timeofday, float totaltime, bool use_weather) {
+weather::humidity_t BiomeManager::calcBlockHumidity(const v3pos_t &p, uint64_t seed,
+		float timeofday, float totaltime, bool use_weather)
+{
 
-	auto humidity = NoiseFractal2D(&(mapgen_params->bparams->np_humidity), p.X, p.Z, seed);
+	auto humidity =
+			NoiseFractal2D(&(mapgen_params->bparams->np_humidity), p.X, p.Z, seed);
 	// auto humidity = calcHumidityAtPoint(p);
 
 	if (use_weather) {
@@ -344,7 +353,8 @@ weather::humidity_t BiomeManager::calcBlockHumidity(const v3pos_t &p, uint64_t s
 		seasonv /= 86400 * weather_humidity_days; // bad weather change speed (2 days)
 		seasonv += (f32)p.Z / weather_humidity_width;
 		humidity += weather_humidity_season * sin(seasonv * M_PI);
-		humidity += weather_humidity_daily * (sin(cycle_shift(timeofday, -0.1) * M_PI) - 0.5);
+		humidity +=
+				weather_humidity_daily * (sin(cycle_shift(timeofday, -0.1) * M_PI) - 0.5);
 	}
 	humidity += p.Y / weather_humidity_height; // upper=dry, lower=wet, 3c per 1000