Subset MODIS-Aqua (Clorofila-a)

Pegando o gancho com uma postagem [1] que já tinha guardada de um tempo e do pacote pyMODIS [2], resolvi fazer o meu subset de imagens MODIS em Python.

Para esse código, foram testadas apenas imagens L3 de Clorofila-a, já processadas pela agência espacial norte americana (NASA).

Para executar esse código, você precisará de uma imagem: Imagem MODIS de 9km.

Dê um bunzip2 no seu diretório, onde estão as imagens .bz2, para que o código abaixo identifique somente as imagens sem a extensão (.bz2).

Um detalhe, é que abaixo segue como fazer o processo em batelada, mas para o exemplo de uma imagem, já é possível visualizar a coisa.

Diferentemente do post [1], não encontrei dificuldades em utilizar os dados. Logo não fiz uma conversão em ascii e nem utilizei o HDFview.

[1] http://www.trondkristiansen.com/?page_id=479
[2] https://bitbucket.org/tylere/pymodis/overview

#/usr/bin/env python
# -*- coding: utf-8 -*-

'''
subset_MODIS

Programa para Fazer um subset de imagens MODIS
Dentro dos limites Lat e Long determinados.

Author =  Arnaldo Russo
e-mail =  arnaldorusso@gmail.com
'''
from pyhdf.SD import *
import numpy as np
import matplotlib.pyplot as plt
import scipy.io
import glob

#setting limits to be cut
LATLIMS = ([-40, -30])
LONLIMS = ([-40, -30])

indir = '/DATA/database/Images/L3/9km/'
outdir = '/DATA/database/Images/L3/resize/'

filelist = glob.glob(indir+'A*')
nfiles = len(filelist)
files = []
for path in filelist:
  files.append(path[len(indir):]) #remove path name

names = []
multi_img = []
for i in range(len(filelist)):
    A = SD(filelist[i])
    a = A.attributes()
    for k in xrange(0,len(a.keys())):
        nm = a.keys()[k]
        names.append(nm.replace(" ","_"))
    pin = dict(zip(names,a.values()[:]))
    lon = np.arange(pin['Westernmost_Longitude'],pin['Easternmost_Longitude'],pin['Longitude_Step'])
    lat= np.arange(pin['Northernmost_Latitude'],pin['Southernmost_Latitude'],-pin['Latitude_Step'])
    #Get the indices needed for the area of interest
    ilt = np.int(np.argmin(np.abs(lat-max(LATLIMS)))) #argmin catch the indices
    ilg = np.int(np.argmin(np.abs(lon-min(LONLIMS)))) #of minor element
    ltlm = np.int(np.fix(np.diff(LATLIMS)/pin['Latitude_Step']+0.5))
    lglm = np.int(np.fix(np.diff(LONLIMS)/pin['Longitude_Step']+0.5))
    #retrieve data SDS
    d = A.datasets()
    sds_name = d.keys()[0] #name of sds. Dictionary method.
    sds = A.select(sds_name)
    #load the subset of data needed for the map limits given
    P = sds[ilt:(ilt+ltlm),ilg:(ilg+lglm)]
    P[P==-32767]=np.nan #Rrs_670:bad_value_scaled = -32767s ;
    P=np.double(P)
    P=(pin['Slope']*P+pin['Intercept'])
    LT=lat[ilt+np.arange(0,ltlm-1)]
    LG=lon[ilg+np.arange(0,lglm-1)]
    Plg,Plt = np.meshgrid(LG,LT) #Further plots
    P = np.log(P) #chlorophyll mapping
    multi_img.append(P)


multi_img = np.asarray(multi_img)
plt.imshow(P)
plt.show()

### Save matrix
np.save("multi_img",multi_img)
# If you want to share with your Matlab friends Users =]
scipy.io.savemat(outdir+'multi_matrix.mat', mdict={'multi_img': multi_img})

Resultando nessa imagem de Clorofila-a