Structure from motion using ArcGIS orthomapping Python Api

Share This Post

In my previous articles on photogrammetry, I introduced structure from motion photogrammetry, explained block adjustment, and demonstrated how to generate photogrammetry products for urban monitoring and biodiversity conservation using ArcGIS PRO’s reality mapping feature. ArcGIS PRO provides a user-friendly interface and powerful tools for preparing and processing UAV, airborne, or satellite data, generating various results, visualizing and analyzing the data, and performing in-depth exploration. However, in some cases, automating the orthomapping or reality mapping process is necessary. This can be achieved using the ArcGIS Python API, which also allows developers to incorporate the tools into their scripts, workflows, or projects.

In this article, I will demonstrate how to use the ArcGIS Orthomapping API to generate orthophotos, Digital Terrain Models (DTM), and Digital Surface Models (DSM) from drone photos. For the purpose of this demonstration, I will use five out of the 28 photos of the Esri Headquarters building, which are used in the reality mapping documentation tutorial. It’s important to note that, at the time of writing this article, orthomapping with the ArcGIS Python API is only possible with ArcGIS Enterprise configured with a raster analysis server.

Environment Setup¶

The first step is to set up the environment. The Orthomapping API requires the ArcGIS API for Python (version >= 1.4.2) to be installed, and the ArcGIS Enterprise must have a raster analysis server with an image server license privilege. The ArcGIS API can be set up using package managers (conda, pip), ArcGIS PRO, Docker, Google Collab, or an offline installation. To use VS Code with an installation of the ArcGIS API on the ArcGIS PRO default environment, you need to clone the environment, navigate to the cloned environment, and copy the address. Then, open VS Code, press Ctrl + Shift + P (or Cmd + Shift + P on macOS) to open the Command Palette, click on the option “Enter interpreter path…” at the bottom of the dropdown menu, and enter the path to your Python environment (something like “C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3_clone\python.exe”) in the input field.

The Orthomapping API workflow can be broadly categorized into three steps: creating an image collection, adjusting images, and generating products. These steps are demonstrated in the following sections.

Step 1 : Creating Image Collection¶

ArcGIS Orthomapping works with image collections. An image collection is created using the arcgis.raster.analytics.create_image_collection() method, which takes multiple image items as input and creates a mosaic dataset in a registered raster store. To get started, let’s import the required modules for this task

In [ ]:

import cv2
import matplotlib.pyplot as plt
import os
import glob
import ntpath
import arcgis
from arcgis.gis import GIS
import exifread
from arcgis.raster.analytics import create_image_collection
from arcgis.raster.orthomapping import compute_sensor_model
from arcgis.raster.orthomapping import generate_report
from arcgis.raster.orthomapping import match_control_points
from arcgis.raster.orthomapping import edit_control_points
from arcgis.raster.orthomapping import query_control_points
from arcgis.raster.orthomapping import compute_seamlines
from arcgis.raster.orthomapping import color_correction
from arcgis.raster.orthomapping import generate_orthomosaic
from arcgis.raster.orthomapping import generate_dem
from arcgis.raster.orthomapping import get_processing_states

Next, we connect to GIS and confirm that we have the necessary support to perform the orthomapping processing.

In [ ]:

portal_url = "https://yourportal.domain.com/portal"
portal_username = "your username"
portal_password = "your password"

gis = GIS(url=portal_url, username=portal_username, password=portal_password)

In [ ]:

print(f'raster anlytics supported: {arcgis.raster.analytics.is_supported(gis)}')
print(f'orthomapping supported: {arcgis.raster.orthomapping.is_supported(gis)}')

raster anlytics supported: True
orthomapping supported: True

Now let’s create a project folder where all our items will be stored.

In [ ]:

# Create a project folder
prj_folder_name = "Esri_Headquaters5"
prj_folder_item = gis.content.create_folder(folder=prj_folder_name, owner=portal_username)

The next cell contains code to upload images to an ArcGIS portal, extract GPS information from the images, set camera properties and essentially define the raster_type_params variable that will serve as an input for the creat_image_collection() function. The code crawls through a folder containing JPEGs, creates an image item template, uploads each image to the project folder, and extracts GPS information (latitude, longitude, and altitude) from each image’s EXIF header file. Information about the camera used to capture the images is specified in the camera_properties’ variable.

In [ ]:

# Upload images to the project folder
image_folder_path = r"\\folder\path\to\images\data_EsriHeadquater5\Images"
image_list = glob.glob(os.path.join(image_folder_path, '*.JPG'))
image_item_list = []
item_prop_template = {"type": "Image"}

