I am working on a script that automatically calculates something for me. To be sure, I checked the result of the calculation with QGIS. Here are the statistics:
MAXIMUM = 0.9896427989006 MEAN = -0.21678630311319 MINIMUM = -0.69955635070801 STDDEV = 0.24564237169808
MAXIMUM = 1 MEAN = -0.4204168536468 MINIMUM = -0.69955635070801 STDDEV = 0.17747956958404
I created my script based on this: https://gist.github.com/jgomezdans/1336359. Here it is:
def filter2 (band2, band5): ds2 = gdal.Open(band2) data2 = ds2.ReadAsArray() ds5 = gdal.Open(band5) data5 = ds5.ReadAsArray() check = np.logical_and(data2 > 0, data5 > 0) filter2 = np.where(check, (data2 - data5) / (data2 + data5), 0) return filter2
In QGIS I used the raster calculator:
(band2 - band5) / (band2 + band5)
Bands 2 and 5 originate from a landsat 7 image and are converted to reflectance values (so no DN).
I think I am doing something wrong with the check and np.where lines but when I just calculate the ratio in the script as:
filter2 = (data2 - data5) / (data2 + data5)
the answer is NaN.
I am not sure if I should worry about the differences in the statistics, but I think I should. And also, I'd like to know why I get NaN when I do not use the check.