Index: trunk/zoo-project/zoo-services/gdal/ndvi/cgi-env/ExtractNDVI.zcfg
===================================================================
--- trunk/zoo-project/zoo-services/gdal/ndvi/cgi-env/ExtractNDVI.zcfg (revision 334)
+++ trunk/zoo-project/zoo-services/gdal/ndvi/cgi-env/ExtractNDVI.zcfg (revision 334)
@@ -0,0 +1,57 @@
+[ExtractNDVI]
+ Title = Normalized Difference Vegetation Index.
+ Abstract = Create raster file containing the Normalized Difference Vegetation Index.
+ Profile = urn:ogc:wps:1.0.0:buffer
+ processVersion = 2
+ storeSupported = true
+ statusSupported = true
+ serviceProvider = ndvi
+ serviceType = Python
+
+ title = Demo
+
+
+ [raster]
+ Title = The raster file to be used
+ Abstract = URI to a raster file that contains the mentionned two bands.
+ minOccurs = 1
+ maxOccurs = 1
+
+
+ mimeType = image/tiff
+
+
+ [red]
+ Title = The red band number
+ Abstract = The RED band number to use.
+ minOccurs = 0
+ maxOccurs = 1
+
+ DataType = int
+
+ value = 1
+
+
+ [nir]
+ Title = The NIR band number
+ Abstract = The NIR band number to use.
+ minOccurs = 0
+ maxOccurs = 1
+
+ DataType = int
+
+ value = 2
+
+
+
+
+ [raster]
+ Title = The resulting raster
+ Abstract = A single band raster containing the NDVI of the original raster using formula : (NIR - RED) / (NIR + RED) .
+
+
+ mimeType = image/tiff
+ extension = tif
+
+
+
Index: trunk/zoo-project/zoo-services/gdal/ndvi/cgi-env/ndvi.py
===================================================================
--- trunk/zoo-project/zoo-services/gdal/ndvi/cgi-env/ndvi.py (revision 334)
+++ trunk/zoo-project/zoo-services/gdal/ndvi/cgi-env/ndvi.py (revision 334)
@@ -0,0 +1,102 @@
+#! /usr/bin/env python
+
+import sys, os, struct
+
+import osgeo.gdal as gdal
+
+# Calculate and output NDVI from raster bands
+def ExtractNDVI(conf, inputs, outputs):
+
+ # Open the input dataset
+ gdal.FileFromMemBuffer('/vsimem//temp.tif', inputs["raster"]["value"])
+ dataset = gdal.Open( '/vsimem//temp.tif')
+ if dataset is None:
+ conf["lenv"]["message"]="The dataset could not be openned properly"
+ return 4
+
+ # Create the output dataset
+ driver = gdal.GetDriverByName( "GTiff" )
+ # Get the spatial information from the input file
+ geoTransform=None
+ geoProjection=None
+ print >> sys.stderr,dir(dataset)
+ try:
+ geoTransform = dataset.GetGeoTransform()
+ except:
+ print >> sys.stderr, "Unable to load geotransform"
+ try:
+ geoProjection = dataset.GetProjection()
+ except:
+ print >> sys.stderr, "Unable to load projection"
+
+ # Create an output file of the same size as the inputted
+ # image but with only 1 output image band.
+ newDataset = driver.Create("/vsimem//output"+conf["lenv"]["sid"], \
+ dataset.RasterXSize, dataset.RasterYSize,1, \
+ gdal.GDT_Float32)
+ # Set spatial informations of the new image.
+ if geoTransform:
+ newDataset.SetGeoTransform(geoTransform)
+ if geoProjection:
+ newDataset.SetProjection(geoProjection)
+ if newDataset is None:
+ conf["lenv"]["message"]='Could not create output image'
+ return 4
+
+ # Get the RED and NIR image bands of the image
+ red_id=int(inputs["red"]["value"])
+ nir_id=int(inputs["nir"]["value"])
+ red_band = dataset.GetRasterBand(red_id) # RED BAND
+ nir_band = dataset.GetRasterBand(nir_id) # NIR BAND
+
+ # Loop through each line in turn.
+ numLines = red_band.YSize
+ for line in range(numLines):
+ # Define variable for output line.
+ outputLine = ''
+ # Read in data for the current line from the
+ # image band representing the red wavelength
+ red_scanline = red_band.ReadRaster( 0, line, red_band.XSize, 1, \
+ red_band.XSize, 1, gdal.GDT_Float32 )
+ # Unpack the line of data to be read as floating point data
+ red_tuple = struct.unpack('f' * red_band.XSize, red_scanline)
+
+ # Read in data for the current line from the
+ # image band representing the NIR wavelength
+ nir_scanline = nir_band.ReadRaster( 0, line, nir_band.XSize, 1, \
+ nir_band.XSize, 1, gdal.GDT_Float32 )
+ # Unpack the line of data to be read as floating point data
+ nir_tuple = struct.unpack('f' * nir_band.XSize, nir_scanline)
+
+ # Loop through the columns within the image
+ for i in range(len(red_tuple)):
+ # Calculate the NDVI for the current pixel.
+ ndvi_lower = (nir_tuple[i] + red_tuple[i])
+ ndvi_upper = (nir_tuple[i] - red_tuple[i])
+ ndvi = 0
+ # Becareful of zero divide
+ if ndvi_lower == 0:
+ ndvi = 0
+ else:
+ ndvi = ndvi_upper/ndvi_lower
+ # Add the current pixel to the output line
+ outputLine = outputLine + struct.pack('f', ndvi)
+ # Write the completed line to the output image
+ newDataset.GetRasterBand(1).WriteRaster(0, line, red_band.XSize, 1, \
+ outputLine, buf_xsize=red_band.XSize,
+ buf_ysize=1, buf_type=gdal.GDT_Float32)
+
+ # Delete the output line following write
+ del outputLine
+ print >> sys.stderr,'NDVI Calculated and Outputted to File'
+ print >> sys.stderr,dir(newDataset)
+ newDataset.FlushCache()
+ vsiFile=gdal.VSIFOpenL("/vsimem//output"+conf["lenv"]["sid"],"r")
+ i=0
+ while gdal.VSIFSeekL(vsiFile,0,os.SEEK_END)>0:
+ i+=1
+ fileSize=gdal.VSIFTellL(vsiFile)
+ gdal.VSIFSeekL(vsiFile,0,os.SEEK_SET)
+ outputs["raster"]["value"]=gdal.VSIFReadL(fileSize,1,vsiFile)
+ gdal.Unlink("/vsimem//output"+conf["lenv"]["sid"])
+ return 3