Got area selection calculations working

requirements.txt

@@ -4,3 +4,5 @@ geojson==3.1.0
 ipyleaflet==0.18.2
 ipywidgets==8.1.2
 leafmap==0.31.9
+rioxarray
+geopandas

toy_app/pages/03_cog.py

@@ -1,5 +1,6 @@
 from ipyleaflet import Polygon, WidgetControl
 from ipywidgets import widgets
+from ipywidgets.widgets.widget import jsondumps
 import leafmap
 import logging
 import solara
@@ -17,6 +18,11 @@ log = logging.getLogger(__name__)
 libya_fire_url = (
     "https://github.com/opengeos/data/releases/download/raster/Libya-2023-07-01.tif"
 )
+richness = "https://data.source.coop/cboettig/mobi/species-richness-all/SpeciesRichness_All.tif"
+deforest = "https://data.source.coop/vizzuality/lg-land-carbon-data/deforest_carbon_100m_cog.tif"
+irrecoverable = (
+    "https://data.source.coop/cboettig/carbon/cogs/irrecoverable_c_total_2018.tif"
+)
 libya_fire_name = "Libya Fire"
 libya_layer_visible = solara.reactive("True")
 test_poly = {
@@ -42,6 +48,8 @@ test_poly = {
         ],
     },
 }
+geometry = solara.reactive(jsondumps(test_poly))
+lost_carbon_count = solara.reactive(0)
 
 
 class Map(leafmap.Map):
@@ -64,6 +72,13 @@ class Map(leafmap.Map):
         self.add_cog_layer(
             url=libya_fire_url, name=libya_fire_name
         )  # Just display the cog for now, the layer toggle isn't working as expected
+        self.add_cog_layer(
+            url=irrecoverable,
+            name="irrecoverable c",
+            palette="reds",
+            zoom_to_layer=False,
+            opacity=0.8,
+        )
         # self.toggle_libya_layer()
         # user_poly = map_utils.get_user_polygon(self)
         # if user_poly:
@@ -71,82 +86,39 @@ class Map(leafmap.Map):
         # self.add_click_detector_widget()
         self.add_button()
 
-    def add_click_detector_widget(self, position="bottomright"):
-        output = widgets.Output()
-        control = WidgetControl(widget=output, position=position)
-        self.add(control)
-
-        def update_user_click(**kwargs):
-            with output:
-                print(kwargs)
-
-        self.on_interaction(update_user_click)
-
     def add_button(self):
-        button = widgets.Button(description="Get ROI")
+        button = widgets.Button(description="Calculate")
         output = widgets.Output()
 
         def on_button_click(_):
             if self.user_roi is not None:
-                with output:
-                    output.clear_output()
-                    print(map_utils.calculate_area(self.user_roi))
+                geometry.set(self.user_roi)
+                selected_polygon = map_utils.read_polygon(self.user_roi)
+                selected_on_carbon = map_utils.extract_geom(
+                    selected_polygon, irrecoverable
+                ).fillna(0)
+                area = round(float(map_utils.area_hectares(selected_polygon)))
+                carbon_total = round(float(selected_on_carbon.mean()) * area)
+                lost_carbon_count.set(carbon_total)
 
         button.on_click(on_button_click)
 
         widget = widgets.VBox([button, output])
         self.add_widget(widget)
 
-    def toggle_libya_layer(self):
-        if libya_layer_visible.value == "True":
-            # The layer should be visible, if you don't find it, add it
-            log.warning("Layer should be visible, checking it needs to be added")
-            for layer in self.layers:
-                l_name = layer.name
-                if l_name == libya_fire_name:
-                    return
-            else:
-                log.warning("Layer not found, adding new cog layer")
-                self.add_cog_layer(url=libya_fire_url, name=libya_fire_name)
-        else:
-            # The layer should not be visible, if you find it, remove it
-            log.warning("Layer should not be visible, checking to see if it's there")
-            for layer in self.layers:
-                l_name = layer.name
-                if l_name != libya_fire_name:
-                    continue
-                log.warning("Layer found, removing it")
-                self.remove_layer(l_name)
-
 
 zoom = solara.reactive(2)
 center = solara.reactive((20, 0))
-drawing = solara.reactive(Polygon())
 
 
 @component
 def Page():
-    l_map = Map()
-
-    # def toggle_map():
-    #     log.warning("Attempting to toggle the layer")
-    #     if libya_layer_visible.value == "True":
-    #         libya_layer_visible.set("False")
-    #     else:
-    #         libya_layer_visible.set("True")
-    #     log.warning(f"New layer visibility is set to {libya_layer_visible.value}")
-    def calc_area():
-        log.info("Trying to calculate the area")
-        # user_polygon = l_map.user_roi
-        user_polygon = drawing.value
-        log.warning(f"The polygon is: {user_polygon}")
-        log.warning(f"The user_roi is: {l_map.user_rois}")
-
-        print("Not yet implemented")
-
     with solara.Column(style={"min_width": "500px"}):
-        # Map.element(  # type: ignore
-        l_map.element(  # type: ignore
+        with solara.Row():
+            solara.Text(
+                f"Average Carbon Lost within the selected area in 2018: {lost_carbon_count.value:,} Tonnes"
+            )
+        Map.element(  # type: ignore
             zoom=zoom.value,
             on_zoom=zoom.set,
             center=center.value,
@@ -154,14 +126,9 @@ def Page():
             scroll_wheel_zoom=True,
             toolbar_ctrl=False,
             draw_control=True,
-            on_draw=drawing.set,
-            # data_ctrl=False,
             height="780px",
         )
         with solara.Row():
             solara.Text(f"Center: {center.value}")
             solara.Text(f"Zoom: {zoom.value}")
-
-    with solara.Row():
-        solara.Button(label="Calculate Polygon Area", on_click=calc_area)
-        # solara.Button(label="Toggle Fire Image", on_click=toggle_map)  # Disable toggle for now since it's not working as expected
+            solara.Text(f"UserPoly: {geometry.value}")

toy_app/src/map_utils.py

@@ -7,6 +7,9 @@ import typing
 import leafmap
 import geojson
 from area import area as area
+import rioxarray
+import json
+import geopandas as gpd
 
 
 log = logging.getLogger(__name__)
@@ -21,3 +24,22 @@ def calculate_area(poly: geojson.GeoJSON) -> float:
         raise TypeError("The given GeoJSON did not contain a Polygon as was expected")
 
     return output
+
+
+def extract_geom(gdf, cog):
+    x = rioxarray.open_rasterio("/vsicurl/" + cog, masked=True).rio.clip(
+        gdf.geometry.values, gdf.crs, from_disk=True
+    )
+    return x
+
+
+def read_polygon(polygon):
+    geojson_str = json.dumps(polygon)
+    gdf = gpd.read_file(geojson_str, driver="GeoJSON")
+    gdf.set_crs("epsg:4326")
+    return gdf
+
+
+def area_hectares(gdf):
+    area = gdf.to_crs("EPSG:9822").area / 10000.0
+    return area