for image_full_path in image_list:
    image_name = ntpath.split(image_full_path)[1]
    item_prop_template["title"] = image_name
    item_prop_template["tags"] = image_name
    item_prop_template["description"] = image_name

    image_item = gis.content.add(item_properties=item_prop_template, data=image_full_path, 
                                owner=portal_username, folder=prj_folder_name)
    image_item_list.append(image_item)

# Get Images GPS information

def dms_to_dd(dms):
    degrees = dms[0].num / dms[0].den
    minutes = dms[1].num / dms[1].den / 60
    seconds = dms[2].num / dms[2].den / 3600
    dd = degrees + minutes + seconds
    return dd

def get_image_location(image_path):
    with open(image_path, 'rb') as file:
        # Read the EXIF tags
        tags = exifread.process_file(file, details=False)
        
        # Check if GPS information is available
        if 'GPS GPSLatitude' in tags and 'GPS GPSLongitude' in tags:
            # Extract latitude and longitude values
            lat = tags['GPS GPSLatitude']
            lon = tags['GPS GPSLongitude']
            
            # Extract latitude reference (N/S) and convert to a sign (+/-)
            lat_ref = tags['GPS GPSLatitudeRef']
            lat_sign = 1 if lat_ref.values == 'N' else -1
            
            # Extract longitude reference (E/W) and convert to a sign (+/-)
            lon_ref = tags['GPS GPSLongitudeRef']
            lon_sign = 1 if lon_ref.values == 'E' else -1
            
            # Convert latitude and longitude to decimal degrees
            lat_dd = dms_to_dd(lat.values) * lat_sign
            lon_dd = dms_to_dd(lon.values) * lon_sign
            
            # Extract altitude if available
            if 'GPS GPSAltitude' in tags:
                alt = tags['GPS GPSAltitude']
                altitude = float(alt.values[0].num) / float(alt.values[0].den)
            else:
                altitude = None
            
            # Return image filename, latitude, longitude, and altitude as a list
            return [ntpath.split(image_path)[1], lat_dd, lon_dd, altitude]
    
    # Return None if GPS information is not available or not found
    return None

gps2 = [get_image_location(image) for image in image_list]

# Set camera properties and raster params

camera_properties = {"maker":"DJI","model":"FC6510","focallength":8.8,"columns":5472,"rows":3648,"pixelsize":0.0024}

raster_type_params = {
    "gps": gps2, "cameraProperties": camera_properties, "isAltitudeFlightHeight": "False",
    "averagezdem":{
        "url":"https://elevation3d.arcgis.com/arcgis/rest/services/WorldElevation3D/Terrain3D/ImageServer"
    }
}

raster_type_params

Out[ ]:

{'gps': [['DJI_0006.JPG', 34.05874605555555, -117.19607425000001, 407.935],
  ['DJI_0007.JPG', 34.05866786111111, -117.19616886111112, 407.735],
  ['DJI_0008.JPG', 34.05859125, -117.19626125, 407.735],
  ['DJI_0020.JPG', 34.05867922222222, -117.19652513888889, 407.835],
  ['DJI_0021.JPG', 34.05875980555555, -117.19642794444445, 407.835]],
 'cameraProperties': {'maker': 'DJI',
  'model': 'FC6510',
  'focallength': 8.8,
  'columns': 5472,
  'rows': 3648,
  'pixelsize': 0.0024},
 'isAltitudeFlightHeight': 'False',
 'averagezdem': {'url': 'https://elevation3d.arcgis.com/arcgis/rest/services/WorldElevation3D/Terrain3D/ImageServer'}}

With the necessary preparations complete, we can now call the create_image_collection() function. The function parameters include the image collection name, the raster items, the raster type name, the raster type parameters, and the folder where the image collection item will be created in the ArcGIS portal.

In [ ]:

# Create an image collection Imagery Layer Item

image_collection_name = "EsriHead5ImageCollection"
image_collection_item = create_image_collection(image_collection=image_collection_name,
                                                input_rasters=image_item_list,
                                                raster_type_name="UAV/UAS",
                                                raster_type_params=raster_type_params,
                                                folder=prj_folder_name)

After executing the code, we can observe the created image collection layer.

In [ ]:

image_collection_item.layers[0]

Out[ ]:

Step 2: Image Adjustment¶

Once the image collection layer is created, we need to perform block adjustment, including adding ground control points, and compute seamlines and color correction.

2.1 Compute and apply initial block adjustment¶

