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

Last change on this file since 426 was 426, checked in by djay, 11 years ago

Add header to ndvi service.

  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-python
File size: 5.2 KB
Line 
1#
2# Author : Gérald FENOY
3#
4# Copyright 2009-2013 GeoLabs SARL. All rights reserved.
5#
6# Permission is hereby granted, free of charge, to any person obtaining a
7# copy of this software and associated documentation files (the
8# "Software"), to deal in the Software without restriction, including with
9# out limitation the rights to use, copy, modify, merge, publish,
10# distribute, sublicense, and/or sell copies of the Software, and to
11# permit persons to whom the Software is furnished to do so, subject to
12# the following conditions:
13#
14# The above copyright notice and this permission notice shall be included
15# in all copies or substantial portions of the Software.
16#
17# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
18# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
19# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
20# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
21# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
22# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
23# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
24#
25
26
27import sys, os, struct
28
29import osgeo.gdal as gdal
30
31# Calculate and output NDVI from raster bands
32def ExtractNDVI(conf, inputs, outputs):
33   
34    # Open the input dataset
35    gdal.FileFromMemBuffer('/vsimem//temp.tif', inputs["raster"]["value"])
36    dataset = gdal.Open( '/vsimem//temp.tif')
37    if dataset is None:
38        conf["lenv"]["message"]="The dataset could not be openned properly"
39        return 4
40
41    # Create the output dataset
42    driver = gdal.GetDriverByName( "GTiff" )
43    # Get the spatial information from the input file
44    geoTransform=None
45    geoProjection=None
46    print >> sys.stderr,dir(dataset) 
47    try:
48        geoTransform = dataset.GetGeoTransform()
49    except:
50        print >> sys.stderr, "Unable to load geotransform"
51    try:
52        geoProjection = dataset.GetProjection()
53    except:
54        print >> sys.stderr, "Unable to load projection"
55
56    # Create an output file of the same size as the inputted
57    # image but with only 1 output image band.
58    newDataset = driver.Create("/vsimem//output"+conf["lenv"]["sid"], \
59                               dataset.RasterXSize, dataset.RasterYSize,1, \
60                               gdal.GDT_Float32)
61    # Set spatial informations of the new image.
62    if geoTransform:
63        newDataset.SetGeoTransform(geoTransform)
64    if geoProjection:
65        newDataset.SetProjection(geoProjection)
66    if newDataset is None:
67        conf["lenv"]["message"]='Could not create output image'
68        return 4
69   
70    # Get the RED and NIR image bands of the image
71    red_id=int(inputs["red"]["value"])
72    nir_id=int(inputs["nir"]["value"])
73    red_band = dataset.GetRasterBand(red_id) # RED BAND
74    nir_band = dataset.GetRasterBand(nir_id) # NIR BAND
75
76    # Loop through each line in turn.
77    numLines = red_band.YSize
78    for line in range(numLines):
79        # Define variable for output line.
80        outputLine = ''
81        # Read in data for the current line from the
82        # image band representing the red wavelength
83        red_scanline = red_band.ReadRaster( 0, line, red_band.XSize, 1, \
84                                          red_band.XSize, 1, gdal.GDT_Float32 )
85        # Unpack the line of data to be read as floating point data
86        red_tuple = struct.unpack('f' * red_band.XSize, red_scanline)
87           
88        # Read in data for the current line from the
89        # image band representing the NIR wavelength
90        nir_scanline = nir_band.ReadRaster( 0, line, nir_band.XSize, 1, \
91                                          nir_band.XSize, 1, gdal.GDT_Float32 )
92        # Unpack the line of data to be read as floating point data
93        nir_tuple = struct.unpack('f' * nir_band.XSize, nir_scanline)
94
95        # Loop through the columns within the image
96        for i in range(len(red_tuple)):
97            # Calculate the NDVI for the current pixel.
98            ndvi_lower = (nir_tuple[i] + red_tuple[i])
99            ndvi_upper = (nir_tuple[i] - red_tuple[i])
100            ndvi = 0
101            # Becareful of zero divide
102            if ndvi_lower == 0:
103                ndvi = 0
104            else:
105                ndvi = ndvi_upper/ndvi_lower
106            # Add the current pixel to the output line
107            outputLine = outputLine + struct.pack('f', ndvi)
108        # Write the completed line to the output image
109        newDataset.GetRasterBand(1).WriteRaster(0, line, red_band.XSize, 1, \
110                                         outputLine, buf_xsize=red_band.XSize, 
111                                         buf_ysize=1, buf_type=gdal.GDT_Float32)
112   
113        # Delete the output line following write
114        del outputLine
115    print >> sys.stderr,'NDVI Calculated and Outputted to File'
116    print >> sys.stderr,dir(newDataset)
117    newDataset.FlushCache()
118    vsiFile=gdal.VSIFOpenL("/vsimem//output"+conf["lenv"]["sid"],"r")
119    i=0
120    while gdal.VSIFSeekL(vsiFile,0,os.SEEK_END)>0:
121        i+=1
122    fileSize=gdal.VSIFTellL(vsiFile)
123    gdal.VSIFSeekL(vsiFile,0,os.SEEK_SET)
124    outputs["raster"]["value"]=gdal.VSIFReadL(fileSize,1,vsiFile)
125    gdal.Unlink("/vsimem//output"+conf["lenv"]["sid"])
126    return 3
Note: See TracBrowser for help on using the repository browser.

Search

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