tag:blogger.com,1999:blog-49982108226567971462024-03-05T23:20:14.754-03:00ciclotux - pedalando livreLinux, Debian, Python, Remote Sensing, Sensoriamento Remoto, numpy, scipy, matplotlib, Basemap, Ocean Color, MODIS, SeaWiFS, Foucault, Sartre, BicicletaArnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.comBlogger14125tag:blogger.com,1999:blog-4998210822656797146.post-56525962436079031452013-11-28T12:45:00.003-02:002013-11-28T12:45:31.789-02:00Ciclotux.blogspot.com está desativado!<br />
<br />
Deixo <a href="http://www.ciclotux.org/">aqui</a> o redirecionamento para a continuidade desse bloBs, Python Powered!<br />
<br />
<h4>
<span style="color: #990000; font-family: Verdana, sans-serif;"><b style="background-color: white;"><a href="http://ciclotux.org/">CicloTux.org</a></b></span></h4>
<br />
Se você mantém algum feed desse blog, por favor modifique o seu apontador.<br />
<br />
Até mais.Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-85243691909629517642013-11-18T02:34:00.000-02:002013-11-18T08:51:11.415-02:00O amor tambem fala PythonEncontrei uma maneira elegante para dizer " Baixinha, eu te amo!"<br />
Claro que é um tanto Pythonico, mas como o amor pode ser expresso em qualquer linguagem, por que não em Python?<br />
<pre class="prettyprint">In [3]: Baixinha(eu[te <3 ])
Out [3]: 'E eh inf!'
</pre>
Eu digo, "Baixinha eu te amo"
<p>Como eu sou muito amigo do Python... ele já sabe disso e complementa.
<b>"E eh infinito"!</b>
<p>Como o mínimo de amor que ela aceita é um <b><3</b> , escrevi uma função Python para isso.
<pre class="prettyprint">import numpy as np
eu = 'eu'
te = 'te'
def Baixinha(val):
if val < '<3':
raise ReferenceError
else:
return 'E eh %f!' % np.inf
</pre>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh8aoStbRKmbgS3Lk7QTXn7xZLzdqPKze6QpOPC4R1wm8y64doDEdstam9h2pTDiL1jNQc7qryzwC5EhlA01h5x9Tdp7DNycQzaoyo2zmnjDU8iWfhDke8sZRndEjEm5kUTPYQkeFD3QDZT/s1600/github_pythocat.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh8aoStbRKmbgS3Lk7QTXn7xZLzdqPKze6QpOPC4R1wm8y64doDEdstam9h2pTDiL1jNQc7qryzwC5EhlA01h5x9Tdp7DNycQzaoyo2zmnjDU8iWfhDke8sZRndEjEm5kUTPYQkeFD3QDZT/s320/github_pythocat.png" /></a>Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com2tag:blogger.com,1999:blog-4998210822656797146.post-49317507637747615772013-01-28T14:54:00.000-02:002013-01-28T14:54:05.351-02:00Data Frame - R and PythonPara acessar alguns dados com <b>R</b>, e poder realizar operações, utilizamos a função <b>data.frame</b> podendo agrupando fatores entre outras operações. Em Python a coisa pode ser bastante parecida e com diversas funcionabilidades do pacote <b><a href="http://pandas.pydata.org/" target="_blank">Pandas</a> [1]</b> para Python.<br />
<br />
Digamos uma tabela da seguinte maneira:<br />
<pre class="prettyprint">ID sp1 sp2 sp3
C1 1 3 2
C1 3 2 0
C1 2 1 1
C2 2 4 0
C2 1 3 1
C3 0 3 4
C3 2 1 2
</pre>
<br />
<pre class="prettyprint">#! /usr/bin/R
> dados = read.table('dados.txt', header=T)
> dados
ID sp1 sp2 sp3
1 C1 1 3 2
2 C1 3 2 0
3 C1 2 1 1
4 C2 2 4 0
5 C2 1 3 1
6 C3 0 3 4
7 C3 2 1 2
> attach(dados)
# transformando a coluna ID em fator.
> dados$ID <- data-blogger-escaped-colunas="" data-blogger-escaped-dados.="" data-blogger-escaped-dados="" data-blogger-escaped-das="" data-blogger-escaped-de="" data-blogger-escaped-factor="" data-blogger-escaped-uma="" data-blogger-escaped-ver=""> dados$sp1
[1] 1 3 2 2 1 0 2
# para fazer uma média, por exemplo, utilizando o fator ID.
> xsp1 <- data-blogger-escaped-mean="" data-blogger-escaped-sp1="" data-blogger-escaped-tapply="" id=""> xsp1
C1 C2 C3
2.0 1.5 1.0
<!-----><!-----></pre>
<br />
Para fazer isso em Python, você poderia usar o <b>numpy</b>, simplesmente, ou usar as elegantes ferramentas do "pandas". Para quem vem do <b>R</b>, o <b><a href="http://pandas.pydata.org/" target="_blank">Pandas</a></b> é uma forma fácil de se entender com dados tabelados e cheios de fatores.<br />
<pre class="prettyprint">>>> import pandas as pd
>>> dados = pd.read_table('dados.txt',sep=" ")
>>> dados
ID sp1 sp2 sp3
0 C1 1 3 2
1 C1 3 2 0
2 C1 2 1 1
3 C2 2 4 0
4 C2 1 3 1
5 C3 0 3 4
6 C3 2 1 2
# Acessar somente uma das colunas.
>>> dados['sp1']
0 1
1 3
2 2
3 2
4 1
5 0
6 2
Name: sp1
# Para fazer a média, com a função groupby
>>> dados.groupby('ID').sp1.mean()
ID
C1 2.0
C2 1.5
C3 1.0
Name: sp1
# Ainda é possível fazer várias coisas com o objeto dados (que agora é um DataFrame).
# Por exemplo, pegar a média e a variância, ordenados por ID.
>>> dados.groupby('ID').sp1.agg(['mean', 'var'])
mean var
ID
C1 2.0 1.0
C2 1.5 0.5
C3 1.0 2.0
</pre>
<br />
Até mais.<br />
[1] <a href="http://pandas.pydata.org/">http://pandas.pydata.org/</a>Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-53852313349312933502012-12-10T22:45:00.001-02:002012-12-10T22:45:45.870-02:00Latitude e Longitude - Ice cover - NSIDCAproveitando o tempo de recesso sem postagens, vou falar um pouquinho de como fiz para construir as matrizes de lat e long, dos dados de gelo do NSIDC, com as dicas da Claudinha Parise.<br />
<br />
Os dados estão em formato binário, como integer e armazenados de 4 em 4 bytes (em Little endian)[1]. Não sei se o matlab chegaria a dar problemas com isso, mas em Python é bom especificar. ("<b><span style="color: blue;"><</span></b>" Little Endian; "<b><span style="color: blue;">></span></b>" Big Endian); o "i" é de integer, "u" seria de string(unicode)[2]. Essa foi a dica do Irber; pois eu tava usando dtype = 'uint32', pensando somente na questão dos 4 bytes e esquecendo do "little endian".<br />
<br />
Dessa forma, baixa-se o dado de latitude para o hemisferio sul [<a href="ftp://sidads.colorado.edu/pub/DATASETS/brightness-temperatures/polar-stereo/tools/geo-coord/grid/pss25lats_v3.dat" target="_blank">aqui</a>] e divide-se por 100000 para ter os graus decimais. Mesmo procedimento para a Longitude.<br />
<br />
<pre class="prettyprint">import numpy as np
f = open('pss25lats_v2.dat','rb')
lat = np.fromfile(f, dtype='<i4')/100000.
lat = lat.reshape(332,316)
f.close()
</pre>
<br />
<br />
[1] <a href="http://nsidc.org/data/polar_stereo/tools_geo_pixel.html">http://nsidc.org/data/polar_stereo/tools_geo_pixel.html</a><br />
<br />
[2] <a href="http://docs.scipy.org/doc/numpy/user/basics.byteswapping.html">http://docs.scipy.org/doc/numpy/user/basics.byteswapping.html</a><br />
<br />Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-49803113592151170842012-11-01T09:33:00.001-02:002012-11-01T09:37:18.186-02:00Loops 'tradicionais'Meu camarada <a href="http://blog.luizirber.org/" target="_blank">Irber</a>, relembrando como fazer um loop 'tradicionalista' onde se pegam os índices dos ítens de uma lista ou de um numpy array, usando o enumerate, para pegar o valor e seu "índice".<br />
<br />
Buscando, dessa forma, um loop tradicionalista, porém mais pythonico (eu acho).<br />
<pre class="prettyprint">valores = ([1,3,5,7,9])
for i, value in enumerate(valores):
print i, value
Out[2]:
0 1
1 3
2 5
3 7
4 9
</pre>Sendo o jeito um tanto mais "conservador" e talvez menos pythonico:<br />
<pre class="prettyprint">for i in xrange(len(valores)):
print i, valores[i]
Out[4]:
0 1
1 3
2 5
3 7
4 9
</pre>Até mais.Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-75835373185767410182012-09-28T09:26:00.002-03:002012-09-28T09:26:59.089-03:00Livros II - Oceanic Sciences and HydrologyUm livro interessante, sobre Python aplicado às ciências atmosféricas e oceânicas.<br />
<a href="http://www.johnny-lin.com/pyintro/"></a><a href="http://www.johnny-lin.com/pyintro/" target="_blank">A Hands-On Introduction to Using Python in the Atmospheric and Oceanic Sciences</a><br />
Outro livro bacana é sobre a utilizacão do Python em hidrologia.<br />
<a href="http://greenteapress.com/pythonhydro/pythonhydro.pdf" target="_blank">Python in Hydrology</a><br />
Boas leituras.Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-496596626957826402012-07-31T11:20:00.000-03:002012-07-31T11:20:00.003-03:00Livro: Estudos OceanográficosEstá disponível para download o livro:<br />
<a href="http://www.oceano.furg.br/sistema/upload_php/estudos_oceanograficos.pdf"><b>Estudos Oceanográficos - do instrumental ao prático</b></a><br />
Organizado pelo Professor Danilo Calazans - FURG.<br />
<br />
Outros livros estão disponíveis no site do <a href="http://www.cdmb.furg.br/"><b>Ciências do Mar</b></a>.<br />
<br />
Boas leituras!Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-244103288153856462012-06-08T10:37:00.000-03:002012-06-08T10:41:24.640-03:00Papers - Oceanography and Python<br />
Dando parabéns aos Oceanógrafos pelo "dia do Oceanógrafo", seguem algumas leituras do dia de hoje.<br />
<br />
[1] - <a href="http://www.sciencedirect.com/science/article/pii/S0098300401000863" target="_blank">Geophysical data analysis using Python </a>- Computers & Geosciences (2002)<br />
[2] - <a href="http://www.sciencedirect.com/science/article/pii/S1364815210003221" target="_blank">OOFɛ: A Python engine for automating regional and coastal ocean forecasts</a> - Environmental Modelling & Software (2011)<br />
<br />
Encontrei esse trabalho[1] dos Bascos autores do <a href="http://www.pyclimate.org/" target="_blank">PyClimate</a> (Jon Sáenz, Juan Zubillaga, Jesús Fernández), falando sobre o pacote PyClimate e a implementação de ferramentas básicas para a oceanografia.<br />
<br />
Funções de distribuição de probabilidades (PDF's), método de EOF e SVD, uso de netCDF, como trabalhar com dias julianos nas convenções do Cooperative Ocean/Atmosphere Research Data Service (COARDS), filtros multivariados e outros detalhes.<br />
<br />
O segundo trabalho[2] é um pacote de ferramentas para trabalhar com o ROMS automatizar previsões costeiras entre outras coisas. Uma demostração do que já é aplicado no Brasil está aqui: <a href="http://oceano.fis.ufba.br/oof/main/forecast.php">http://oceano.fis.ufba.br/oof/main/forecast.php</a><br />
Ainda não disponibilizado o código (não sei o porquê ainda), mas ele referencia outras ferramentas da modelagem que integram o FORTRAN, como o a implementação do modelo atmosférico Neelin-Zeng QTCM1 (feito já em python <a href="http://www.johnny-lin.com/py_pkgs/qtcm/">http://www.johnny-lin.com/py_pkgs/qtcm/</a> ).<br />
<br />
Anexos:<br />
<a href="https://docs.google.com/open?id=0B8iBXJOZmSx-UldiMUVYME1ON1U" target="_blank">[1]</a><br />
<a href="https://docs.google.com/open?id=0B8iBXJOZmSx-QkliYS1EalZUOTA" target="_blank">[2]</a><br />
<br />
Boas leituras.Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com1tag:blogger.com,1999:blog-4998210822656797146.post-62918225664032482112012-05-05T21:03:00.000-03:002012-05-05T21:03:49.135-03:00Python pra OceanógrafosAqui vão alguns módulos para trabalhar com dados oceanográficos.<br />
<br />
<b>oceans</b> # http://pypi.python.org/pypi/oceans/<br />
<b>seawater</b> # http://pypi.python.org/pypi/seawater/2.0.1<br />
<b>pyclimate</b> # http://www.pyclimate.org/<br />
<b>pyhdf</b> # http://pysclint.sourceforge.net/pyhdf/<br />
<b>pupynere</b> # para trabalhar com netcdf (http://pypi.python.org/pypi/pupynere/1.0.15 )<br />
<b>puppy</b> # http://pypi.python.org/pypi/Puppy/0.1.4<br />
<b>pydap</b> # http://pydap.org/<br />
<b>pandas</b> # http://pandas.pydata.org/<br />
<b>rpy2</b> # http://rpy.sourceforge.net/<br />
<b>f2py</b> # http://www.scipy.org/F2py<br />
<b>GDAL</b> # http://www.gdal.org/<br />
<b>pyModis</b> # https://github.com/lucadelu/pyModis<br />
<br />
E vai o link do wiki_livros, para você que já é experiente no uso... ajuda lá! =]<br />
<a href="http://pt.wikibooks.org/wiki/Python_para_ocean%C3%B3grafos">http://pt.wikibooks.org/wiki/Python_para_oceanógrafos</a>Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-44602864744477769882012-05-05T20:34:00.001-03:002012-05-05T20:34:09.590-03:00Problemas para compilar. Python.hEstava encontrando erro (Python.h: No such file or directory), para vários arquivos distintos, dependendo do pacote que ia instalar com o <b>pip</b>.<br />
<br />
Fiz uma busca e não encontrei respostas para os pacotes que apresentavam problemas na compilação.<br />
Em específico para a instalação do Pyclimate em uma das máquinas.<br />
<pre class="prettyprint">src/JDTime_wrap.c:44:20: fatal error: Python.h: Arquivo ou diretório não encontrado
</pre><br />
Para resolver isso:<br />
<pre class="prettyprint">#!/bin/bash
aptitude install python-dev
</pre>Instalará os headers e você poderá compilar normalmente.<br />
Espero que ajude!Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-6039789687148162002012-04-28T09:21:00.001-03:002012-04-28T09:21:38.731-03:00Subset MODIS-Aqua (Clorofila-a)Pegando o gancho com uma postagem <a href="http://www.trondkristiansen.com/?page_id=479">[1]</a> que já tinha guardada de um tempo e do pacote pyMODIS <a href="https://bitbucket.org/tylere/pymodis/overview">[2]</a>, resolvi fazer o meu subset de imagens MODIS em Python.<br />
<br />
Para esse código, foram testadas apenas imagens L3 de Clorofila-a, já processadas pela agência espacial norte americana (NASA).<br />
<br />
Para executar esse código, você precisará de uma imagem: <a href="http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/A20120012012031.L3m_MO_CHL_chlor_a_9km.bz2">Imagem MODIS de 9km.</a><br />
<br />
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).<br />
<br />
Um detalhe, é que abaixo segue como fazer o processo em batelada, mas para o exemplo de uma imagem, já é possível visualizar a coisa. <br />
<br />
Diferentemente do post <a href="http://www.trondkristiansen.com/?page_id=479">[1]</a>, não encontrei dificuldades em utilizar os dados. Logo não fiz uma conversão em ascii e nem utilizei o HDFview.<br />
<br />
[1] <a href="http://www.trondkristiansen.com/?page_id=479">http://www.trondkristiansen.com/?page_id=479</a><br />
[2] <a href="https://bitbucket.org/tylere/pymodis/overview">https://bitbucket.org/tylere/pymodis/overview</a><br />
<br />
<pre class="prettyprint">#/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})
</pre><br />
Resultando nessa imagem de Clorofila-a<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgWDjuLG0TVGGaKldaqKJ5L-Wu0MftQx2djepYBmsRePLhXVFmceY94_lOmcE1yTiVY7wJZZPfqvZnrZ8dQ94WsZJvUHB8KYzN-eUzI2FsI6trfvs0hRnJ76H322gDid2KD27sgFiUyH81t/s1600/chla.png" imageanchor="1" style="margin-left:1em; margin-right:1em"><img border="0" height="320" width="500" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgWDjuLG0TVGGaKldaqKJ5L-Wu0MftQx2djepYBmsRePLhXVFmceY94_lOmcE1yTiVY7wJZZPfqvZnrZ8dQ94WsZJvUHB8KYzN-eUzI2FsI6trfvs0hRnJ76H322gDid2KD27sgFiUyH81t/s320/chla.png" /></a></div><br />Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-41616717794421174722012-03-13T10:46:00.000-03:002012-03-13T11:11:00.163-03:00Ler HDF-EOS de Imagens MODIS com PythonLer HDF de Imagens MODIS com Python.<br />
Read HDF-EOS with Python<br />
<br />
Enfrentei esse problema por um tempo, por conta de como era feito o acesso às variáveis e como elas seriam carregadas.<br />
Usualmente no laboratório as pessoas utilizam o "Matlab's" para estas finalidades e para listar os atributos, dão o comando hdfinfo e para carregar o hdfread.<br />
<br />
Pelo que li, as bibliotecas incorporadas para essa finalidade, são desenvolvidas em C. Já que o HDF-EOS é uma extensão da NCSA (National Center for Supercomputing Applications). Dessa forma, o uso do HDF-EOS (Earth Observing System) é uma das formas científicas de se guardar os dados hierarquicamente.<br />
<br />
Em <b>Python</b>, existem algumas bibliotecas para essa finalidade.<br />
Testei o <b>gdal</b>(Geospatial Data Abstraction Library)[1] e o <b>pyhdf</b>[2,3]. Ambos são muito simples de serem instalados.<br />
<br />
Antes de mostrar como fazer para carregar os dados, programas no bash shell já tem uma importante função na visualização do HDF. O uso do gdalinfo é habilitado pela intalação no proprio linux é o ncdump é uma ferramenta para acesso das informações de NetCDF, mas que funcionou muito bem para o HDF-EOS.<br />
<pre class="prettyprint">#bash
# gdalinfo para listar todas os atributos
yepan@waiapi:$ gdalinfo A20100322010059.L3m_MO_CHL_chlor_a_4km.hdf
#ncdump
yepan@waiapi:$ ncdump -h A20100322010059.L3m_MO_CHL_chlor_a_4km.hdf
</pre>Para acesso ao HDF, fiz o seguinte, primeiramente com o gdal.<br />
<pre class="prettyprint">#!/usr/bin/env python
import gdal
filename = 'A20100322010059.L3m_MO_CHL_chlor_a_4km.hdf'
ds = gdal.Open(filename)
sds_md = ds.GetMetadata() #lista os atributos com os valores, guardando
#em um dicionário.
for i in sds_md.iterkeys():
print i,sds_md[i]
data = ds.ReadAsArray()
data
Out:
array([[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
...,
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.]], dtype=float32)
</pre>Se o HDF tiver algum "Subdataset" é possível verificar desse jeito:<br />
<pre class="prettyprint">#!/usr/bin/env python
"""
neste caso é um HDF de coordenadas de gelo
http://www.iup.uni-bremen.de:8084/amsredata/asi_daygrid_swath/l1a/s6250/grid_coordinates/LongitudeLatitudeGrid-s6250-Antarctic.hdf
"""
import gdal
import gdalnumeric
filename = 'LongitudeLatitudeGrid-s6250-Antarctic.hdf'
ds = gdal.Open(filename)
sds_md = ds.GetMetadata('SUBDATASETS')
{'SUBDATASET_1_DESC': '[1328x1264] Longitudes (32-bit floating-point)',
'SUBDATASET_1_NAME': 'HDF4_SDS:UNKNOWN:"LongitudeLatitudeGrid-s6250-Antarctic.hdf":0',
'SUBDATASET_2_DESC': '[1328x1264] Latitudes (32-bit floating-point)',
'SUBDATASET_2_NAME': 'HDF4_SDS:UNKNOWN:"LongitudeLatitudeGrid-s6250-Antarctic.hdf":1'}
datakeys = {}
datasets = ['SUBDATASET_1','SUBDATASET_2']
datanames = ['Longitudes', 'Latitudes']
for (j,i) in enumerate(datasets):
this = {}
this['name'] = sds_md[i + '_NAME']
this['description'] = sds_md[i + '_DESC']
this['data'] = gdalnumeric.LoadFile(this['name'])
datakeys[datanames[j]] = this.copy()
#para ver os dados:
datakeys['Longitudes']['data']
datakeys['Latitudes']['data']
</pre>Quanto ao pyhdf, não me aprofundei ainda nos estudos, mas posso adiantar que ele funcionou para abrir a imagem do MODIS e as coordenadas de gelo também, mas de forma muito mais simples. Vamos lá!<br />
<pre class="prettyprint">#!/usr/bin/env python
from pyhdf.SD import *
filename = 'A20100322010059.L3m_MO_CHL_chlor_a_4km.hdf'
ds = SD(filename)
ds.attributes() #mostra os atributos
ds.datasets()
Out : {'l3m_data': (('fakeDim0', 'fakeDim1'), (4320, 8640), 5, 0)}
data = ds.select('l3m_data')
data[:]
Out :
array([[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
...,
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.],
[-32767., -32767., -32767., ..., -32767., -32767., -32767.]], dtype=float32)
</pre>Para acessar as coordenadas de dados de gelo:<br />
<pre class="prettyprint">#!/usr/bin/env python
from pyhdf.SD import *
filename = 'LongitudeLatitudeGrid-s6250-Antarctic.hdf'
ds = SD(filename)
ds.datasets()
Out:
{'Latitudes': (('fakeDim2', 'fakeDim3'), (1328, 1264), 5, 1),
'Longitudes': (('fakeDim0', 'fakeDim1'), (1328, 1264), 5, 0)}
lat = ds.select('Latitudes')
lat.dimensions()
Out: {'fakeDim2': 1328, 'fakeDim3': 1264}
#para ver os dados lat[:]
#Existe também a função "get", que permite pegar uma porção do dado, para datasets multidimensionais
</pre><br />
[1] <a href="http://www.gdal.org/">http://www.gdal.org/</a><br />
[2] <a href="http://hdfeos.org/software/pyhdf.php">http://hdfeos.org/software/pyhdf.php</a><br />
[3] <a href="http://pysclint.sourceforge.net/pyhdf/">http://pysclint.sourceforge.net/pyhdf/</a>Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0tag:blogger.com,1999:blog-4998210822656797146.post-315895461153900792012-02-08T18:11:00.039-02:002012-02-16T19:34:55.709-02:00Ler dados binários de gelo marinho com PythonEstava com erros ao ler arquivos binários de médias mensais de gelo marinho, conseguidos através do site:<br />
<a href="ftp://sidads.colorado.edu/DATASETS/seaice/polar-stereo/bootstrap/final-gsfc/south/monthly">ftp://sidads.colorado.edu/DATASETS/seaice/polar-stereo/bootstrap/final-gsfc/south/monthly</a><br />
<span style="background-color: white;"><br />
</span><br />
<span style="background-color: white;"></span><br />
<pre class="prettyprint"><pre class="brush: python">import numpy as np
import matplotlib.pyplot as plt
from glob import glob
from scipy.io import loadmat
# Lista de arquivos binários
d = sorted(glob('bt*'))
# pegando somente um dos arquivos para leitura
f = open(d[0],'rb')
map1 = np.fromfile(f, dtype=np.uint16).reshape(316,332)
f.close()
im = map1.T/10.
SI3[:,:,i] = im
# Exemplo de linha do arquivo. Essa é a menor delas.
len(a) = 70.
a = f.readlines()[1]
Out [2]: '\x03\x18\x03,\x03C\x03G\x03R\x03i\x03d\x03q\x03\x88\x03\x87\x03\x9a\</pre>
<pre class="brush: python">x03\xbd\x03\xbd\x03\xb8\x03\xbd\x03\xc0\x03\xbf\x03\xc5\x03\xb5\x03\xa2\x03\x85</pre>
<pre class="brush: python">\x03[\x03N\x03B\x03\x1a\x03\xeb\x02\xc4\x02\x95\x02\x94\x02p\x02w\x02X\x02H\x02.\x02\n'</pre>
</pre>
<span style="background-color: white;"><br />
</span><br />
<span style="background-color: white;">Sabia que estava errando alguma coisa, pois ao plotar em um mapa, esses dados, eles ficavam assim:</span><br />
<span style="background-color: white;"></span><br />
<pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiLEdHVUiZGuoHtZSUJEz_3k_guQ2z1x_9z5H81HGYngDFI0zwnTexBmCE_r7pRImUKxqFSout3I8juxb1-pyYkEyP4O-kwvAUrSaNs3EW9PmtcUAsszocWgxr1orFqLsjx_u_NGVUSzaIq/s1600/testepy.png" imageanchor="1" style="margin-left: auto; margin-right: auto; white-space: pre;"><span style="color: black;"><img border="0" height="481" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiLEdHVUiZGuoHtZSUJEz_3k_guQ2z1x_9z5H81HGYngDFI0zwnTexBmCE_r7pRImUKxqFSout3I8juxb1-pyYkEyP4O-kwvAUrSaNs3EW9PmtcUAsszocWgxr1orFqLsjx_u_NGVUSzaIq/s640/testepy.png" width="640" /></span></a></td></tr>
<tr><td class="tr-caption" style="font-size: 10px;"><span style="font-size: small;"><span style="white-space: pre;">Figura tosca dos dados de Gelo.</span></span></td></tr>
</tbody></table>
</pre>
<span style="background-color: white;">Depois de muito quebrar a cabeça, analisando pedaços da matriz e comparando com um arquivo que tinha certo, vi que parecia estar faltando algo da codificação.</span><br />
<span style="background-color: white;">De fato, isso se dava ao reshape que fazia ao final de fazer a leitura com o <b>np.fromfile</b>.</span><br />
<span style="background-color: white;">O que acontece, é que por falta de informações sobre esses dados, não sabia como eles estava encapsulados. A única informação que tinha é que deveriam ser lidos de 2 em 2 bytes ( o que é justificado pelo <b>dtype = np.uint16</b> ).</span><br />
<span style="background-color: white;">Mas, como sabia que estava errado, varri a rede buscando alguma informação mais efetiva sobre isso.</span><br />
<span style="background-color: white;">Muito falaram sobre o módulo struct [<a href="http://stackoverflow.com/questions/2865996/reading-a-binary-file-in-python-into-a-struct" target="_blank">link da dica 1</a>], mas que não consegui progredir por esse caminho.</span><br />
<span style="background-color: white;">Comecei a verificar detalhadamente cada parte do código, que tinha, e acabei descobrindo de onde vinha o problema! =]). Tudo estava centrado na questão de que o NumPy faz um reshape nas matrizes de acordo com referências em C e ou FORTRAN.</span><br />
<br />
<pre class="prettyprint"><span style="background-color: white;"> reshape(a, newshape, order='C'):
"""
Gives a new shape to an array without changing its data.
Parameters
----------
a : array_like
Array to be reshaped.
newshape : int or tuple of ints
The new shape should be compatible with the original shape. If
an integer, then the result will be a 1-D array of that length.
One shape dimension can be -1. In this case, the value is inferred
from the length of the array and remaining dimensions.
order : {'C', 'F'}, optional
Determines whether the array data should be viewed as in C
(row-major) order or FORTRAN (column-major) order.</span>
<span style="background-color: white;"></span></pre>
<span style="background-color: white;">Li isso tudo muito rápido e fui testar... Coloquei um order=True e FUNCIONOU! (Acho que por sorte mesmo.) Agora está mudado para order='F', onde o meu dado "passa bem, sem remédios nem cirurgia".</span><br />
<pre style="background-attachment: initial; background-clip: initial; background-color: white; background-image: initial; background-origin: initial;"><pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"><pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"><div style="font-family: 'Times New Roman'; white-space: normal;">
A saída que obtinha era:</div>
<pre class="prettyprint"><div style="font-family: 'Times New Roman'; white-space: normal;">
<span style="background-color: white;"></span></div>
In [5]: map1[:,123]
Out [6]:
array([1200, 1200, 1200, 1200, 1200, 1200, 712, 870, 928, 958, 957,
974, 969, 983, 984, 979, 976, 978, 976, 977, 978, 979,
982, 986, 988, 988, 993, 988, 990, 987, 997, 993, 991,
995, 997, 990, 986, 992, 989, 993, 990, 975, 985, 979,
979, 988, 976, 984, 970, 958, 962, 956, 941, 907, 841,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 829, 845,
774, 583, 486, 474, 557, 681, 774, 824, 840, 857, 861,
862, 848, 838, 821, 784, 747, 689, 605, 522, 529, 500,
468, 469, 470, 437, 414, 361, 360, 288, 236, 205, 151,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 191,
326, 590, 832, 914, 915, 950, 935, 1200, 951, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 695, 890, 945, 973, 960, 956, 972,
976, 976])</pre>
<pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"><pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"><pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"></pre>
<pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"><span style="background-color: white; font-family: 'Times New Roman'; white-space: normal;">Essa saída, é muito diferente da saída correta, que se consegue com o processamento, incluindo o order='F'.</span></pre>
<pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"><pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"><pre class="prettyprint"><pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;">In [7]: map1[:,123]
Out [8]:
array([1200, 1200, 1200, 1200, 0, 1200, 1200, 1200, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 104, 228,
301, 389, 1200, 519, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 651, 703, 748, 742, 739, 777, 831,
915, 958, 965, 978, 988, 984, 988, 991, 989, 986, 986,
986, 987, 983, 990, 989, 980, 985, 986, 980, 980, 974,
971, 968, 980, 988, 983, 991, 991, 992, 990, 982, 978,
973, 987, 989, 993, 987, 980, 977, 980, 982, 981, 984,
976, 978, 965, 948, 937, 920, 881, 846, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200, 1200,
1200, 852, 845, 852, 806, 657, 630, 554, 573, 621, 682,
738, 778, 801, 804, 764, 711, 677, 662, 650, 628, 581,
544, 517, 483, 459, 422, 395, 394, 400, 409, 380, 373,
352, 279, 186, 111, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0], dtype=uint16)</pre>
</pre>
</pre>
</pre>
</pre>
</pre>
E aqui a figura tosca, mas com os dados passando bem. (Sea Ice, binary data)</pre>
</pre>
<pre style="background-attachment: initial; background-clip: initial; background-color: white; background-image: initial; background-origin: initial;"><pre style="background-attachment: initial; background-clip: initial; background-image: initial; background-origin: initial;"><div>
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="float: left; margin-right: 1em; text-align: left;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgUH5B4ETmQAjvRlFn2DtNWMc3MhiACaNhHBiKSJAA6DGQN4n7nK5XlGQVVaFgH-lC_IhzG2PVh8aS1_J9RlpCMA37eY8JK6wW39ZBAj-IYsSd9dj_3g8tvcfqDX4tCoyL2sx4fzY2Zwnlx/s1600/teste2.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="482" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgUH5B4ETmQAjvRlFn2DtNWMc3MhiACaNhHBiKSJAA6DGQN4n7nK5XlGQVVaFgH-lC_IhzG2PVh8aS1_J9RlpCMA37eY8JK6wW39ZBAj-IYsSd9dj_3g8tvcfqDX4tCoyL2sx4fzY2Zwnlx/s640/teste2.png" width="640" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Figura tosca, mas com os dados funcionando.</td></tr>
</tbody></table>
<span style="font-family: 'Times New Roman';"><span style="white-space: normal;">
</span></span></div>
</pre>
</pre>
</pre>Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com1tag:blogger.com,1999:blog-4998210822656797146.post-35000756027095923962012-02-08T15:54:00.000-02:002012-03-17T20:29:58.585-03:00E com Python digo:<pre style="background: #ffffff; color: black;"><span style="color: maroon; font-weight: bold;">print</span> 'Oi Mundo<span style="color: #808030;">!</span>' </pre>
<pre style="background: none repeat scroll 0% 0% rgb(255, 255, 255); color: black;"> </pre>
<br />
<pre style="background: none repeat scroll 0% 0% rgb(255, 255, 255); color: black;"> </pre>Arnaldinhohttp://www.blogger.com/profile/09766843820336991926noreply@blogger.com0