GDAL

In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib

GDAL (Geospatial Data Abstraction Library) possui um módulo Python:

In [2]:
import gdal

Podemos utilizá-lo para manipular arquivos GeoTIFF:

In [3]:
bahia_data = gdal.Open('data/2010-07-03T142001_RE4_3A-NAC_6683686_113276.tif')

Informações do dataset

In [4]:
bahia_data.GetMetadata()
Out[4]:
{'AREA_OR_POINT': 'Area',
 'DATE_OF_CREATION': '06_54 29-Jun-11',
 'DATE_OF_UPDATE': '06_54 29-Jun-11',
 'FILE_ID': 'IMG_0000-000-374-194__imagery.pix',
 'GENERATING_FACILITY': 'RapidEye AG',
 'SOFTWARE': 'DPS 8.1'}
In [5]:
bahia_data.GetDriver().ShortName
Out[5]:
'GTiff'
In [6]:
bahia_data.GetDriver().LongName
Out[6]:
'GeoTIFF'
In [7]:
print bahia_data.RasterXSize
print bahia_data.RasterYSize
print bahia_data.RasterCount
5000
5000
5
In [8]:
bahia_data.GetProjection()
Out[8]:
'PROJCS["WGS 84 / UTM zone 23S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32723"]]'
In [9]:
geotransform = bahia_data.GetGeoTransform()
print 'Origin', (geotransform[0], geotransform[3])
print 'Pixel Size', (geotransform[1], geotransform[5])
Origin (403500.0, 8656500.0)
Pixel Size (5.0, -5.0)

Dados como um array NumPy

In [10]:
B = bahia_data.ReadAsArray()
In [71]:
B.shape
Out[71]:
(5, 5000, 5000)

Da documentação do RapidEye, sabemos o significado das 5 bandas:

In [12]:
band_name = {0:'Blue', 1:'Green', 2:'Red', 3:'Red Edge', 4:'Near-Infrared'}

E seus comprimentods de onda:

In [14]:
band_range = {0:(440,510), 1:(520,590), 2:(630,685), 3:(690,730), 4:(760,850)}
In [15]:
bands = range(5)
for b in bands:
    print 'Banda %d, %s, %s' % (b, band_name[b], band_range[b])
Banda 0, Blue, (440, 510)
Banda 1, Green, (520, 590)
Banda 2, Red, (630, 685)
Banda 3, Red Edge, (690, 730)
Banda 4, Near-Infrared, (760, 850)
In [16]:
figsize(20,30)
for b in bands:
    subplot(3,2,b+1)
    imshow(B[b], interpolation='nearest', cmap=cm.hot)
    colorbar()
    title(band_name[b])