Initial v2 draft of litho1.0 implementation.

MFraters · MFraters · commit 87c889b0c4c2 · 2025-08-19T17:38:16.000+02:00
diff --git a/include/world_builder/types/one_of.h b/include/world_builder/types/one_of.h
@@ -43,6 +43,13 @@ namespace WorldBuilder
         OneOf(const Interface &type_1,
               const Interface &type_2);
 
+        /**
+         * Constructor for the declaration
+         */
+        OneOf(const Interface &type_1,
+              const Interface &type_2,
+              const Interface &type_3);
+
         /**
          * Constructor for cloning an array.
          */
diff --git a/include/world_builder/world.h b/include/world_builder/world.h
@@ -257,6 +257,11 @@ namespace WorldBuilder
        */
       std::vector<Point<2> > cross_section;
 
+      /**
+       * The path to the Litho1.0 data directory
+       */
+      std::string litho_1_0_path;
+
       /**
        * Todo
        */
diff --git a/source/world_builder/features/continental_plate_models/composition/uniform.cc b/source/world_builder/features/continental_plate_models/composition/uniform.cc
@@ -62,10 +62,14 @@ namespace WorldBuilder
                             "Uniform compositional model. Sets constant compositional field.");
 
           // Declare entries of this plugin
-          prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))),
+          prm.declare_entry("min depth", Types::OneOf(Types::Double(0),
+                                                      Types::Array(Types::ValueAtPoints(0., 2.)),
+                                                      Types::String("")),
                             "The depth in meters from which the composition of this feature is present.");
 
-          prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits<double>::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits<double>::max(), 2.))),
+          prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits<double>::max()),
+                                                      Types::Array(Types::ValueAtPoints(std::numeric_limits<double>::max(), 2.)),
+                                                      Types::String("")),
                             "The depth in meters to which the composition of this feature is present.");
 
           prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0),
diff --git a/source/world_builder/features/continental_plate_models/temperature/linear.cc b/source/world_builder/features/continental_plate_models/temperature/linear.cc
@@ -64,10 +64,12 @@ namespace WorldBuilder
                             "Linear temperature model. Can be set to use an adiabatic temperature at the boundaries.");
 
           // Declare entries of this plugin
-          prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))),
+          prm.declare_entry("min depth", Types::OneOf(Types::Double(0),
+                                                      Types::Array(Types::ValueAtPoints(0., 2.)),
+                                                      Types::String("")),
                             "The depth in meters from which the composition of this feature is present.");
 
