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
Nenhum comentário:
Postar um comentário