Geospatial Image Processing

This section provides some basic recipes for image processing with Geomatica’s python functions. Included in this section are recipes for pansharpening and atmospheric correction.

Pansharpen satellite images

This recipe provides a simple way to pansharpen the multispectral and panchromatic bands of a satellite image. Specifically, this recipe is pansharpening Landsat-8 ortho images downloaded directly from the Earth Explorer website:

Note: the images must be coregistered

from pci.pansharp2 import pansharp2

ms_file = r"C:\landsat8\LC80150292013111LGN01\LC80150292013111LGN01_MTL.txt-MS"
pan_file = r"C:\landsat8\LC80150292013111LGN01\LC80150292013111LGN01_MTL.txt-PAN"
pansharp_file = r"C:\landsat8\pansharp\psh_landsat8.pix"

pansharp2(fili=ms_file, fili_pan=pan_file, filo=pansharp_file)

Pansharpen satellite images - Batch Processing

This recipe provides a simple way to batch pansharpen a directory of Landsat Images. Full Tutorial:

from pci.pansharp2 import pansharp2
from pci.exceptions import *
import os
import fnmatch

# Set the working directory
working_dir = r'I:\Data\Sensors\Landsat8'
output_dir = os.path.join(working_dir, 'Pan_Images')
if not os.path.isdir(output_dir):

input_files = []

for r, d, f in os.walk(working_dir):
    for inFile in fnmatch.filter(f, '*_MTL.txt'):
        input_files.append(os.path.join(r, inFile))

# Iterate through list of files and run
for image in input_files:
        pansharp2(fili='-'.join([image, 'MS']),
                  fili_pan='-'.join([image, 'PAN']),
                  filo=os.path.join(output_dir, '.'.join([os.path.basename(image).split('_MTL')[0], 'pix'])))
    except PCIException, e:
        print e
    except Exception, e:
        print e

Level 2 to Level 3 Airphoto Processing (Pansharpening)

This recipe provides a simple way to pansharpen level 2 aerial imagery in order to convert it to level 3 imagery.

from pci.ihr import ihr
from import link
from pci.setpro2 import setpro2
from pci.pansharp2 import pansharp2
from pci.fexport import fexport

#create pix live link file for MS and PAN images
pan_level2 = r"E:\Esri_cloud_Aerial\demo_DryRun\Level2\00025\Lvl02-00025-Pan.tif"
pan_link = r"E:\Esri_cloud_Aerial\demo_DryRun\working\Lvl02-00025-Pan.pix"

link(fili=pan_level2, filo=pan_link)

ms_level2 = r"E:\Esri_cloud_Aerial\demo_DryRun\Level2\00025\Lvl02-00025-Col.tif"
ms_link = r"E:\Esri_cloud_Aerial\demo_DryRun\working\Lvl02-00025-Col.pix"

link(fili=ms_level2, filo=ms_link)

#get extents of pan image
header_info = []

ihr(file=pan_link, imstat=header_info)

'''ihr returns the extents of the imagery to the output list 'header_info':
header_info[11] = Upper Left X
header_info[12] = Upper Left Y
header_info[13] = Lower Right X
header_info[14] = Lower Right Y'''

#set meter projection for pan (have to flip the ULY and LRY values for meter projection
setpro2(file=pan_link, project="METRE", llbound="N", ulx=str(header_info[11]), uly=str(header_info[14]),
        lrx=str(header_info[13]), lry=str(header_info[12]))

#set meter projection for MS image (use same extents as PAN image and change resolution). Using same values
setpro2(file=ms_link, project="METRE", llbound="N", ulx=str(header_info[11]), uly=str(header_info[14]),
        lrx=str(header_info[13]), lry=str(header_info[12]))

#run Pansharp2 on MS and PAN imagery
level3_image = r"E:\Esri_cloud_Aerial\demo_DryRun\output_level3\00025-PSH.pix"

pansharp2(fili=ms_link, dbic=[1,2,3,4], dbic_ref=[1,2,3,4], fili_pan=pan_link, dbic_pan=[1], filo=level3_image)

#set Level3 (pansharpened image) back to pixel projection (no projection)
setpro2(file=level3_image, project="PIXEL", llbound="N", ulx=str(header_info[11]), uly=str(header_info[12]),
        lrx=str(header_info[13]), lry=str(header_info[14]))

#optionally convert it to a TIF image
level3_image_tif = r"E:\Esri_cloud_Aerial\demo_DryRun\output_level3\00025-PSH.tif"

fexport(fili=level3_image, filo=level3_image_tif, dbic=[1,2,3,4], ftype="TIF")

Cloud Masking & Haze Removal with optional Atmospheric Correction

This recipe outlines an ATCOR workflow. Depending on your requirements, you can simply run masking and hazerem to produce a haze free image.

import os
from pci.masking import masking
from pci.hazerem import hazerem
from pci.atcor import atcor
from pci.tersetup import tersetup
from pci.illumcast import illumcast

working_dir = r'C:\ATCOR\Landsat8'
ms_image = os.path.join(working_dir, 'Raw', 'LC80180302014247LGN00_MTL.txt-MS')
dem = os.path.join(working_dir, 'DEM', 'dem_30m.pix')
output_dir = r'C:\ATCOR\Mask_Haze_Atcor'
masks = os.path.join(output_dir, "LC80090272016206LGN00-Masks.pix")
hazefree_ms = os.path.join(output_dir, 'LC80090272016206LGN00_MS_HazeFree.pix')
ter_file = os.path.join(output_dir, 'dem_30m_Tersetup.pix')
illum_file = os.path.join(output_dir, 'dem_30m_Illumcast.pix')
atcor_ms = os.path.join(output_dir, 'LC80090272016206LGN00_MS_Atcor.pix')

# Create cloud, water and haze masks
masking(fili=ms_image, hazecov=[50], clthresh=[18, 22, 1], filo=masks)

# Remove the haze from both the multispectral and panchromatic bands
hazerem(fili=ms_image,  maskfili=masks, hazecov=[50], filo=hazefree_ms)

# Calculate the slope, aspect, and skyview rasters in preparation for atmospheric correction
tersetup(filedem=dem, terfile=ter_file)
illumcast(fili=ter_file, filref=ms_image, filo=illum_file)

# Calculate atmospheric correction for multispectral image.
atcor(fili=hazefree_ms, maskfili=masks, terfile=ter_file, illufile=illum_file, atmdef="rural", atmcond="Fall",
      outunits="Percent_Reflectance", filo=atcor_ms)

Calculate Top-of-atmosphere (TOA) Reflectance image

This recipe performs a top-of-atmosphere (TOA) calculation for satellite imagery. The output image is a TOA reflectance image, where pixels will be provided a value between 0 and 100%.

from pci.dn2toa import dn2toa

#Run Top of Atmosphere Reflectance computation
raw_dn_image = r"D:\K3_20131003174133_07365_20171339_L1R_Bundle\K3_20131003174133_07365_20171339_L1R_Aux.xml-MS"
channels = [1,2,3,4]
toa_image = r"D:\K3\toa\kompsat3_toa.pix"

dn2toa(fili=raw_dn_image, dbic=channels, outunits="Percent_Reflectance", filo=toa_image)