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