The cell below calls the compute_sensor_model() function which performs a bundle block adjustment on the image collection item and applies a frame transform to the images based on the raster type information. The mode parameter is set to ‘Quick’, indicating that the adjustment should be done at reduced image resolution for fast or quick estimation.

In [ ]:

compute_sensor_model(image_collection=image_collection_item, mode='Quick', location_accuracy='Low')

Let’s check the processing report to get information about the quality of the adjusted images.

In [ ]:

# Generate realitymapping report
generate_report(image_collection=image_collection_item, report_format="PDF")

2.2 Add ground control points¶

Adding ground control points is an optional step but recommended for accurate georeferencing of the images. If control points were acquired through field surveying, you need to determine the pixel location of the GCP target in at least one of the images. This can be done using the image information window in ArcGIS PRO or any other convenient method. Once you have the input control points defined, you can use the arcgis.raster.orthomapping.match_control_points() method to find matching tie points on the image items of the image collection. In cases where control points are not surveyed in the field, the arcgis.raster.orthomapping.compute_control_points() function can be used to compute control point sets between the reference image and the mosaicked image items within the image collection.

In [ ]:

input_control_points = [{"pointId": 1,"x": 481870.2391,"y": 3768697.345,"z": 392.046,
                       "spatialReference": {"wkid":26911},
                       "xyAccuracy": "0.017","zAccuracy": "0.057",
                       "imagePointSpatialReference": {"ics":"", "topup":False}, 
                       "imagePoints": [{
                           "imageID": 2,
                           "x": 5332,
                           "y": -1809,
                           "u":5332,
                           "v":-1809
                       },
                       {
                           "imageID": 5,
                           "x": 2175,
                           "y": -1876,
                           "u": 2175,
                           "v": -1876
                       }] 
                     }]

In [ ]:

control_point_sets = match_control_points(image_collection=image_collection_item, 
                                          control_points=input_control_points, 
                                          similarity="High")

In [ ]:

control_point_sets

Out[ ]:

[{'pointID': 1,
  'x': 481870.23913892277,
  'y': 3768697.3450019867,
  'z': 392.0460000000021,
  'spatialReference': {'wkid': 26911, 'latestWkid': 26911},
  'xyAccuracy': 0.017,
  'ZAccuracy': 0.057,
  'imagePointSpatialReference': {'wkid': 32611, 'latestWkid': 32611},
  'imagePoints': [{'imageID': 2,
    'x': 481873.31900000013,
    'y': 3768694.0404000003,
    'u': 5332,
    'v': -1809},
   {'imageID': 5,
    'x': 481873.8250000002,
    'y': 3768693.6229999997,
    'u': 2175,
    'v': -1876},
   {'imageID': 1,
    'x': 481873.33299999963,
    'y': 3768693.7498000003,
    'u': 5337.953758484229,
    'v': -679.599729284758},
   {'imageID': 3,
    'x': 481873.3822999997,
    'y': 3768694.4681,
    'u': 5340.696405420868,
    'v': -2940.4835710592693},
   {'imageID': 4,
    'x': 481873.92080000043,
    'y': 3768693.8454,
    'u': 2204.6979602035144,
    'v': -747.9752159571483}]}]

In [ ]:

control_point_sets=[{"type":2, 'x': 481870.23913892277, 'y': 3768697.3450019867, 'z': 392.0460000000021,
                     "pointId":1,"xyAccuracy":"0.017","zAccuracy":"0.057",
                       'imagePoints': [{'imageID': 2,
                                    'x': 481873.31900000013,
                                    'y': 3768694.0404000003,
                                    'u': 5332,
                                    'v': -1809},
                                {'imageID': 5,
                                    'x': 481873.8250000002,
                                    'y': 3768693.6229999997,
                                    'u': 2175,
                                    'v': -1876},
                                {'imageID': 1,
                                    'x': 481873.33299999963,
                                    'y': 3768693.7498000003,
                                    'u': 5337.953758484229,
                                    'v': -679.599729284758},
                                {'imageID': 3,
                                    'x': 481873.3822999997,
                                    'y': 3768694.4681,
                                    'u': 5340.696405420868,
                                    'v': -2940.4835710592693},
                                {'imageID': 4,
                                    'x': 481873.92080000043,
                                    'y': 3768693.8454,
                                    'u': 2204.6979602035144,
                                    'v': -747.9752159571483}]}]

In [ ]:

edit_control_points(image_collection=image_collection_item, control_points=control_point_sets)

