quinta-feira, 27 de abril de 2017

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

domingo, 23 de abril de 2017

Como criar uma crosstable utilizando Python Pandas?

Um exemplo de como utilizar o Pandas para formatar os dados criando uma tabela de referência cruzada:

import pandas as pd

dados = pd.read_csv("arquivo.csv")

dados.describe()

nfocos_ref  nfocos_atnpp  nf_aqua_terra  julianday  ano
count  366.000000  366.000000  366.000000  366.000000  366.0
mean  6.344262  52.486339  11.292350  183.500000  2016.0
std  11.279588  88.889148  18.364328  105.799338  0.0
min  0.000000  0.000000  0.000000  1.000000  2016.0
25%  0.000000  0.000000  0.000000  92.250000  2016.0
50%  0.000000  8.500000  2.000000  183.500000  2016.0
75%  8.000000  62.750000  15.000000  274.750000  2016.0
max  68.000000  546.000000  117.000000  366.000000  2016.0

pv=dados.pivot_table(index='ano', columns='julianday', values='nfocos_atnpp',aggfunc='sum')

ano/julianday  1  2  3  4  5  6  7  8  9  10
1998    NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN
1999    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
2009    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  2.0  2.0 
 
O links para a documentação oficial é:
http://pandas.pydata.org/pandas-docs/stable/generated/pandas.crosstab.html