-          prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits<double>::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits<double>::max(), 2.))),
+          prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits<double>::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits<double>::max(), 2.)),Types::String("")),
                             "The depth in meters to which the composition of this feature is present.");
 
           prm.declare_entry("top temperature", Types::Double(293.15),
diff --git a/source/world_builder/parameters.cc b/source/world_builder/parameters.cc
@@ -18,6 +18,7 @@
 */
 
 
+#include "world_builder/assert.h"
 #include "world_builder/features/continental_plate_models/composition/interface.h"
 #include "world_builder/features/continental_plate_models/velocity/interface.h"
 #include "world_builder/features/continental_plate_models/grains/interface.h"
@@ -517,6 +518,7 @@ namespace WorldBuilder
     // 2. One double provided: use the default value everywhere. Return first with one value and second with size 0.
     // 3. One value in a double array and no points provided: use that value everywhere. Return first with one value and second with size 0.
     // 4. Other: fill the vectors with the default value and addition points and then add new point.
+    // 5. A string provided: for example litho1 sediment 3, and a path provided?
     //    If a value without points is encountered, the additional points are used.
     std::pair<std::vector<double>,std::vector<double>> result;
 
@@ -525,9 +527,11 @@ namespace WorldBuilder
     // start with adding the additional points with the default value
     // to do this we need the default value
     double default_value = 0;
-    bool is_array = true;
+    enum InputTypes {DOUBLE_TYPE,ARRAY_TYPE,STRING_TYPE,MAX_INPUT_TYPES};
+    InputTypes input_type = InputTypes::MAX_INPUT_TYPES;
     if (Pointer((strict_base + "/" + name).c_str()).Get(parameters) != nullptr && Pointer((strict_base + "/" + name).c_str()).Get(parameters)->IsArray())
       {
+        input_type = InputTypes::ARRAY_TYPE;
         const std::string value_def_path = get_full_json_schema_path() + "/" + name + "/oneOf/1/items/items/anyOf/0/default value";
         Value *value_def = Pointer(value_def_path.c_str()).Get(declarations);
         WBAssertThrow(value_def != nullptr,
@@ -538,9 +542,15 @@ namespace WorldBuilder
         // So no try/catch needed.
         default_value = value_def->GetDouble();
       }
+    else if (Pointer((strict_base + "/" + name).c_str()).Get(parameters) != nullptr && Pointer((strict_base + "/" + name).c_str()).Get(parameters)->IsString())
+      {
+        input_type = InputTypes::STRING_TYPE;
+        // I don't know if a default value makes sense here
+
+      }
     else
       {
-        is_array = false;
+        input_type = InputTypes::DOUBLE_TYPE;
         Value *value_def = Pointer((get_full_json_schema_path() + "/" + name + "/oneOf/0/default value").c_str()).Get(declarations);
         WBAssertThrow(value_def != nullptr,
                       "internal error: could not retrieve the default value at: "
@@ -556,8 +566,8 @@ namespace WorldBuilder
     // check if there is a user defined value
     if (Pointer((strict_base + "/" + name).c_str()).Get(parameters) != nullptr)
       {
-        // there is a user defined value, so either case 2, 3 or 4.
-        if (is_array)
+        // there is a user defined value, so either case 2, 3 4, or 5.
+        if (input_type == InputTypes::ARRAY_TYPE)
           {
             Value *array = Pointer((strict_base  + "/" + name).c_str()).Get(parameters);
 
@@ -690,6 +700,147 @@ namespace WorldBuilder
                   }
               }
           }
+        else if (input_type == InputTypes::STRING_TYPE)
+          {
+            // case 5: a string. Get the string and decode it
+            std::string value = "";
+
+            try
+              {
+                value = Pointer((strict_base + "/" + name).c_str()).Get(parameters)->GetString();
+              }
+            catch (...)
+              {
+                WBAssertThrow(false, "Could not convert values of " << strict_base << "/" << name << " into a String. "
+                              << "The provided value was \"" <<  Pointer((strict_base + "/" + name).c_str()).Get(parameters)->GetString() << "\".");
+              }
+            if (value.substr(0,8) == "Litho1.0")
+              {
+                enum LayerType
+                {
+                  // ASTHENOSPHERE,
+                  LITHOSPHERE,
+                  CRUST3,
+                  CRUST2,
+                  CRUST1,
+                  SEDIMENT3,
+                  SEDIMENT2,
+                  SEDIMENT1,
+                  ICE,
+                  WATER,
+                  MAX_LAYERTYPE
+                };
+                size_t colon_location = value.find(':');
+                WBAssertThrow(colon_location!= std::string::npos,"when chosing Litho1.0 you need to specify a subunit with the colon (:)." );
+                std::string rest_of_string = value.substr(colon_location); // TODO: improve and make more robust, look at gwb-dat example
+                while ((!rest_of_string.empty()) && (rest_of_string[0] == ' '))
+                  rest_of_string.erase(rest_of_string.begin());
+                while ((!rest_of_string.empty()) && (rest_of_string[0] == ':'))
+                  rest_of_string.erase(rest_of_string.begin());
+                while ((!rest_of_string.empty()) && (rest_of_string[0] == ' '))
+                  rest_of_string.erase(rest_of_string.begin());
+                while ((!rest_of_string.empty()) && (rest_of_string[rest_of_string.size() - 1] == ' '))
+                  rest_of_string.erase(rest_of_string.end() - 1);
+                LayerType layer_type = LayerType::MAX_LAYERTYPE;
+                if (rest_of_string == "lithosphere")
+                  {
+                    layer_type = LayerType::LITHOSPHERE;
+                  }
+                else if (rest_of_string == "crust 3")
+                  {
+                    layer_type = LayerType::CRUST3;
+                  }
+                else if (rest_of_string == "crust 2")
+                  {
+                    layer_type = LayerType::CRUST2;
+                  }
+                else if (rest_of_string == "crust 1")
+                  {
+                    layer_type = LayerType::CRUST1;
+                  }
+                else if (rest_of_string == "sediment 3")
+                  {
+                    layer_type = LayerType::SEDIMENT3;
+                  }
+                else if (rest_of_string == "sediment 2")
+                  {
+                    layer_type = LayerType::SEDIMENT2;
+                  }
+                else if (rest_of_string == "sediment 1")
+                  {
+                    layer_type = LayerType::SEDIMENT1;
+                  }
+                else if (rest_of_string == "ice")
+                  {
+                    layer_type = LayerType::ICE;
+                  }
+                else if (rest_of_string == "water")
+                  {
+                    layer_type = LayerType::WATER;
+                  }
+                WBAssertThrow(layer_type != LayerType::MAX_LAYERTYPE, "Could not find litho1.0 layer type " << rest_of_string);
+                std::string litho_1_0_path = "";
+                Value *value_pointer = Pointer("/Litho1.0 path").Get(parameters);
+                WBAssertThrow(value_pointer != nullptr, "You are using the Litho1.0 dataset, but you have not specified the path to the dataset. Please add below \"version\": \"Litho1.0 path\". ");
+                try
+                  {
+                    litho_1_0_path = value_pointer->GetString();
+                  }
+                catch (...)
+                  {
+                    WBAssertThrow(false, "Could not convert values of /Litho1.0 path into a String. ");
+                  }
+                //std::cout << "flag 160" << std::endl;
+                //WBAssertThrow(false,"Got it! path = " << litho_1_0_path);
+                constexpr int n1 = 40962;
+                constexpr int n_layers = 9;
+                float location_data[40962 * 2];
+                float depth_data[40962 * 9];
+                char tessfile[200];
+                sprintf(tessfile, "%s/litho_depth_data.bin", litho_1_0_path.c_str());
+                FILE *fptrb = nullptr;
+                if ((fptrb = fopen(tessfile, "rb")) == nullptr)
+                  {
+                    WBAssertThrow(false, "ERROR: Could not open file " << tessfile);
+                  }
+                fread(&depth_data, sizeof(float), n1 * n_layers, fptrb);
+                fclose(fptrb);
+
+                sprintf(tessfile, "%s/Icosahedron_Level7_LatLon_mod.txt",
+                        litho_1_0_path.c_str());
+                FILE *fp = nullptr;
+                if ((fp = fopen(tessfile, "r")) == nullptr)
+                  {
+                    WBAssertThrow(false, "ERROR: Could not open file " << tessfile);
+                  }
+
+                // fread(&depth_data, sizeof(float), n1 * n_layers, fp);
+                //std::cout << std::endl << litho_1_0_path << std::endl;
+                int counter_location_data = 0;
+                float latitude, glatitude, longitude = 0;
+                while (fscanf(fp, "%f %f %f", &latitude, &glatitude, &longitude) != EOF)
+                  {
+                    const size_t global_index = (counter_location_data/2)*n_layers+layer_type;
+                    result.first.emplace_back(depth_data[global_index]);
+                    result.second.emplace_back(longitude * Consts::PI / 180.);
+                    result.second.emplace_back(latitude * Consts::PI / 180.);
+                    std::cout << counter_location_data << "long:lat:depth = " << longitude << ":"
+                              << latitude << ", " << ", "
+                              << ":" << depth_data[global_index] << " (node:global index:local index = " << counter_location_data/2 << ":" << global_index << ":" << layer_type << ")"<< std::endl;
+
+                    counter_location_data += 2;
+                  }
+
+                fclose(fp);
+                // float depth_data[40962 * 9];
+              }
+            else
+              {
+                std::cout << "could not find \"" << value << "\"" << std::endl;
+                assert(false);
+              }
+            assert(result.first.size() > 0);
+          }
         else
           {
             // case 2: there one value, not an array
@@ -711,7 +862,8 @@ namespace WorldBuilder
         // there is no user defined value. Case one: return the default value and no points
         result.first.emplace_back(default_value);
       }
-
+    assert(result.first.size() > 0);
+    WBAssertThrow(result.first.size() > 0, "Internal error: no entries returned when fetching value at points.");
     return result;
   }
 
diff --git a/source/world_builder/types/one_of.cc b/source/world_builder/types/one_of.cc
@@ -33,6 +33,17 @@ namespace WorldBuilder
 
     }
 
+    OneOf::OneOf(const Interface &type_1,
+                 const Interface &type_2,
+                 const Interface &type_3)
+    {
+      this->type_name = Types::type::OneOf;
+      inner_types_ptr.emplace_back(type_1.clone());
+      inner_types_ptr.emplace_back(type_2.clone());
+      inner_types_ptr.emplace_back(type_3.clone());
+
+    }
+
     OneOf::OneOf(OneOf const &other)
     {
       this->type_name = Types::type::OneOf;
diff --git a/source/world_builder/world.cc b/source/world_builder/world.cc
@@ -113,6 +113,7 @@ namespace WorldBuilder
       prm.declare_entry("$schema", Types::String(""),"The optional filename or https address to a JSON schema file");
 
       prm.declare_entry("cross section", Types::Array(Types::Point<2>(),2,2),"This is an array of two points along where the cross section is taken");
+      prm.declare_entry("Litho1.0 path", Types::String(""), "Path to the litho1.0 files");
 
       prm.declare_entry("potential mantle temperature", Types::Double(1600),
                         "The potential temperature of the mantle at the surface in Kelvin.");
@@ -197,6 +198,7 @@ namespace WorldBuilder
     prm.get_unique_pointers<Features::Interface>("features",prm.features);
 
     const bool set_cross_section = prm.check_entry("cross section");
+    litho_1_0_path = prm.get<std::string>("Litho1.0 path");
 
     const CoordinateSystem coordinate_system = prm.coordinate_system->natural_coordinate_system();
 
