From 9c740cbf1bf97f90d2cc35d10d84eecfb203b25d Mon Sep 17 00:00:00 2001
From: Luca Delucchi <lucadeluge@gmail.com>
Date: Tue, 18 Feb 2025 15:49:21 +0100
Subject: [PATCH] added capability to generate vector polygon of the raster
 extent

---
 opendm/orthophoto.py | 48 +++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 47 insertions(+), 1 deletion(-)

diff --git a/opendm/orthophoto.py b/opendm/orthophoto.py
index 319a504c..7014f415 100644
--- a/opendm/orthophoto.py
+++ b/opendm/orthophoto.py
@@ -14,6 +14,7 @@ from opendm import io
 from opendm.tiles.tiler import generate_orthophoto_tiles
 from opendm.cogeo import convert_to_cogeo
 from osgeo import gdal
+from osgeo import ogr
 
 
 def get_orthophoto_vars(args):
@@ -84,7 +85,50 @@ def generate_kmz(orthophoto_file, output_file=None, outsize=None):
 
     system.run('gdal_translate -of KMLSUPEROVERLAY -co FORMAT=PNG "%s" "%s" %s '
                '--config GDAL_CACHEMAX %s%% ' % (orthophoto_file, output_file, bandparam, get_max_memory()))    
-    
+
+def generate_extent_polygon(orthophoto_file, output_file=None):
+    """Function to return the orthophoto extent as a polygon into a gpkg file
+
+    Args:
+        orthophoto_file (str): the path to orthophoto file
+        output_file (str, optional): the path to the gpkg file. Defaults to None.
+    """
+    if output_file is None:
+        base, ext = os.path.splitext(orthophoto_file)
+        output_file = base + '_extent.gpkg'
+    # read geotiff file
+    gtif = gdal.Open(orthophoto_file)
+    srs =  gtif.GetSpatialRef()
+    geoTransform = gtif.GetGeoTransform()
+    # calculate the coordinates
+    minx = geoTransform[0]
+    maxy = geoTransform[3]
+    maxx = minx + geoTransform[1] * gtif.RasterXSize
+    miny = maxy + geoTransform[5] * gtif.RasterYSize
+    # create polygon in wkt format
+    poly_wkt = "POLYGON ((%s %s, %s %s, %s %s, %s %s, %s %s))" % (minx, miny, minx, maxy, maxx, maxy, maxx, miny, minx, miny)
+    # set up the shapefile driver
+    driver = ogr.GetDriverByName("GPKG")
+
+    # create the data source
+    ds = driver.CreateDataSource(output_file)
+
+    # create one layer
+    layer = ds.CreateLayer("extent", srs, ogr.wkbPolygon)
+    layer.CreateField(ogr.FieldDefn("id", ogr.OFTInteger))
+
+    # create the feature and set values
+    featureDefn = layer.GetLayerDefn()
+    feature = ogr.Feature(featureDefn)
+    feature.SetGeometry(ogr.CreateGeometryFromWkt(poly_wkt))
+    feature.SetField("id", 1)
+    # add feature to layer
+    layer.CreateFeature(feature)
+    # save and close everything
+    feature = None
+    ds = None
+    return True
+
 def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_tiles_dir, resolution):
     if args.crop > 0 or args.boundary:
         Cropper.crop(bounds_file_path, orthophoto_file, get_orthophoto_vars(args), keep_original=not args.optimize_disk_space, warp_options=['-dstalpha'])
@@ -104,6 +148,8 @@ def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_ti
     if args.cog:
         convert_to_cogeo(orthophoto_file, max_workers=args.max_concurrency, compression=args.orthophoto_compression)
 
+    generate_extent_polygon(orthophoto_file)
+
 def compute_mask_raster(input_raster, vector_mask, output_raster, blend_distance=20, only_max_coords_feature=False):
     if not os.path.exists(input_raster):
         log.ODM_WARNING("Cannot mask raster, %s does not exist" % input_raster)