source: trunk/zoo-project/zoo-services/gdal/ndvi/cgi-env/ndvi.py @ 781

Last change on this file since 781 was 781, checked in by djay, 4 years ago

Fix issue #140

  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-python
File size: 4.3 KB
Line 
1#
2# Python Scripting for ArcGIS
3# by Dr Peter Bunting is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License -
4# See more at:
5# http://www.landmap.ac.uk/index.php/Learning-Materials/Python-Scripting/9.4-Calculate-NDVI-using-GDAL
6#
7
8import sys, os, struct
9
10import osgeo.gdal as gdal
11
12# Calculate and output NDVI from raster bands
13def ExtractNDVI(conf, inputs, outputs):
14   
15    # Open the input dataset
16    gdal.FileFromMemBuffer('/vsimem//temp.tif', inputs["raster"]["value"])
17    dataset = gdal.Open( '/vsimem//temp.tif')
18    if dataset is None:
19        conf["lenv"]["message"]="The dataset could not be openned properly"
20        return 4
21
22    # Create the output dataset
23    driver = gdal.GetDriverByName( "GTiff" )
24    # Get the spatial information from the input file
25    geoTransform=None
26    geoProjection=None
27    print >> sys.stderr,dir(dataset) 
28    try:
29        geoTransform = dataset.GetGeoTransform()
30    except:
31        print >> sys.stderr, "Unable to load geotransform"
32    try:
33        geoProjection = dataset.GetProjection()
34    except:
35        print >> sys.stderr, "Unable to load projection"
36
37    # Create an output file of the same size as the inputted
38    # image but with only 1 output image band.
39    newDataset = driver.Create("/vsimem//output"+conf["lenv"]["sid"], \
40                               dataset.RasterXSize, dataset.RasterYSize,1, \
41                               gdal.GDT_Float32)
42    # Set spatial information of the new image.
43    if geoTransform:
44        newDataset.SetGeoTransform(geoTransform)
45    if geoProjection:
46        newDataset.SetProjection(geoProjection)
47    if newDataset is None:
48        conf["lenv"]["message"]='Could not create output image'
49        return 4
50   
51    # Get the RED and NIR image bands of the image
52    red_id=int(inputs["red"]["value"])
53    nir_id=int(inputs["nir"]["value"])
54    red_band = dataset.GetRasterBand(red_id) # RED BAND
55    nir_band = dataset.GetRasterBand(nir_id) # NIR BAND
56
57    # Loop through each line in turn.
58    numLines = red_band.YSize
59    for line in range(numLines):
60        # Define variable for output line.
61        outputLine = ''
62        # Read in data for the current line from the
63        # image band representing the red wavelength
64        red_scanline = red_band.ReadRaster( 0, line, red_band.XSize, 1, \
65                                          red_band.XSize, 1, gdal.GDT_Float32 )
66        # Unpack the line of data to be read as floating point data
67        red_tuple = struct.unpack('f' * red_band.XSize, red_scanline)
68           
69        # Read in data for the current line from the
70        # image band representing the NIR wavelength
71        nir_scanline = nir_band.ReadRaster( 0, line, nir_band.XSize, 1, \
72                                          nir_band.XSize, 1, gdal.GDT_Float32 )
73        # Unpack the line of data to be read as floating point data
74        nir_tuple = struct.unpack('f' * nir_band.XSize, nir_scanline)
75
76        # Loop through the columns within the image
77        for i in range(len(red_tuple)):
78            # Calculate the NDVI for the current pixel.
79            ndvi_lower = (nir_tuple[i] + red_tuple[i])
80            ndvi_upper = (nir_tuple[i] - red_tuple[i])
81            ndvi = 0
82            # Becareful of zero divide
83            if ndvi_lower == 0:
84                ndvi = 0
85            else:
86                ndvi = ndvi_upper/ndvi_lower
87            # Add the current pixel to the output line
88            outputLine = outputLine + struct.pack('f', ndvi)
89        # Write the completed line to the output image
90        newDataset.GetRasterBand(1).WriteRaster(0, line, red_band.XSize, 1, \
91                                         outputLine, buf_xsize=red_band.XSize, 
92                                         buf_ysize=1, buf_type=gdal.GDT_Float32)
93   
94        # Delete the output line following write
95        del outputLine
96    print >> sys.stderr,'NDVI Calculated and Outputted to File'
97    print >> sys.stderr,dir(newDataset)
98    newDataset.FlushCache()
99    vsiFile=gdal.VSIFOpenL("/vsimem//output"+conf["lenv"]["sid"],"r")
100    i=0
101    while gdal.VSIFSeekL(vsiFile,0,os.SEEK_END)>0:
102        i+=1
103    fileSize=gdal.VSIFTellL(vsiFile)
104    gdal.VSIFSeekL(vsiFile,0,os.SEEK_SET)
105    outputs["raster"]["value"]=gdal.VSIFReadL(fileSize,1,vsiFile)
106    gdal.Unlink("/vsimem//output"+conf["lenv"]["sid"])
107    return 3
Note: See TracBrowser for help on using the repository browser.

Search

Context Navigation

ZOO Sponsors

http://www.zoo-project.org/trac/chrome/site/img/geolabs-logo.pnghttp://www.zoo-project.org/trac/chrome/site/img/neogeo-logo.png http://www.zoo-project.org/trac/chrome/site/img/apptech-logo.png http://www.zoo-project.org/trac/chrome/site/img/3liz-logo.png http://www.zoo-project.org/trac/chrome/site/img/gateway-logo.png

Become a sponsor !

Knowledge partners

http://www.zoo-project.org/trac/chrome/site/img/ocu-logo.png http://www.zoo-project.org/trac/chrome/site/img/gucas-logo.png http://www.zoo-project.org/trac/chrome/site/img/polimi-logo.png http://www.zoo-project.org/trac/chrome/site/img/fem-logo.png http://www.zoo-project.org/trac/chrome/site/img/supsi-logo.png http://www.zoo-project.org/trac/chrome/site/img/cumtb-logo.png

Become a knowledge partner

Related links

http://zoo-project.org/img/ogclogo.png http://zoo-project.org/img/osgeologo.png