Skip To Content ArcGIS for Developers Sign In Dashboard

ArcGIS API for Python

Guide to Orthomapping (Part 2)

Continued from Part 1

In Part 1 we have learnt the concept of Orthomapping, the pre-requisite configuration environment for orthomapping tools to run on, and how to organize imageries and create imagery collection layer item.

Next, Part 2 of this guide will cover how to apply block adjustments, manipulate control points, and compute seamlines and color correction in Step 2, and procedures in generating orthomosaic, DEM, and Vegetation Indices through orthomapping tools in Step 3.

Let's get started! First, get connected to your GIS instance.

In [1]:
import arcgis
from arcgis.gis import GIS

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)

The Image Collection Layer item has already been created in Part 1, so let us reuse the reference to the layer item here.

In [3]:
image_collection_item = gis.content.search("omImageCollection")[0]
image_collection_item.url
Out[3]:
'https://yourportal.domain.com/server/rest/services/Hosted/omImageCollection/ImageServer'

Step 2. Getting Adjusted

In Part 1 of this guide, we explained what an ortho image collection is and how to create one for digital aerial imagery. At this point, the imagery has been organized and managed so that we can access all the necessary metadata, information, tools and functionality to work with our imagery, but we haven’t yet performed a bundle block adjustment [1].

In order to understand bundle block adjustment better, let's first rewind on some of these terms:

  • Tie Points "A tie point is a point that may not have known ground coordinates, but is visually recognizable in the overlap area between two or more images. The corresponding positions of the tie point in each overlapped image are identified and measured. The figure below shows tie points depicted as black diamonds." [4] Tie points are used to minimize the misalignment between two or more images.
  • Ground Control Points "Ground control points, shown as blue stars in the figure below, are identifiable features located on earth surfaces that have known ground coordinates x,y, and z." Generally, the GCPs can be prepared in three ways: (1) Survey points that come in a form of x,y,z. You can input the coordinates into the control point table manually. (2) Ortho rectified image chips. If you have image chips for some features or locations within your block, create a mosaic dataset from these image chips, and compute control points from your mosaic dataset using the Compute Control Points tool, then add it to the control point table. And (3) Ortho rectified image service. You can also compute control points based on a single ortho rectified image or image service, and add them to the control point table check points. "Check points are known ground control points, but instead of them being used in computing adjustment, they are used in checking the accuracy of adjustment computation." [4] GCPs are used to georeference the images to the ground.
  • Block "Mosaicking satellite images or aerial photographs for an area or a project", is called a block. For example, the nine overlapping images or photographs are for a single area, which can be defined as a block.
  • Block Adjustment A Process that consists of the three key parts - Tie points, GCPs, and Triangulation (which computes the transformation by minimizing and distributing the errors among images and control points).

Fig 1. Inside each mosaic item, the blue stars represent ground control points (GCP), and the black diamonds represent the tie points (Source: [4])

2.1 Compute and apply initial block adjustment

The block adjustment process is defined as "using the ground control and tie point information, a bundle adjustment computation calculates the exterior orientation for each image, such that they are consistent with neighboring images. The orientation for the whole block of images is then adjusted to fit the ground. This block adjustment process produces the best statistical fit between images, for the whole contiguous block, minimizing errors with the ground control. The adjusted transformation for each image item comprising the block is stored in the workspace for the orthomosaic" [2].

arcgis.raster.orthomapping.compute_sensor_model() computes the bundle block adjustment for the image collection and applies the frame transform to the images, which is based on the raster type information we used when creating the image collection. It also generates the control points feature class, solution table, solution points feature class, and flight path feature class. (These feature classes and tables will be saved in raster store and will not be created in the ArcGIS Enterprise as separated items.)

For preliminary quality purposes, we can perform block adjustment at reduced image resolution using "Quick" mode. This API also supports "Full" and "Refine" modes. Running under these two modes will make it ready for another block adjustment run with more ground control points and matching tie points later on. Once the image collection is in "full" block adjustment states, it will use the existing control point table in the further block adjustments. "mode" parameter will be optional for block adjustment at this time.

It also offers an adjustment report which can be accessed using arcgis.raster.orthomapping.generate_report().

2.1.1 Compute sensor model

In [19]:
from arcgis.raster.orthomapping import compute_sensor_model
In [1]:
compute_sensor_model(image_collection=image_collection_item, mode='Quick', location_accuracy='Low')
Out[1]:
'https://yourportal.domain.com/server/rest/services/Hosted/omImageCollection/ImageServer'
In [1]:
image_collection_item.layers[0]
Out[1]:

To help us compare the results of orthomapping procedures better, let's save the first image we obtained from compute_sensor_model as img1:

In [ ]:
img1 = image_collection_item.layers[0].export_image(size=[1200,450], f='image', save_folder='.', save_file='img1.jpg')

2.1.2 Generate orthomapping report (optional)

The report would contain information about the quality of the adjusted images, the distribution of the control points etc.

In [2]:
from arcgis.raster.orthomapping import generate_report
generate_report(image_collection=image_collection_item, report_format="PDF")
Out[2]:
'https://yourportal.domain.com/server/rest/directories/arcgisjobs/system/orthomappingtools_gpserver/jfe3fdb2cf00141ab812637b524b6d813/scratch/orthoreport_20190819105129.pdf'

2.2 Add ground control points (optional)

In photogrammetry industry, the orthomapping workflow for drone imagery may need to run block adjustment twice. The initial block adjustment is computed to generate the tie points among overlapping images for preliminary results. The second block adjustment is computed with ground control points for refined results.

Note that, it is optional to compute/match/edit/query control points, or recompute sensor model for refined results, as stated in sections 2.2.1 to 2.2.5.

2.2.1 Compute control points

arcgis.raster.orthomapping.compute_control_points() tool will compute the control point sets between the reference image and the mosaicked image items within the image collection.

Each control point set will be identified with the same "PointID" in the whole control point feature class.

Each control point set will consist of one ground control point and matching tie points on all available mosaicked image items of the image collection.

These new control point sets will be added to the control point feature class of the image collection automatically for computing block adjustment later on.

A georeferenced imagery data can be used as reference image layer in order to compute the control points. For example, a previously georeferenced data covering the same area, or ArcGIS Online World Imagery. The tool uses advanced image matching algorithm to find ground control points on reference image and its tie points on image collection item. It requires the reference image and the image collection having similar spatial resolution and texture for better image matching results.

By default, the mode option is set to Quick which means to compute tie points and adjustment at 8x of the source imagery resolution. The image_location_accuracy option allows users to specify the GPS location accuracy level of the source image. It determines how far the underline tool will search for neighboring matching images, then calculate tie points and compute adjustments. When set to Low, it means that GPS accuracy is of 20 to 50 meters, and the tool uses a maximum of 4 by 12 images.

In [24]:
from arcgis.raster.orthomapping import compute_control_points
In [4]:
compute_control_points(image_collection=image_collection_item, 
                       reference_image="https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer",
                       image_location_accuracy='Low')
Out[4]:
'https://yourportal.domain.com/server/rest/services/Hosted/omImageCollection/ImageServer'

2.2.2 Match control points

Alternatively, we can run match control points first, then run edit control points to update the control point feature class.

With given ground control points, arcgis.raster.orthomapping.match_control_points() will help us find all matching tie points on image items of the image collection.

For each given ground control point, it is required to set at least another matching tie point on image item as input. Then, this tool will help find all remaining tie points on all other image items of the image collection.

The result of this tool will be an array of control point sets. It can be used as input of arcgis.raster.orthomapping.edit_control_points() tool of orthomapping workflows.

In [26]:
from arcgis.raster.orthomapping import match_control_points
In [27]:
input_control_points = [{"pointId": 1,"x": -117.0926538,"y": 34.00704253,"z": 634.2175,
                       "spatialReference": {"wkid":4326},
                       "xyAccuracy": "0.008602325","zAccuracy": "0.015",
                       "imagePoints": [{
                           "imageID": 13,
                           "x": 4932.0895715254455,
                           "y": -1833.8401744114286
                       }] 
                     }]
In [28]:
control_point_sets = match_control_points(image_collection=image_collection_item, 
                                          control_points=input_control_points, 
                                          similarity="High")
In [29]:
control_point_sets
Out[29]:
[{'pointID': 1,
  'x': -117.09265379999133,
  'y': 34.007042530296786,
  'z': 634.2174999999988,
  'spatialReference': {'wkid': 4326, 'latestWkid': 4326},
  'xyAccuracy': 0.008602325,
  'ZAccuracy': 0.015,
  'imagePointSpatialReference': {'ics': '', 'topup': True},
  'imagePoints': [{'imageID': 13,
    'x': 4932.087025858724,
    'y': -1833.8399285802177,
    'u': 948.756293148009,
    'v': -545.7443838746151}]}]

2.2.3 Edit control points

Append additional control point sets to the image collection’s control point feature class with arcgis.raster.orthomapping.edit_control_points()

In [30]:
from arcgis.raster.orthomapping import edit_control_points
In [31]:
control_point_sets=[{"type":2,"x":-117.0926538,"y":34.00704253,"z":634.2175,
                     "pointId":1,"xyAccuracy":"0.008602325","zAccuracy":"0.015",
                     "imagePoints":[{"imageID":12,"x":5101.184543436052,"y":-2841.9339596208747,"u":545.0939944525426,"v":-468.5878186407665},
                                    {"imageID":2,"x":2020.0648901122495,"y":-2412.6233293180594,"u":3128.136467198136,"v":-2662.994335058943},
                                    {"imageID":13,"x":4932.089326515454,"y":-1833.840551626632,"u":555.357416431812,"v":-1555.0023564815583}]}]
