Processamento de Imagens utilizando python e GDAL

Neste Exemplo são apresentados métodos de leitura de arquivos matriciais fazendo uma varredura em  uma pasta para então processar todos os arquivo. Um ponto importante é que os dados foram lidos como array multidimensional, portanto sem informação espacial, porém o arquivo de saída foi  transformado em um Geotiff.

import glob, os, sys
import numpy as np
import osgeo.gdal as gdal
from osgeo import osr


path_in_aqm = r"c:\Projetos\AQM_Renata\Area_Queimada_modis_americadosul\aqm"  #2005_01_AQM.bin
path_in_foco = r"c:\Projetos\AQM_Renata\Area_Queimada_modis_americadosul\focos" #2005_01_grd_focos.bin
path_out = r"c:\Projetos\AQM_Renata\Area_Queimada_modis_americadosul\aqm_focos_tif"

anos = range(2005,2016)
meses = ["%02d"%a for a in range(1,13)]
#ano = 2014
#mes = "06"

for ano in anos:
    for mes in meses:

        arq_aqm = os.path.join(path_in_aqm,"%s_%s_AQM.bin"%(ano,mes))
        arq_foco = os.path.join(path_in_foco,"%s_%s_grd_focos.bin"%(ano,mes))
        arq_saida = os.path.join(path_out,"%s_%s_aqm_focos.tif"%(ano,mes))

        img_aqm = np.fromfile(arq_aqm,dtype=np.int16,count=-1)
        img_foco = np.fromfile(arq_foco,dtype=np.int16,count=-1)

        img_aqm = np.resize(img_aqm,(6300,5000))
        img_foco = np.resize(img_foco,(6300,5000))

        idx_aqm = np.where(img_aqm==1,20,0)
        idx_foco = np.where(img_foco>0,40,0)

        nova_area = (idx_aqm+idx_foco).astype("uint8")

        #nova_area.tofile(path_in_aqm + "/comparacao_%s_%s_uint8.bin"%(ano,mes))

        # definicioes para gravar um aquivo no formato envi usando gdal
        spatialRef = osr.SpatialReference()
        spatialRef.SetWellKnownGeogCS( "EPSG:4326" )
        geotransform = ( -83.005, 0.01 , 0 , 13.005, 0 , -0.01)

        driver = gdal.GetDriverByName("GTIFF")

        ds = driver.Create(arq_saida, 5000,6300,1,gdal.GDT_Byte , ['COMPRESS=LZW'])
        ds.GetRasterBand(1).WriteArray(nova_area)

        ds.SetProjection(spatialRef.ExportToWkt())
        ds.SetGeoTransform(geotransform)

        ds=None

print "Final do Processamento"
print "-"*10