From 9c740cbf1bf97f90d2cc35d10d84eecfb203b25d Mon Sep 17 00:00:00 2001 From: Luca Delucchi Date: Tue, 18 Feb 2025 15:49:21 +0100 Subject: [PATCH 1/4] 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) From 2215e7ab7a7d4ca45f1d086ed9a19bce200d798a Mon Sep 17 00:00:00 2001 From: Luca Delucchi Date: Wed, 19 Feb 2025 22:17:38 +0100 Subject: [PATCH 2/4] added creation also of a dxf file --- opendm/orthophoto.py | 54 +++++++++++++++++++++++++------------------- 1 file changed, 31 insertions(+), 23 deletions(-) diff --git a/opendm/orthophoto.py b/opendm/orthophoto.py index 7014f415..147fef9a 100644 --- a/opendm/orthophoto.py +++ b/opendm/orthophoto.py @@ -93,9 +93,34 @@ def generate_extent_polygon(orthophoto_file, output_file=None): 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' + def _create_vector(ortho_file, poly, format, output=None): + if output is None: + base, ext = os.path.splitext(ortho_file) + output_file = base + '_extent.' + format.lower() + # set up the shapefile driver + driver = ogr.GetDriverByName(format) + + # create the data source + ds = driver.CreateDataSource(output_file) + + # create one layer + layer = ds.CreateLayer("extent", srs, ogr.wkbPolygon) + if format != "DXF": + 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)) + if format != "DXF": + feature.SetField("id", 1) + # add feature to layer + layer.CreateFeature(feature) + # save and close everything + feature = None + ds = None + return True + # read geotiff file gtif = gdal.Open(orthophoto_file) srs = gtif.GetSpatialRef() @@ -107,26 +132,9 @@ def generate_extent_polygon(orthophoto_file, output_file=None): 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 + # create vector file + _create_vector(orthophoto_file, poly_wkt, "GPKG", output_file) + _create_vector(orthophoto_file, poly_wkt, "DXF", output_file) return True def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_tiles_dir, resolution): From ece421faec223772f903e564270fc3279d673928 Mon Sep 17 00:00:00 2001 From: Luca Delucchi Date: Thu, 27 Feb 2025 06:08:40 +0100 Subject: [PATCH 3/4] keep only DXF to support AutoCAD users to see correctly the output geotiff --- opendm/orthophoto.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/opendm/orthophoto.py b/opendm/orthophoto.py index 147fef9a..923cc491 100644 --- a/opendm/orthophoto.py +++ b/opendm/orthophoto.py @@ -133,7 +133,9 @@ def generate_extent_polygon(orthophoto_file, output_file=None): # 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) # create vector file - _create_vector(orthophoto_file, poly_wkt, "GPKG", output_file) + # just the DXF to support AutoCAD users + # to load the geotiff raster correctly. + # _create_vector(orthophoto_file, poly_wkt, "GPKG", output_file) _create_vector(orthophoto_file, poly_wkt, "DXF", output_file) return True From 131b29c3f3897a914e15e4970706ed37d3e95ce6 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Thu, 27 Feb 2025 01:31:06 -0500 Subject: [PATCH 4/4] Add try/catch, minor tweaks --- opendm/orthophoto.py | 40 +++++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/opendm/orthophoto.py b/opendm/orthophoto.py index 923cc491..54b1e9b7 100644 --- a/opendm/orthophoto.py +++ b/opendm/orthophoto.py @@ -91,7 +91,7 @@ def generate_extent_polygon(orthophoto_file, output_file=None): Args: orthophoto_file (str): the path to orthophoto file - output_file (str, optional): the path to the gpkg file. Defaults to None. + output_file (str, optional): the path to the output file. Defaults to None. """ def _create_vector(ortho_file, poly, format, output=None): if output is None: @@ -119,25 +119,27 @@ def generate_extent_polygon(orthophoto_file, output_file=None): # save and close everything feature = None ds = None - return True - # 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) - # create vector file - # just the DXF to support AutoCAD users - # to load the geotiff raster correctly. - # _create_vector(orthophoto_file, poly_wkt, "GPKG", output_file) - _create_vector(orthophoto_file, poly_wkt, "DXF", output_file) - return True + try: + 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) + # create vector file + # just the DXF to support AutoCAD users + # to load the geotiff raster correctly. + # _create_vector(orthophoto_file, poly_wkt, "GPKG", output_file) + _create_vector(orthophoto_file, poly_wkt, "DXF", output_file) + gtif = None + except Exception as e: + log.ODM_WARNING("Cannot create extent layer for %s: %s" % (ortho_file, str(e))) + def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_tiles_dir, resolution): if args.crop > 0 or args.boundary: