Latitude e Longitude - Ice cover - NSIDC

Aproveitando 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.

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. ("<" Little Endian; ">" 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".

Dessa forma, baixa-se o dado de latitude para o hemisferio sul [aqui] e divide-se por 100000 para ter os graus decimais. Mesmo procedimento para a Longitude.

import numpy as np
f = open('pss25lats_v2.dat','rb')
lat = np.fromfile(f, dtype='<i4')/100000.
lat = lat.reshape(332,316)
f.close()

[1] http://nsidc.org/data/polar_stereo/tools_geo_pixel.html

[2] http://docs.scipy.org/doc/numpy/user/basics.byteswapping.html

Loops 'tradicionais'

Meu camarada Irber, 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".

Buscando, dessa forma, um loop tradicionalista, porém mais pythonico (eu acho).

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

Sendo o jeito um tanto mais "conservador" e talvez menos pythonico:

for i in xrange(len(valores)):
    print i, valores[i]

Out[4]:
0 1
1 3
2 5
3 7
4 9

Até mais.

Livros II - Oceanic Sciences and Hydrology

Um livro interessante, sobre Python aplicado às ciências atmosféricas e oceânicas.
A Hands-On Introduction to Using Python in the Atmospheric and Oceanic Sciences
Outro livro bacana é sobre a utilizacão do Python em hidrologia.
Python in Hydrology
Boas leituras.

Livro: Estudos Oceanográficos

Está disponível para download o livro:
Estudos Oceanográficos - do instrumental ao prático
Organizado pelo Professor Danilo Calazans - FURG.

Outros livros estão disponíveis no site do Ciências do Mar.

Boas leituras!

Papers - Oceanography and Python

Dando parabéns aos Oceanógrafos pelo "dia do Oceanógrafo", seguem algumas leituras do dia de hoje.

[1] - Geophysical data analysis using Python - Computers & Geosciences (2002)
[2] - OOFɛ: A Python engine for automating regional and coastal ocean forecasts - Environmental Modelling & Software (2011)

Encontrei esse trabalho[1] dos Bascos autores do PyClimate (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.

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.

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: http://oceano.fis.ufba.br/oof/main/forecast.php
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 http://www.johnny-lin.com/py_pkgs/qtcm/ ).

Anexos:
[1]
[2]

Boas leituras.

Python pra Oceanógrafos

Aqui vão alguns módulos para trabalhar com dados oceanográficos.

oceans # http://pypi.python.org/pypi/oceans/
seawater # http://pypi.python.org/pypi/seawater/2.0.1
pyclimate # http://www.pyclimate.org/
pyhdf # http://pysclint.sourceforge.net/pyhdf/
pupynere # para trabalhar com netcdf (http://pypi.python.org/pypi/pupynere/1.0.15 )
puppy # http://pypi.python.org/pypi/Puppy/0.1.4
pydap # http://pydap.org/
pandas # http://pandas.pydata.org/
rpy2 # http://rpy.sourceforge.net/
f2py # http://www.scipy.org/F2py
GDAL # http://www.gdal.org/
pyModis # https://github.com/lucadelu/pyModis

E vai o link do wiki_livros, para você que já é experiente no uso... ajuda lá! =]
http://pt.wikibooks.org/wiki/Python_para_oceanógrafos

Problemas para compilar. Python.h

Estava encontrando erro (Python.h: No such file or directory), para vários arquivos distintos, dependendo do pacote que ia instalar com o pip.

Fiz uma busca e não encontrei respostas para os pacotes que apresentavam problemas na compilação.
Em específico para a instalação do Pyclimate em uma das máquinas.

src/JDTime_wrap.c:44:20: fatal error: Python.h: Arquivo ou diretório não encontrado

Para resolver isso:

#!/bin/bash
aptitude install python-dev

Instalará os headers e você poderá compilar normalmente.
Espero que ajude!

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

Ler HDF-EOS de Imagens MODIS com Python

Ler HDF de Imagens MODIS com Python.
Read HDF-EOS with Python

Enfrentei esse problema por um tempo, por conta de como era feito o acesso às variáveis e como elas seriam carregadas.
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.

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.

Em Python, existem algumas bibliotecas para essa finalidade.
Testei o gdal(Geospatial Data Abstraction Library)[1] e o pyhdf[2,3]. Ambos são muito simples de serem instalados.

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.

#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

Para acesso ao HDF, fiz o seguinte, primeiramente com o gdal.

#!/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)

Se o HDF tiver algum "Subdataset" é possível verificar desse jeito:

#!/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']

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á!

#!/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)

Para acessar as coordenadas de dados de gelo:

#!/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

[1] http://www.gdal.org/
[2] http://hdfeos.org/software/pyhdf.php
[3] http://pysclint.sourceforge.net/pyhdf/

Ler dados binários de gelo marinho com Python

Estava com erros ao ler arquivos binários de médias mensais de gelo marinho, conseguidos através do site:
ftp://sidads.colorado.edu/DATASETS/seaice/polar-stereo/bootstrap/final-gsfc/south/monthly

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\
x03\xbd\x03\xbd\x03\xb8\x03\xbd\x03\xc0\x03\xbf\x03\xc5\x03\xb5\x03\xa2\x03\x85
\x03[\x03N\x03B\x03\x1a\x03\xeb\x02\xc4\x02\x95\x02\x94\x02p\x02w\x02X\x02H\x02.\x02\n'

Sabia que estava errando alguma coisa, pois ao plotar em um mapa, esses dados, eles ficavam assim:

Figura tosca dos dados de Gelo.

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.
De fato, isso se dava ao reshape que fazia ao final de fazer a leitura com o np.fromfile.
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 dtype = np.uint16 ).
Mas, como sabia que estava errado, varri a rede buscando alguma informação mais efetiva sobre isso.
Muito falaram sobre o módulo struct [link da dica 1], mas que não consegui progredir por esse caminho.
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.

 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.

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".

A saída que obtinha era:


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])

Essa saída, é muito diferente da saída correta, que se consegue com o processamento, incluindo o order='F'.
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)





E aqui a figura tosca, mas com os dados passando bem. (Sea Ice, binary data)


Figura tosca, mas com os dados funcionando.

E com Python digo:

print 'Oi Mundo!'