Leer HDF4 data (productos modis) utilizando PyTHON con la librería GDAL

En un post anterior ya se había señalado el soporte GDAL para la lectura de archivos con formato hdf4/hdf5 de productos MODIS. Estos productos están constituidos, a su vez, por una serie de sub datasets que pueden ser accedidos por separado, por ejemplo, en un script de PyTHON.

El objetivo de este post es acceder a los dataset de un producto modis ‘MOD13Q1’ (Índices de Vegetación), específicamente el primero (valores no escalados de NDVI), para luego imprimir las coordenadas (x,y) de un value específico (256) escogido arbitrariamente. El código propuesto es el siguiente:

from osgeo import gdal
import struct

nameraster = "MOD13Q1.A2005193.h10v08.005.2008215173619.hdf"

hdf_file = gdal.Open(nameraster)

subDatasets = hdf_file.GetSubDatasets()

dataset = gdal.Open(subDatasets[0][0])
geotransform = dataset.GetGeoTransform()
band = dataset.GetRasterBand(1)

fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 
            'Int32':'i', 'Float32':'f', 'Float64':'d'}

print "rows = %d columns = %d" % (band.YSize, band.XSize)

BandType = gdal.GetDataTypeName(band.DataType)

print "Data type = ", BandType

print "Executing with %s" % nameraster

print "test_value = 256"

X = geotransform[0] #x coordinate
Y = geotransform[3] #y coordinate

for y in range(band.YSize):

    scanline = band.ReadRaster(0, 
                               y, 
                               band.XSize, 
                               1, 
                               band.XSize, 
                               1, 
                               band.DataType)

    values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)

    for value in values:

        if(value == 256):       
            print "%.4f %.4f %.2f" % (X, Y, value)
        X += geotransform[1] #x pixel size
    X = geotransform[0]
    Y += geotransform[5] #y pixel size

dataset = None

cuya ejecución en la bash console produce este resultado:

rows = 4800 columns = 4800
Data type =  Int16
Executing with MOD13Q1.A2005193.h10v08.005.2008215173619.hdf
test_value = 256
-8188820.6083 1048476.6775 256.00
-8188820.6083 1048245.0211 256.00
-7870988.0847 1031797.4197 256.00
-8567115.4413 1028554.2307 256.00
-8526343.9223 981296.3336 256.00
-8739236.1155 975968.2374 256.00
-8180249.3230 972493.3920 256.00
-8690819.9366 957899.0414 256.00
-8771668.0057 954192.5397 256.00
-8703329.3800 915274.2715 256.00
-8585416.2936 915042.6151 256.00
-7988669.5147 722999.4941 256.00
-7986584.6075 717208.0852 256.00
-7994924.2364 699833.8583 256.00
-7992839.3292 697748.9511 256.00
-7992607.6728 697748.9511 256.00
-7993534.2983 696822.3257 256.00
-8177932.7594 553658.6963 256.00
-8177701.1031 553658.6963 256.00
-7904114.9440 532114.6549 256.00
-8148975.7146 505474.1737 256.00
-8148744.0583 505242.5174 256.00
-8595145.8607 460532.8402 256.00
-8595145.8607 460301.1839 256.00
-8410979.0558 408873.4723 256.00
-8410747.3995 408873.4723 256.00
-8410979.0558 408641.8160 256.00
-8410747.3995 408641.8160 256.00
-8408199.1795 407715.1905 256.00
-8429743.2209 372040.1114 256.00
-8441094.3824 325477.1834 256.00
-8442484.3206 325013.8706 256.00
-8658156.3901 293045.2932 256.00
-8764486.6586 153819.8219 256.00
-8764255.0022 153819.8219 256.00
-8764486.6586 153588.1655 256.00
-8764255.0022 153588.1655 256.00
-8317621.5435 146638.4748 256.00
-8377620.5403 85481.1962 256.00
Esta entrada fue publicada en Código Python, GDAL, modis, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s