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 | |
---|
27 | import sys, os, struct |
---|
28 | |
---|
29 | import osgeo.gdal as gdal |
---|
30 | |
---|
31 | # Calculate and output NDVI from raster bands |
---|
32 | def 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 |
---|