Geospatial Image Processing¶

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)

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):
    os.mkdir(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:
    try:
        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

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"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")

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)

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)