Improve get_max_pixel_coords function

leafmap/common.py

@@ -14870,7 +14870,15 @@ def read_geojson(data: str, **kwargs: Any) -> Dict[str, Any]:
     return requests.get(data, **kwargs).json()
 
 
-def get_max_pixel_coords(geotiff_path, band_idx=1, roi=None, dst_crs="EPSG:4326"):
+def get_max_pixel_coords(
+    geotiff_path,
+    band_idx=1,
+    roi=None,
+    dst_crs="EPSG:4326",
+    output=None,
+    return_gdf=True,
+    **kwargs,
+):
     """
     Find the geographic coordinates of the maximum pixel value in a GeoTIFF.
 
@@ -14879,6 +14887,8 @@ def get_max_pixel_coords(geotiff_path, band_idx=1, roi=None, dst_crs="EPSG:4326"
         band_idx (int): Band index to use (default is 1).
         roi (str): Path to a vector dataset containing the region of interest (default is None).
         dst_crs (str): Desired output coordinate system in EPSG format (e.g., "EPSG:4326").
+        output (str): Path to save the output GeoDataFrame (default is None).
+        return_gdf (bool): Whether to return a GeoDataFrame (default is True).
 
     Returns:
         dict: Maximum pixel value and its geographic coordinates in the specified CRS.
@@ -14893,7 +14903,16 @@ def get_max_pixel_coords(geotiff_path, band_idx=1, roi=None, dst_crs="EPSG:4326"
     with rasterio.open(geotiff_path) as dataset:
         # If ROI is provided, handle potential CRS differences
         if roi:
-            gdf = gpd.read_file(roi)
+            if isinstance(roi, str):
+                gdf = gpd.read_file(roi)
+            elif isinstance(roi, gpd.GeoDataFrame):
+                gdf = roi
+            elif isinstance(roi, dict):
+                gdf = gpd.GeoDataFrame.from_features([roi])
+            else:
+                raise ValueError(
+                    "Invalid ROI input. Must be a file path or a GeoDataFrame."
+                )
             roi_geojson = gdf.__geo_interface__
 
             # Reproject ROI to match the raster's CRS if necessary
@@ -14931,4 +14950,19 @@ def get_max_pixel_coords(geotiff_path, band_idx=1, roi=None, dst_crs="EPSG:4326"
         src_crs = dataset.crs
         x, y = transform(src_crs, dst_crs, [original_coords[0]], [original_coords[1]])
 
-    return {"max_value": max_value, "coordinates": (x[0], y[0]), "crs": dst_crs}
+        if return_gdf:
+            x_coords = [x[0]]
+            y_coords = [y[0]]
+            # Create a DataFrame
+            df = pd.DataFrame({"x": x_coords, "y": y_coords})
+
+            # Convert the DataFrame to a GeoDataFrame
+            gdf = gpd.GeoDataFrame(
+                df, geometry=gpd.points_from_xy(df.x, df.y), crs=dst_crs
+            )
+
+            if output:
+                gdf.to_file(output, **kwargs)
+
+        else:
+            return {"max_value": max_value, "coordinates": (x[0], y[0]), "crs": dst_crs}