notebook.community

Edit and run

Lesson 13, NDWI and reprojection in python
Team Rython
Dainius Masiliunas & Tim Weerman 20 January 2016
Apache license 2.0



In [2]:

    
# Import modules
import urllib
import os
from osgeo import gdal
from osgeo.gdalconst import GA_ReadOnly, GDT_Float32
import numpy as np

# Change working directory
os.getcwd()
#######################################################
#os.chdir("Yourworkingdirectory/geoscripting/Lesson13")
#######################################################









    Out[2]:





'/home/tim/geoscripting/Lesson13'

The function download() lets you download the two necessary bands from dropbox. If the files are already present in your data directory it won't download the files from dropbox.



In [4]:

    
filelist = ["LC81980242014260LGN00_sr_band4.tif", "LC81980242014260LGN00_sr_band5.tif"]

def download():
    if not os.path.isfile("data/" +filelist[0]):
        # Download the file
        print "The file is downloading, this can take a while"
        urllib.urlretrieve("https://www.dropbox.com/s/zb7nrla6fqi1mq4/LC81980242014260-SC20150123044700.tar.gz?dl=1", "data/Netherlands.tar.gz")
        # Untar the file
        for filename in filelist:
           os.system("tar -zxvf data/Netherlands.tar.gz -C data " +filename)

The calculateNDWI() calculates the NDWI from the 2 given files. The output file is written to "output".



In [5]:

    
def calculateNDWI(filename4, filename5, output):
    # Set
    dataSource4 = gdal.Open("data/" +filename4, GA_ReadOnly)
    dataSource5 = gdal.Open("data/" +filename5, GA_ReadOnly)

    # Save bands 4 & 5
    band4 = dataSource4.GetRasterBand(1).ReadAsArray(0,0,dataSource4.RasterXSize, dataSource4.RasterYSize).astype(np.float32)
    band5 = dataSource5.GetRasterBand(1).ReadAsArray(0,0,dataSource5.RasterXSize, dataSource5.RasterYSize).astype(np.float32)

    # Derive NDWI
    # Set np.errstate to avoid warning of invalid values (i.e. NaN values) in the divide 
    with np.errstate(divide='ignore'):
        ndwi = ((band4-band5)/(band4+band5))
        
    ndwi[np.isnan(ndwi)] = -99
    # Save the image
    driver = gdal.GetDriverByName('GTiff')
    outDataSet=driver.Create(output, dataSource4.RasterXSize, dataSource4.RasterYSize, 1, GDT_Float32)
    outBand = outDataSet.GetRasterBand(1)
    outBand.WriteArray(ndwi,0,0)
    outBand.SetNoDataValue(-99)

    # Set the projection and extent information of the dataset
    outDataSet.SetProjection(dataSource4.GetProjection())
    outDataSet.SetGeoTransform(dataSource4.GetGeoTransform())

    # Current projection is...
    print "\nInformation about " +output
    print '\nProjection is: ', outDataSet.GetProjection()

    outBand.FlushCache()
    outDataSet.FlushCache()

ReprojectNDWI() reprojects the file to a projection given by the user.



In [6]:

    
def reprojectNDWI(projinput, projoutput, projection):
    # Reprojection
    os.system("gdalwarp -t_srs " +projection +" " +projinput +" " +projoutput)

    # Current projection is...
    reprojected = gdal.Open(projoutput, GA_ReadOnly)
    print "\nInformation about the new reprojected file " +projoutput
    print '\nProjection is: ', reprojected.GetProjection()

Here is an example of how to run the script correctly.



In [7]:

    
download()
calculateNDWI(filelist[0], filelist[1], "data/ndwi.tif")
reprojectNDWI("data/ndwi.tif", "data/ndwi_ll.tif", "'EPSG:4326'")









    



Information about data/ndwi.tif

Projection is:  PROJCS["WGS 84 / UTM zone 31N",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",3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32631"]]

Information about the new reprojected file data/ndwi_ll.tif

Projection is:  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"]]



In [8]:

    
# End of the script
print "This is the end of the script, your image has been created."









    



This is the end of the script, your image has been created.