In [5]:
edit_control_points(image_collection=image_collection_item, control_points=control_point_sets)
Out[5]:
'https://yourportal.domain.com/server/rest/services/Hosted/omImageCollection/ImageServer'

2.2.4 Query control points

An image collection has corresponding control points feature class after computing its sensor model. arcgis.raster.orthomapping.query_control_points() can help query among control point sets which have ground control points inside.

In [33]:
from arcgis.raster.orthomapping import query_control_points
In [ ]:
control_points = query_control_points(image_collection=image_collection_item, query="OBJECTID>0")
In [6]:
control_points
Out[6]:
[{'pointID': 9461,
  'x': -117.09269999992102,
  'y': 34.006999999284744,
  'z': 634.2174999999988,
  'spatialReference': {'wkid': 32611, 'latestWkid': 32611},
  'xyAccuracy': 0.008602325,
  'ZAccuracy': 0.015,
  'imagePointSpatialReference': {'wkid': 32611, 'latestWkid': 32611},
  'imagePoints': [{'imageID': 12,
    'x': 491448.85840000026,
    'y': 3762905.8003000002,
    'u': 545.093994453,
    'v': -468.587818641},
   {'imageID': 2,
    'x': 491425.1074999999,
    'y': 3762931.9886000007,
    'u': 3128.1364672,
    'v': -2662.99433506},
   {'imageID': 13,
    'x': 491447.1553999996,
    'y': 3762935.4902999997,
    'u': 555.357416432,
    'v': -1555.00235648}]}]

2.2.5 Recompute sensor model

After adding the ground control point sets to the image collection's control point table with steps above, we can rerun compute sensor model to refine block adjustment result.

Here, mode='Refine' represents to adjust the image at 1x of the source imagery resolution, while location_accuracy='Medium' means that GPS accuracy is of 10 to 20 meters, and the tool uses a maximum of 4 by 6 images.

In [7]:
compute_sensor_model(image_collection=image_collection_item, mode='Refine', location_accuracy='Medium')
Out[7]:
'https://yourportal.domain.com/server/rest/services/Hosted/omImageCollection/ImageServer'
In [6]:
generate_report(image_collection=imagecollectionItem, report_format="PDF")
Out[6]:
'https://yourportal.domain.com/server/rest/directories/arcgisjobs/system/orthomappingtools_gpserver/jf8aafbaf829a4ba88b860f4a1e818912/scratch/orthoreport_20190924124938.pdf'

Comparing the two ortho mapping reports (the previous report can be seen from step 2.2.1), we can see the differences in tie points, solution points, and the mean projection error.

Fig 2. Left column represents results in step 2.2.1, while the right column shows results obtained from step 2.2.5

Now that after steps we have done with adding, matching, editing, and querying control points, and recomputing the sensor model, let's save the output imagery layer as a static image, so we can compare with the previous image saved:

In [ ]:
img2 = image_collection_item.layers[0].export_image(size=[1200,450], f='image', save_folder='.', save_file='img2.jpg')

2.2.6 Compare two images

The compare_images method below is to show the original and the enhanced images comparatively, and also labeling with the Mean Squared Error and the Structural Similarity index between the two images:

In [1]:
import matplotlib.pyplot as plt
import cv2

""" compare two images, and label the mean squared error and 
    structural similarity index for the images
"""
def compare_images(imageA_path, imageB_path, title, m, s):
    # load the images -- the original, and the enhanced
    original = cv2.imread(imageA_path)
    enhanced = cv2.imread(imageB_path)
    # setup the figure
    fig = plt.figure(title, figsize=(20, 5))
    plt.suptitle("MSE: %.2f, SSIM: %.2f" % (m, s), fontsize=18)
 
    # show first image
    ax = fig.add_subplot(1, 2, 1)
    plt.imshow(original, cmap = plt.cm.gray)
    plt.axis("off")
 
    # show the second image
    ax = fig.add_subplot(1, 2, 2)
    plt.imshow(enhanced, cmap = plt.cm.gray)
    plt.axis("off")
 
    # show the images
    plt.show()

The Mean Squared Error (MSE) and the Structural Similarity index (SSIM) between the two images are computed at the back end using skimage.measure module. As shown below, the difference from step 2.2.1 to 2.2.5 is visualized:

In [2]:
compare_images("./img1.jpg", "./img2.jpg", "Original vs. Enhanced", 490.40, 0.68)

Or using the PIL library to perform ImageChop from one image to another, and then display the difference image:

In [3]:
from PIL import Image, ImageChops

def display_images_diff(imageA_path, imageB_path, diff_out_path):
    im1 = Image.open(imageA_path)
    im2 = Image.open(imageB_path)
    diff = ImageChops.difference(im1, im2)
    diff.save(diff_out_path)
    return diff_out_path
In [4]:
out_path = display_images_diff("./img1.jpg", "./img2.jpg", "./diff_img1_img2.png")
In [5]:
from IPython import display
display.Image(out_path)
Out[5]: