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.

Many more to come...

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 USGS Global Visual Viewer (glovis)

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"
ms_channels = [] #takes all input channels by default
pan_channel = []
pansharp_file = r"C:\landsat8\pansharp\psh_landsat8.pix"


pansharp2(ms_file, ms_channels, [], pan_file, pan_channel, pansharp_file,
          [], "", "", "", "", "")

Pansharpen satellite images - Batch Processing

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

#UNDER CONSTRUCTION

Co-register MS with PAN image and then pansharpen

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

#UNDER CONSTRUCTION

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 pci.link 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"C:\Level2\00025\Lvl02-00025-Pan.tif"
pan_link = r"C:\Lvl02-00025-Pan.pix"
link(pan_level2, pan_link, [])

ms_level2 = r"C:\Level2\00025\Lvl02-00025-Col.tif"
ms_link = r"C:\Lvl02-00025-Col.pix"
link(ms_level2, ms_link, [])

#get extents of pan image
header_info = []
ihr( pan_link, 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(pan_link, "Metre", [], "", "", "N", str(header_info[11]), str(header_info[14]), str(header_info[13]),
        str(header_info[12]), "", "", "", "", "", "", "", "", "", "", "", "", "", "", [], [])

#set meter projection for MS image (use same extents as PAN image and change resolution). Using same values
setpro2(ms_link, "Metre", [], "", "", "N", str(header_info[11]), str(header_info[14]), str(header_info[13]),
        str(header_info[12]), "", "", "", "", "", "", "", "", "", "", "", "", "", "", [], [])

#run Pansharp2 on MS and PAN imagery
level3_image = r"C:\output_level3\00025-PSH.pix"
pansharp2(ms_link, [1,2,3,4], [1,2,3,4], pan_link, [1], level3_image, [], "", "", "", "PIX", "")

#set Level3 (pansharpened image) back to pixel projection (no projection)
setpro2(level3_image, "PIXEL", [], "", "", "N", str(header_info[11]), str(header_info[12]), str(header_info[13]),
        str(header_info[14]), "", "", "", "", "", "", "", "", "", "", "", "", "", "", [], [])

#optionally convert it to a TIF image
level3_image_tif = r"C:\output_level3\00025-PSH.tif"
channels = [1,2,3,4]

fexport(level3_image, level3_image_tif, [], channels, [], [], [], [], "TIF", "")

Cloud Masking & Haze Removal

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

#UNDER CONSTRUCTION

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"C:\K3\RAW\K3_20131003174133_07365_20171339_L1R_Bundle\K3_20131003174133_07365_20171339_L1R_Aux.xml-MS"
channels = [1,2,3,4]
toa_image = r"C:\K3\toa\kompsat3_toa.pix"

dn2toa(raw_dn_image, channels, "", "", "", [], "Reflectance", toa_image, "", "")