After adding the ground control point sets to the image collection’s control point table, you can rerun the compute sensor model to refine the block adjustment result.

In [ ]:

compute_sensor_model(image_collection=image_collection_item, mode='Refine', location_accuracy='Medium')

You can compare the processing report again to see how the projection error has changed.

In [ ]:

generate_report(image_collection= image_collection_item, report_format="PDF")

2.3 Compute seamlines and color correction¶

The next part of the process involves computing seamlines and applying color correction. Seamlines are used to create a visually seamless mosaic display when working with overlapped images in an image collection. They help eliminate visible seams or discontinuities that may appear when combining multiple images together. By applying seamlines, the images can be blended seamlessly, resulting in a more cohesive and visually appealing mosaic display. Color correction is important for improving the visual appearance, especially for sensors like satellites where cross-strip images may have different acquisition conditions.

In [ ]:

compute_seamlines(image_collection=image_collection_item, seamlines_method="DISPARITY",
                  context={"minRegionSize":100, "blendType":"Both", "blendWidth":None, "blendUnit":"Pixels", 
                           "requestSizeType":"Pixels", "requestSize":1000, "minThinnessRatio":0.07, "maxSilverSize":20})

In [ ]:

color_correction(image_collection=image_collection_item, color_correction_method="DODGING", dodging_surface_type="SECOND_ORDER",
                 context={"skipRows": 10, "skipCols": 10, "reCalculateStats": "OVERWRITE"})

Step 3: Generating products¶

After the adjustment process, you can generate orthomosaic, DTM, and DSM products. The code snippets below demonstrate how to generate these products.

3.1 Generate orthomosaic¶

In [ ]:

ortho_mosaic_name = "EsriHeadQ5_OrthoMosaic"
ortho_mosaic_item = generate_orthomosaic(image_collection=image_collection_item, out_ortho=ortho_mosaic_name,
                                         regen_seamlines=True, recompute_color_correction=True, folder=prj_folder_name)

In [ ]:

ortho_mosaic_item.url

In [ ]:

ortho_mosaic_item.layers[0]

Out[ ]:

3.2 Generate DEM¶

In [ ]:

dtm_name = "EsriHeadQ5DTM"
dtm_item=generate_dem(image_collection=image_collection_item, out_dem=dtm_name, cell_size=0.1314245599999937,
                     surface_type="DTM", matching_method="ETM",
                     context={"maxObjectSize":15,"minAngle":10,"maxAngle":70,"minOverlap":0.6,"maxGSDDif":2,
                              "numImagePairs":8,"adjQualityThreshold":0.2,"method":"TRIANGULATION",
                              "smoothingMethod":"GAUSS5x5","applyToOrtho":True,"regenPointCloud":False},
                     folder=prj_folder_name)

In [ ]:

dtm_item.url

In [ ]:

dtm_item.layers[0]

Out[ ]:

In [ ]:

dsm_name = "EsriHeadQ5DSM"
dsm_item=generate_dem(image_collection=image_collection_item, out_dem=dsm_name, cell_size=0.1314245599999937,
                     surface_type="DSM", matching_method="ETM",
                     context={"maxObjectSize":15,"minAngle":10,"maxAngle":70,"minOverlap":0.6,"maxGSDDif":2,
                              "numImagePairs":8,"adjQualityThreshold":0.2,"method":"TRIANGULATION",
                              "smoothingMethod":"GAUSS5x5","applyToOrtho":True,"regenPointCloud":False},
                     folder=prj_folder_name)

In [ ]:

dsm_item.url

In [ ]:

dsm_item.layers[0]

Out[ ]:

Summary¶

This article introduces the use of the ArcGIS Orthomapping Python API for performing Structure from Motion (SfM) photogrammetry. The article walks through the process of creating an image collection, adjusting images, and generating orthophotos, Digital Terrain Models (DTM), and Digital Surface Models (DSM) from drone photos. By following this guide, readers can harness the power of the ArcGIS Orthomapping Python API for accurate and efficient photogrammetry workflows.

More To Explore

Software Engineering

STEM at Your Fingertips: Public Libraries as Hubs of Innovation

In an age where knowledge is power and access to educational resources is key, public libraries stand out as beacons of learning and enlightenment. Beyond

Naheem March 31, 2024

ArcGIS

Structure from motion using ArcGIS orthomapping Python Api

In my previous articles on photogrammetry, I introduced structure from motion photogrammetry, explained block adjustment, and demonstrated how to generate photogrammetry products for urban monitoring

Naheem July 6, 2023