Various bug fixes and improvements

.github/scripts/gen_perf_comment.py

@@ -3,7 +3,7 @@
 from scipy.stats import t
 import html, re
 
-THRESH = 0.03
+THRESH = 0.05
 ALPHA = 0.05
 
 

README.md

@@ -177,6 +177,7 @@ For Windows users, we recommend using our program through Windows Subsystem for
 - If running one of the commands starting with `.venv/bin/cmake` gives you an error, it is likely that a dependency was not installed correctly. Rerun the appropriate commands above to install the required dependencies and try again. If it still fails, make sure you have the virtual environment activated by running `source .venv/bin/activate` in your terminal, and then try again.
 - If you get an error which mentions permission issues, try running the command that gave you the error with `sudo` prefixed. Alternatively, you may follow the next instruction.
 - If you still get permission issues or VScode's `CMake: Install` does not work, make sure you own the relevant directories (i.e. `/usr/local/bin` and the working directory). You may assign ownership to your account with `sudo chown -R $(whoami) .`, replacing `.` with the directory of choice.
+- Use `cartogram --help` to see all available options. If you are not getting the intended output, review the appropriate flags to check whether any of them are relevant to your use case. The default behavior may not always be what you expect, and explicitly specifying a flag can improve performance. For example, using the `--id` flag will speed up the program; otherwise, it will attempt to automatically detect the identifying property, which may take longer.
 
 ### Benchmarking
 
