In [1]:
#Team Big Coders- Elke and Joy
#import modules
from osgeo import gdal
from osgeo.gdalconst import GA_ReadOnly, GDT_Float32
import numpy as np
# open file and print info about the file
# the "../" refers to the parent directory of my working directory
!wget https://www.dropbox.com/s/zb7nrla6fqi1mq4/LC81980242014260-SC20150123044700.tar.gz?dl=1 -O LC81980242014260-SC20150123044700.tar.gz
import tarfile
tar = tarfile.open("LC81980242014260-SC20150123044700.tar.gz", "r:gz")
tar.extractall(path="data")
tar.close()
filename4 = './data/LC81980242014260LGN00_sr_band4.tif'
filename5 = './data/LC81980242014260LGN00_sr_band5.tif'
dataSource4 = gdal.Open(filename4, GA_ReadOnly)
dataSource5 = gdal.Open(filename5, GA_ReadOnly)
In [2]:
print "\nInformation about " + filename4 + filename5
print "Driver4: ", dataSource4.GetDriver().ShortName,"/", \
dataSource4.GetDriver().LongName
print "Driver4: ", dataSource5.GetDriver().ShortName,"/", \
dataSource5.GetDriver().LongName
print "Size4 is ",dataSource4.RasterXSize,"x",dataSource5.RasterYSize, \
'x',dataSource4.RasterCount
print "Size5 is ",dataSource5.RasterXSize,"x",dataSource5.RasterYSize, \
'x',dataSource5.RasterCount
In [3]:
print '\nProjection is: ', dataSource4.GetProjection()
print '\nProjection is: ', dataSource5.GetProjection()
In [3]:
In [3]:
In [4]:
print "\nInformation about the location of the image and the pixel size:"
geotransform4 = dataSource4.GetGeoTransform()
geotransform5 = dataSource5.GetGeoTransform()
print geotransform4
print geotransform5
In [5]:
if not geotransform4 is None:
print 'Origin = (',geotransform4[0], ',',geotransform4[3],')'
print 'Pixel Size = (',geotransform4[1], ',',geotransform4[5],')'
In [6]:
if not geotransform5 is None:
print 'Origin = (',geotransform5[0], ',',geotransform5[3],')'
print 'Pixel Size = (',geotransform5[1], ',',geotransform5[5],')'
In [7]:
# Read data into an array
band4Arr = dataSource4.ReadAsArray(0,0,dataSource4.RasterXSize, dataSource4.RasterYSize)
band5Arr = dataSource5.ReadAsArray(0,0,dataSource5.RasterXSize, dataSource5.RasterYSize)
print type(band4Arr)
print type(band5Arr)
In [8]:
# set the data type
band4Arr=band4Arr.astype(np.float32)
band5Arr=band5Arr.astype(np.float32)
In [9]:
# Derive the NDVI
mask = np.greater(band4Arr+band5Arr,0)
# set np.errstate to avoid warning of invalid values (i.e. NaN values) in the divide
with np.errstate(invalid='ignore') and np.errstate(divide='ignore'):
ndwi = np.choose(mask,(-99,(band4Arr-band5Arr)/(band4Arr+band5Arr)))
print "NDWI min and max values", ndwi.min(), ndwi.max()
# Check the real minimum value
print ndwi[ndwi>-99].min()
In [10]:
# Write the result to disk
driver = gdal.GetDriverByName('GTiff')
outDataSet=driver.Create('./data/ndwi.tif', dataSource4.RasterXSize, dataSource4.RasterYSize, 1, GDT_Float32)
outBand = outDataSet.GetRasterBand(1)
outBand.WriteArray(ndwi,0,0)
outBand.SetNoDataValue(-99)
Out[10]:
In [10]:
In [10]:
In [11]:
# set the projection and extent information of the dataset
outDataSet.SetProjection(dataSource4.GetProjection())
outDataSet.SetGeoTransform(dataSource4.GetGeoTransform())
# Finally let's save it... or like in the OGR example flush it
outDataSet.FlushCache()
In [12]:
!gdalinfo ./data/ndwi.tif
In [13]:
# via the Shell
!gdalwarp -t_srs "EPSG:4326" ./data/ndwi.tif ./data/ndwi_ll.tif
In [14]:
# Let's check what the result is
!gdalinfo ./data/ndwi_ll.tif
In [15]:
# Notebook magic to select the plotting method
# Change to inline to plot within this notebook
#%matplotlib inline
%matplotlib inline
from osgeo import gdal
import matplotlib.pyplot as plt
In [16]:
# Open image
dsll = gdal.Open("./data/ndwi_ll.tif")
# Read raster data
ndvi = dsll.ReadAsArray(0, 0, dsll.RasterXSize, dsll.RasterYSize)
# Now plot the raster data using gist_earth palette
plt.imshow(ndwi, interpolation='nearest', vmin=0, cmap=plt.cm.gist_earth)
plt.show()
dsll = None
In [16]: