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.
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.
image_collection_item = gis.content.search("omImageCollection")[0]
image_collection_item.url'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
from arcgis.raster.orthomapping import compute_sensor_modelcompute_sensor_model(image_collection=image_collection_item, mode='Quick', location_accuracy='Low')'https://yourportal.domain.com/server/rest/services/Hosted/omImageCollection/ImageServer'
image_collection_item.layers[0]
To help us compare the results of orthomapping procedures better, let's save the first image we obtained from compute_sensor_model as img1:
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.
from arcgis.raster.orthomapping import generate_report
generate_report(image_collection=image_collection_item, report_format="PDF")'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.
from arcgis.raster.orthomapping import compute_control_pointscompute_control_points(image_collection=image_collection_item,
reference_image="https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer",
image_location_accuracy='Low')'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.
from arcgis.raster.orthomapping import match_control_pointsinput_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
}]
}]control_point_sets = match_control_points(image_collection=image_collection_item,
control_points=input_control_points,
similarity="High")control_point_sets[{'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()
from arcgis.raster.orthomapping import edit_control_pointscontrol_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}]}]edit_control_points(image_collection=image_collection_item, control_points=control_point_sets)'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.
from arcgis.raster.orthomapping import query_control_pointscontrol_points = query_control_points(image_collection=image_collection_item, query="OBJECTID>0")control_points[{'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.
compute_sensor_model(image_collection=image_collection_item, mode='Refine', location_accuracy='Medium')'https://yourportal.domain.com/server/rest/services/Hosted/omImageCollection/ImageServer'
generate_report(image_collection=imagecollectionItem, report_format="PDF")'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:
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:
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:
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:
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_pathout_path = display_images_diff("./img1.jpg", "./img2.jpg", "./diff_img1_img2.png")from IPython import display
display.Image(out_path)
Though not significant from the difference image shown above, because the accuracy and resolution of the mosaicking and georeferencing process of step 2.2.5 is much finer than that of step 2.2.1, img2 is more closely knitted to the ground truth, and shows fewer exterior skirts, similar to what we see from upsampling.
2.3 Compute seamlines
arcgis.raster.orthomapping.compute_seamlines() task is used to compute seamlines for overlaid images inside the image collection, usually after the image collection has been block adjusted. Seamlines are helpful for generating the seamless mosaicked display for overlapped images in image collection.
from arcgis.raster.orthomapping import compute_seamlinescompute_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})'https://yourportal.domain.com/server/rest/services/Hosted/omImageCollection/ImageServer'
image_collection_item.layers[0]
2.4 Compute color correction
After image collection’s sensor model is computed, arcgis.raster.orthomapping.color_correction() can be called to improve the visual appearance.
This is often needed for satellite data as the cross-strip images are often collected on different days with different atmospheric conditions.
from arcgis.raster.orthomapping import color_correctioncolor_correction(image_collection=image_collection_item, color_correction_method="DODGING", dodging_surface_type="SECOND_ORDER",
context={"skipRows": 10, "skipCols": 10, "reCalculateStats": "OVERWRITE"})'https://yourportal.domain.com/server/rest/services/Hosted/omImageCollection/ImageServer'
image_collection_item.layers[0]
The block adjustment tools allow for an iterative computation, so that you can check on the quality of the adjustment, modify options, add or delete GCPs, or recompute tie points before re-running the adjustment. If you are unsatisfied with the error in the Adjustment Report, try adding GCPs, or try modifying some of the Adjustment Options. You can also change some of your check points back into GCPs, and choose a few other GCPs to be your check points. Re-run the adjustment and see how this impacts the shift. Once you are satisfied with the accuracy of your adjusted imagery, it’s time to make ortho products [1]!
Step 3. Getting Results
3.1 Generate orthomosaic
After the block adjustment, ortho-rectified imagery layer can be created from the image collection with arcgis.raster.orthomapping.generate_orthomosaic()
This function comes with the option to choose whether or not using existing seamlines of the image collection for the orthomosaic imagery layer. Also, it has an option to choose whether or not to apply color correction settings to the output orthomosaic imagery layer. These can be set through the regen_seamlines and recompute_color_correction parameters.
from arcgis.raster.orthomapping import generate_orthomosaicortho_mosaic_name = "omOrthoMosaic"
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)ortho_mosaic_item.url'https://yourportal.domain.com/server/rest/services/Hosted/omOrthoMosaic/ImageServer'
ortho_mosaic_item.layers[0]
Same here, in order to compare the output of this step with previous image saved, let's save into a static image:
img5 = image_collection_item.layers[0].export_image(size=[1200,450], f='image', save_folder='.', save_file='img5.jpg')Then we calculate the MSE and SSIM in the back end, and label in the comparative display:
compare_images("./img1.jpg", "./img5.jpg", "Original vs. Enhanced", 1507.56, 0.66)
Also, using PIL.ImageChops to subtract one image from another, then visualize the difference image:
out_path = display_images_diff("./img1.jpg", "./img5.jpg", "./diff_img1_img5.png")
display.Image(out_path)
Now we have obtained the difference image from step 2.1.1 to step 3.1. Compare the difference image seen in step 2.2.6 to the difference we are seeing here, the latter shows more significant changes, since it not only includes the improvement of accuracy and resolution of the mosaicking and georeferencing process, but also has involved seamlines generation and color correction for visualization purposes.
3.2 Generate DEM
DEM products can be generated from image collections which have been block adjusted with arcgis.raster.orthomapping.generate_dem(). It creates dense point cloud using adjusted image collections, then interpolates the point cloud to generate different DEM surface types using the designated method.
Two types of DEMs that can be produced are Digital Terrain Model (DTM) and Digital Surface Model (DSM) [1]:
DTM — Digital elevation of the earth, not including the elevation of any objects on it, can be used to produce the orthoimage and orthomosaics. Because A DTM is a vector data set composed of regularly spaced points and natural features such as ridges and breaklines, and augments a DEM by including linear features of the bare-earth terrain, they are typically created through stereo photogrammetry. For instance, DTM can be non-continuous and not a surface model. With regularly-space and contour lines, you can interpolate a DTM into a DEM. Overall, a DTM represents distinctive terrain features much better because of its 3D breaklines and regularly spaced 3D mass points [5].
DSM — Digital elevation of the earth, including the elevation of objects on it such as trees and buildings, can be used for classifying features in orthoimages, such as discriminating asphalt pavement and asphalt roofs. It should not be used for image orthorectification unless the source imagery is nadir looking, with no building or feature lean, to produce true orthoimages. A DSM is useful in 3D modeling for telecommunications, urban planning and aviation. Because objects extrude from the Earth, this is particularly useful in examples of runway approach zone encroachmnt, vegetation management, and view obstruction [5].
3.2.1 Generate DTM
from arcgis.raster.orthomapping import generate_demdtm_name = "omDigitalTerrainModel"
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)dtm_item.url'https://yourportal.domain.com/server/rest/services/Hosted/omDigitalTerrainModel/ImageServer'
dtm_item.layers[0]
3.2.2 Generate DSM
dsm_name = "omDigitalSurfaceModel"
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)dsm_item.url'https://yourportal.domain.com/server/rest/services/Hosted/omDigitalSurfaceModel/ImageServer'
dsm_item.layers[0]
3.2.3 Use DEM to visualize scenes
From steps 3.2.1 and 3.2.2, we are getting two products, the DTM and the DSM. Let's visualize with the chimney structure at the study site (with varied elevation sources) to further understand the differences.
Fig 3. Applying different elevation sources to the same scene.
With the DEM generated here, and the orthomosaic image created from previous step, we can use ArcGIS Pro to create a short animation of the flyover of the current scene:
from IPython import display
x = display.Image(filename='../../static/img/guide_orthomapping_scene_dtm.gif')
y = display.Image(filename='../../static/img/guide_orthomapping_scene_dsm.gif')
display.display(x, y)
3.3 Calculate vegetation index (VARI) - persisting output as hosted image service
arcgis.raster.analytics.generate_raster() is used to generate rasters based on any raster function templates you build. The result will be stored inside server's raster store and serving out as hosted image services directly.
from arcgis.raster.analytics import generate_rastervari_raster_function_template = {
"rasterFunction" : "BandArithmetic",
"rasterFunctionArguments" : {
"Method" : 9, # 9 indicates the build-in Visible Atmospherically Resistant Index calculation
"BandIndexes" : "3 2 1", # Band number for Red, Green, Blue data used in VARI calculation.
"Raster" : {
"url" : ortho_mosaic_item.url
}
}
}vari_name = "omVARI"
vari_item = generate_raster(raster_function=vari_raster_function_template,
output_name=vari_name)vari_item.layers[0]
3.4 Calculate vegetation index (VARI) - through on-the-fly processing
In contrast with method above, arcgis.raster.functions module also contains many build-in raster functions that can be applied on imagery layer items on the fly for dynamic analysis.
from arcgis.raster.functions import band_arithmeticband_arithmetic(raster=ortho_mosaic_item.layers[0],
band_indexes="3 2 1",
method=9)
In Part 2, we have gone through processes to create a DTM or DSM, and an orthomosaic with aerial digital imagery, etc. using Python API methods. With these results, you can go on to do all kinds of analysis with elevation data, or use your orthomosaic map for applications in emergency response, resource management, real estate, and more [3].
Last but not least, we will talk about how to handle the processed imageries in Part 3, which includes how to add/delete images in the image collection, and reset/delete the image collection. Please stay tuned!
References
[1] Lenhardt, Ortho Mapping in ArcGIS Pro Part II: Getting Adjusted, https://www.esri.com/arcgis-blog/products/arcgis-pro/imagery/ortho-mapping-arcgis-block-adjustment/ [Online][Assessed on August 13, 2019]
[2] Lenhardt & Liedtke, Ortho Mapping in ArcGIS Pro Part I: Getting Organized with a Workspace, https://www.esri.com/arcgis-blog/products/arcgis-pro/imagery/ortho-mapping-workspace/ [Online][Assessed on August 13, 2019]
[3] Lenhardt, Ortho Mapping in ArcGIS Pro Part III: Getting Results,https://www.esri.com/arcgis-blog/products/arcgis-pro/imagery/ortho-mapping-products-arcgis-pro/ [Online][Assessed on August 13, 2019]
[4] "Block adjustment for mosaic datasets", http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/block-adjustment-for-mosaic-datasets.htm [Online][Assessed on Septermber 23, 2019]
[5] "DEM, DSM, DTM Differences", https://gisgeography.com/dem-dsm-dtm-differences/ [Online][Assessed on Septermber 23, 2019]