@@ -218,7 +219,7 @@ To push changes to production, please follow the the instructions on [go-cart-io
 
 ### Contributing
 
-Contributions are highly encouraged! Please feel free to take a stab at any at any of the open issues and send in a pull request. If you need help getting setup or more guidance contributing, please @ any of the main contributors (@adisidev, @nihalzp, @mgastner) under any of the open issues (or after creating your own issue), and we'll be happy to guide you!
+Contributions are highly encouraged! Please feel free to take a stab at any at any of the open issues and send in a pull request. If you need help getting setup or more guidance contributing, please @ any of the main contributors (@nihalzp, @adisidev, @mgastner) under any of the open issues (or after creating your own issue), and we'll be happy to guide you!
 
 If you'd like to contribute to the project, please run our tests after you make any changes.
 

include/cartogram_info.hpp

@@ -18,6 +18,7 @@ class CartogramInfo
   std::string map_name_;
   std::map<std::string, std::vector<std::string>> unique_properties_map_;
   int id_col_;
+  bool is_projected_;
 
   // TODO: We assume that either all external rings are counterclockwise or
   //       all are clockwise. This dichotomy covers most geospatial boundary

include/constants.hpp

@@ -65,7 +65,4 @@ constexpr double xi_sq = 4.0;
 // ly
 constexpr unsigned int plotted_cell_length = 8;
 
-// String to identify our custom CRS
-constexpr const char *custom_crs = "EPSG:cartesian";
-
 #endif  // CONST_HPP_

src/cartogram_info/cartogram_info.cpp

@@ -10,6 +10,7 @@ CartogramInfo::CartogramInfo(const Arguments args) : args_(args)
   is_world_map_ = args_.world;
   timer.start("Total time");
   crs_ = "+proj=longlat";
+  is_projected_ = false;
 
   if (!args.visual_file_name.empty()) {
 
@@ -322,12 +323,13 @@ void CartogramInfo::rescale_insets()
 
   // Iterate over insets and normalize areas
   for (InsetState &inset_state : inset_states_) {
-    inset_state.normalize_inset_area(cart_initial_total_target_area());
+    inset_state.normalize_inset_area(
+      default_long_grid_length * default_long_grid_length);
 
     // Rescale copy of original map too
     if (args_.redirect_exports_to_stdout) {
       inset_state.normalize_inset_area(
-        cart_initial_total_target_area(),
+        default_long_grid_length * default_long_grid_length,
         false,
         true);
     }

src/cartogram_info/read_geojson.cpp

@@ -202,14 +202,33 @@ static nlohmann::json load_geojson(const std::string &geometry_file_name)
 }
 
 // Read coordinate reference system if it is included in the GeoJSON
-static void extract_crs(const nlohmann::json &j, std::string &crs)
+static void extract_crs(
+  const nlohmann::json &j,
+  std::string &crs,
+  const bool is_projected)
 {
   if (
     j.contains("crs") && j["crs"].contains("properties") &&
     j["crs"]["properties"].contains("name")) {
     crs = j["crs"]["properties"]["name"];
+    std::cerr << "Coordinate reference system found: " << crs << std::endl;
+    return;
+  }
+  std::cerr << "No `crs` field found." << std::endl;
+  if (!is_projected) {
+    std::cerr << "Default coordinate reference system assumed: " << crs
+              << std::endl;
   }
-  std::cerr << "Coordinate reference system: " << crs << std::endl;
+}
+
+static bool is_projected(const nlohmann::json &j)
+{
+  if (
+    j.contains("properties") && j["properties"].contains("projected") &&
+    j["properties"]["projected"].is_boolean()) {
+    return j["properties"]["projected"];
+  }
+  return false;
 }
 
 // Extract unique properties from GeoJSON and return them as a map
@@ -402,14 +421,16 @@ void CartogramInfo::read_geojson()
     std::exit(19);
   }
 
-  extract_crs(j, crs_);
-  // Skip projection, this is an output from our program
-  if (crs_ == custom_crs) {
-    std::cerr << "WARNING: " << custom_crs << " detected. "
+  // If already projected, skip projection
+  if (is_projected(j)) {
+    is_projected_ = true;
+    std::cerr << "WARNING: `projected=true` property detected. "
               << "Applying --skip_projection flag." << std::endl;
     args_.skip_projection = true;
   }
 
+  extract_crs(j, crs_, is_projected_);
+
   // If args_.id_col is specified, use that as the sole unique property
   if (args_.id_col) {
 

src/cartogram_info/shift_insets_to_position.cpp

@@ -7,8 +7,8 @@ void CartogramInfo::reposition_insets(bool output_to_stdout)
   // Warn user about repositoning insets with `--skip_projection` flag
   if (args_.skip_projection && n_insets() > 1) {
     std::cerr << "WARNING: Trying to repostion insets with ";
-    if (crs_ == custom_crs) {
-      std::cerr << "custom coordinate reference system " << custom_crs << ". ";
+    if (is_projected_) {
+      std::cerr << "original input map already projected. ";
     } else {
       std::cerr << "`--skip_projection` flag present. ";
     }
@@ -31,31 +31,10 @@ void CartogramInfo::reposition_insets(bool output_to_stdout)
     bboxes.at(inset_pos) = inset_state.bbox(output_to_stdout);
   }
 
-  // Calculate the width and height of all positioned insets without spacing
-  double width = bboxes.at("C").xmax() - bboxes.at("C").xmin() +
-                 bboxes.at("L").xmax() - bboxes.at("L").xmin() +
-                 bboxes.at("R").xmax() - bboxes.at("R").xmin();
-
-  // Considering edge cases where the width of the insets named "T" or "B"
-  // might be greater than the width of "C", "L", "R" insets combined
-  width = std::max({
-    bboxes.at("T").xmax() - bboxes.at("T").xmin(),  // width of inset T
-    bboxes.at("B").xmax() - bboxes.at("B").xmin(),  // width of inset B
-    width  // width of inset C + L + R
-  });
-
-  // Similarly for height instead of width
-  double height = bboxes.at("C").ymax() - bboxes.at("C").ymin() +
-                  bboxes.at("T").ymax() - bboxes.at("T").ymin() +
-                  bboxes.at("B").ymax() - bboxes.at("B").ymin();
-  height = std::max({
-    bboxes.at("R").ymax() - bboxes.at("R").ymin(),  // height of inset R
-    bboxes.at("L").ymax() - bboxes.at("L").ymin(),  // height of inset L
-    height  // height of inset C + T + B
-  });
+  const double height_C = bboxes.at("C").ymax() - bboxes.at("C").ymin();
+  const double width_C = bboxes.at("C").xmax() - bboxes.at("C").xmin();
 
   // Spacing between insets
-  const double inset_spacing = std::max(width, height) * inset_spacing_factor;
   for (InsetState &inset_state : inset_states_) {
     std::string inset_pos = inset_state.pos();
 
@@ -70,25 +49,43 @@ void CartogramInfo::reposition_insets(bool output_to_stdout)
       x = std::max(
         {bboxes.at("C").xmax(), bboxes.at("B").xmax(), bboxes.at("T").xmax()});
       x += bboxes.at("R").xmax();
+
+      const double width_R = bboxes.at("R").xmax() - bboxes.at("R").xmin();
+      const double inset_spacing = (width_C + width_R) * inset_spacing_factor;
+
       x += inset_spacing;
     } else if (pos == "L") {
       x = std::min(
         {bboxes.at("C").xmin(), bboxes.at("B").xmin(), bboxes.at("T").xmin()});
 
       // At "L", xmin is negative and lies in the 2nd and 3rd quadrant
       x += bboxes.at("L").xmin();
+
+      const double width_L = bboxes.at("L").xmax() - bboxes.at("L").xmin();
+      const double inset_spacing = (width_C + width_L) * inset_spacing_factor;
+
       x -= inset_spacing;
     } else if (pos == "T") {
       y = std::max(
         {bboxes.at("C").ymax(), bboxes.at("R").ymax(), bboxes.at("L").ymax()});
       y += bboxes.at("T").ymax();
+
+      const double height_T = bboxes.at("T").ymax() - bboxes.at("T").ymin();
+      const double inset_spacing =
+        (height_C + height_T) * inset_spacing_factor;
+
       y += inset_spacing;
     } else if (pos == "B") {
       y = std::min(
         {bboxes.at("C").ymin(), bboxes.at("R").ymin(), bboxes.at("L").ymin()});
 
       // At "B", ymin is negative and lies in the 3rd and 4th quadrant
       y += bboxes.at("B").ymin();
+
+      const double height_B = bboxes.at("B").ymax() - bboxes.at("B").ymin();
+      const double inset_spacing =
+        (height_C + height_B) * inset_spacing_factor;
+
       y -= inset_spacing;
     }
 

src/cartogram_info/write_geojson.cpp

@@ -192,13 +192,13 @@ void CartogramInfo::json_to_geojson(
       add_dividers_to_geojson(container[(container.size() - 1)]);
   }
 
-  // Add/replace CRS to custom_crs
-  new_json["crs"] = {{"type", "name"}, {"properties", {{"name", custom_crs}}}};
-
   new_json["properties"]["note"] =
     "Created using cartogram-cpp / go-cart.io with custom projection, not in "
     "EPSG:4326";
 
+  // Write that the map is projected
+  new_json["properties"]["projected"] = true;
+
   // Iterate over GeoDivs and gd_ids in the container. The index
   // container.size()-2 is reserved for the bounding box, and the index
   // container.size()-1 is reserved for the divider lines. Thus, we must

src/inset_state/albers_projection.cpp

@@ -140,7 +140,7 @@ void InsetState::apply_albers_projection()
   const double phi_1 = 0.5 * (phi_0 + max_lat);
   const double phi_2 = 0.5 * (phi_0 + min_lat);
 
-  static thread_local AlbersProjector proj(lambda_0, phi_0, phi_1, phi_2);
+  const AlbersProjector proj(lambda_0, phi_0, phi_1, phi_2);
 
   transform_points([&](const Point &q) {
     return proj(q);

src/inset_state/integrate.cpp

@@ -56,15 +56,15 @@ void InsetState::prepare_for_integration()
   }
 
   // Store initial inset area to calculate area drift,
-  // set area errors based on this initial_area
   store_initial_area();
-  set_area_errors();
 
   // Store initial target area to normalize inset areas
   store_initial_target_area();
 
   // Normalize total target area to be equal to initial area
   normalize_target_area();
+
+  set_area_errors();
 }
 
 void InsetState::cleanup_after_integration()
@@ -147,6 +147,7 @@ bool InsetState::continue_integrating() const
 
 void InsetState::integrate(ProgressTracker &progress_tracker)
 {
+  std::cerr << std::endl << "Integrating inset " << pos_ << std::endl;
 
   timer.start(inset_name_);
 

src/misc/parse_arguments.cpp

@@ -8,7 +8,7 @@ Arguments parse_arguments(const int argc, const char *argv[])
 
   // Create parser for arguments using argparse.
   // From https://github.com/p-ranav/argparse
-  argparse::ArgumentParser arguments("./cartogram", "25.8");
+  argparse::ArgumentParser arguments("./cartogram", "25.9");
 
   // Positional argument accepting geometry file (GeoJSON, JSON) as input
   arguments.add_argument("geometry_file")