// Parallel orthorectification of DigitalGlobe imagery using Polar Geospatial Center Tools on the ABoVE Science Cloud - Linux Environment

Overview

High-resolution DigitalGlobe imagery is available on the ABoVE Science Cloud (ASC) for approved ABoVE Science Team members. Request access here. Most of the DigitalGlobe imagery on the ASC is Basic (1B) imagery. DigitalGlobe describes this type of imagery in this way:
Basic (1B) Imagery is the least processed of the Core Imagery product series and is corrected for radiometric distortions, internal sensor geometry, optical distortions, and sensor distortions. Basic Imagery is neither geo-referenced nor mapped to a cartographic projection. Basic Imagery is provided with sensor models and is intended for sophisticated photogrammetric processing such as orthorectification. Basic Imagery is a scene-based product.”
~DigitalGlobe Core Imagery Product Guide

Depending on your use case, Basic (1B) imagery may require orthorectification in order to be used for location analysis. This guide describes the use of the Polar Geospatial Center’s “Imagery Utils” to orthorectify these data.

An ABoVE Science Cloud Webinar was held on this topic on 2/17/17. Find the webinar here and a link to the slides here.

Prerequisites:

  1. Copy Polar Geospatial Center’s orthorectification tools (https://github.com/PolarGeospatialCenter/imagery_utils) to your home directory:
    • When you log in to the ASC, you are already in your home directory (/att/home/[your user name]). Then use these commands: wget https://github.com/PolarGeospatialCenter/imagery_utils/archive/v1.3.2.zip
      unzip v1.3.2.zip
  2. Make an output folder where orthorectified images will be stored (e.g. backup, home directory) mkdir att/nobackup/[yourusername]/output
  3. Create a simple text file with the locations of the NGA imagery you would like to orthorectify (see Appendix 2 for how to locate NGA imagery on the ASC). Note: Each image location should be placed on a separate line. For example:
    • •[root_location]/NGA/WV02/1B/2010/204/WV02_10300100061FE500_X1BS_052402487010_01/WV02_20100723220557_10300100061FE500_10JUL23220557-M1BS-052402487010_01_P004.ntf
    • •[root_location]/NGA/WV02/1B/2010/204/WV02_10300100061FE500_X1BS_052402487010_01/WV02_20100723220557_10300100061FE500_10JUL23220557-M1BS-052402487010_01_P005.ntf
    • Note: Root locations of files on the ASC are not shared on the website. A version of this file is available in the ASC at [root_location]/NGA_footprints/Documentation/. If you are unsure of where this location, contact support@nccs.nasa.gov or Elizabeth.hoy@nasa.gov for assistance.
  4. Determine the EPSG code you would like to use to project your imagery. These EPSG codes can be found online here: epsg.io or http://spatialreference.org/ref/epsg/.
  5. Have a DEM available which you plan to use in the orthorectification process. On the ASC we have a DEM folder with multiple DEM datasets in it. The two most commonly used by PGC are:

Batch processing of NGA imagery:

  1. Log into Linux VM (e.g. ssh above101)
  2. Start by calling GDAL 2.0.1: source /opt/PGSC-2.0.1/init-gdal.sh
  3. Navigate to the home directory where PGC tools are located: cd ~/imagery_utils-1.3.2/
  4. Call orthorectification python script with desired parameters (see Appendix 1 for parameter options):
    1. Example: python pgc_ortho.py <image_paths_textfile.txt> <output_folder> -p 4326 -d /att/pubrepo/DEM/alaskaned_mosaic_wgs84.tif -t Byte -c rf --parallel-processes 2

Suggestions for script usage on ADAPT

Most users require one of four image processing parameter profiles:
  1. Basic visual analysis over ice (-t Byte -c rf): This profile is optimized to provide imagery for visual analysis over bright surfaces like ice and snow.It outputs a smaller 8-bit image instead of the full 11 bit radiometric range and balances the image colors.
  2. Basic visual analysis over non-ice (-t Byte -c mr): This profile is optimized to provide imagery for visual analysis over darker surfaces like rock and vegetation.It outputs a smaller 8-bit image instead of the full 11 bit radiometric range and balances the image colors, plus it brightens up the low end of the visual spectrum to provide good contrast in dark areas.
  3. Terrain correction only (-t UInt16 -c ns): For users who wish to use the scripts only for terrain correction, this set of parameters leaves the image DN values as is and only applies terrain correction according to the DEM you supply.
  4. Top-of-atmosphere reflectance or radiance (-t Float32 -c rf or -c rd): This set of parameters corrects the DN values to top-of-atmosphere reflectance (-c rf) or radiance (-c rd) and stores the result in a floating point raster.The resulting rasters take up a more space on your filesystem but are more useful for some quantitative analyses.
  5. Height - wgs84 vertical heights are needed for DEM inputs
The ADAPT system is a collection of VMs so the imagery_utils tool options intended for HPC systems using schedulers are not applicable here.  You can ignore the –pbs, --qsubscript, and -l options.  Instead, you should experiment with the parallel-processes option in order to process several images in parallel.

If you experience high I/O loads and low CPU usage even when the –parallel-processes option is near the number of cores on the VM, consider using a working directory that is a local scratch space so that your processing is not limited by network connectivity (--wd option).
 

Error Assessment

The RPC sensor model has well understood error margins, +-4m LE90 and CE90. However, the sensor model error assessment assumes a perfect elevation model. By far the greatest source error in images will come from the uncertainty of the DEM used to correct for terrain distortions, especially with high off-nadir images. Often, error is less than 30m, sometimes less that 5m, and, with a very poor DEM in a mountainous area in Antarctica, it can be as high as 500m.

Appendix 1: Input Parameters of pgc_ortho.py

-f {HFA, JP2OpenJPEG, ENVI, GTiff}
--format {HFA, JP2OpenJPEG, ENVI, GTiff}
Output format (default=GTiff)
HFA: IMAGINE .img format
JP2OpenJPEG: Lossless JPEG2000 format using the OpenJPEG2 driver
ENVI: Binary file with .envi extension
GTiff: GeoTiff format.  Compression can be specified using the --gtiff_compression argument.
--gtiff-compression {jpeg95, lzw} GTiff compression type (default=lzw)
-p EPSG
--epsg EPSG      
EPSG projection code for output files
 
-d DEM
--dem DEM
DEM to use for orthorectification (elevation values should be relative to the WGS84 ellipsoid).  If no values is supplied, an average elevation is taken from the sensor model for each image.
 
-t {Byte, UInt16, Float32}
--outtype {Byte, UInt16, Float32}
Output data type (default=Byte).
-r RESOLUTION
--resolution RESOLUTION
Output pixel resolution in units of the projection
-c {ns, rf, mr, rd}
--stretch {ns, rf, mr, rd}
Stretch type [ns: no stretch, rf: TOA reflectance (default), mr: modified reflectance, rd: absolute TOA radiance]
No stretch: scales the DN values to the output data type. Outtype Byte means values are scaled from 0-255.  UInt16 and Float32 keep the original 11 bit DN values from 0-2047.
Reflectance: calculates top of atmosphere reflectance and scales it to the output data type. Outtype Byte means values are scaled 0-200, UInt16 scales 0-2000, Float32 is unscaled percent reflectance (0-1). Modified reflectance: same as reflectance, but with a histogram stretch applied that brightens the lower end of the dynamic range.  mr stretch is used for non-ice covered areas.
Absolute radiance: calculates absolute top of atmosphere radiance.  UInt16 or Float 32 are the only appropriate data types to use for this stretch because the units are W/m2/micrometers.
  --resample {near, bilinear, cubic, cubicspline, lanczos} Resampling strategy - mimics gdalwarp options
--rgb Output multispectral images as 3 band RGB
--bgrn Output multispectral images as 4 band BGRN (reduce 8 band to 4)
-s
--save-temps
Save temp files. Used for debugging only.
--wd WD              Local working directory (default is the destination directory). Use this if I/O to your storage system is a processing bottleneck
  --skip-warp Skip warping step and only execute the radiometric manipulation step
--skip-dem-overlap-check Skip verification of image-DEM overlap
--no-pyramids Suppress calculation of output image pyramids and stats
--ortho-height ORTHO_HEIGHT Constant elevation to use for orthorectification (value should be in meters above the WGS84 ellipsoid).  This is only used to force orthorectification to a constant elevation other than that included in the image RPC sensor model.
--pbs Submit tasks to PBS. This option should not be used on ADAPT.
--parallel-processes PARALLEL_PROCESSES Number of parallel processes to spawn (default 1).  This should be increased on the ADAPT above nodes to leverage more than one core, but should not exceed the number of cores on the VM.
--qsubscript QSUBSCRIPT Qsub script to use in PBS submission (default is qsub_ortho.sh in the script folder).  Not applicable on ADAPT.
-l L PBS resources requested (mimics qsub syntax). Not applicable on ADAPT.
--dryrun              Print script actions without executing
 

Appendix 2: Locating DG Imagery on the ASC

  1. In [root_location]/NGA_footprints you can locate the most recent DigitalGlobe (DG) footprints files showing footprint locations of all the DG data accessible on the ASC.
    • In [root_location]/NGA_footprints, the first level of subfolders are organized by date of the footprints files such that MM_DD_YYYY defines the month (M)_Day(D)_Year(Y) of the footprints files
    • Within [root_location]/NGA_footprints /MM_DD_YYYY/, individual footprint files are organized by sensor, product level and year such that SSNN_PP_YYYY where SS is the sensor type, NN is the sensor number, PP is the product level (such as 1B) and YYYY is the year of acquisition.
    • Sensor Types:
      • GE – GeoEye
      • IK – Ikonos
      • QB – Quick Bird
      • WV – World View
  2. Users can download the footprint files and view them locally on your own machine or open them in the ASC using: qgis (Linux side of the ASC) or ArcGIS (Windows side of the ASC)
  3. The filepath location of each footprint is stored in the “S_FILEPATH” attribute of the table