1

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:

STATISTICS PYTHON-CODE:

MAXIMUM = 0.9896427989006
MEAN = -0.21678630311319
MINIMUM = -0.69955635070801
STDDEV = 0.24564237169808

STATISTICS QGIS:

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.

closed as off-topic by PolyGeo Sep 14 '16 at 7:52

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Questions seeking help to debug/write/improve code must include the desired behavior, a specific problem or error and the shortest code necessary to reproduce it in the question itself. Providing a clear problem statement and evidence of a code attempt will help others to help you. See: How to create a Minimal, Complete, and Verifiable example." – PolyGeo
If this question can be reworded to fit the rules in the help center, please edit the question.

  • You didn't try filtering the Landsat image in QGIS? (i.e., B2 > 0, B5 > 0) ? How are the statistics after that. Are there negative values in the layers? – reima Jun 1 '15 at 10:35
  • I tried, but it doesn't make a difference – Johan Jun 1 '15 at 11:45
  • To be sure, you should test it on something where you know the result. Construct a small raster with known values where you can hand-compute the result you want. Then test. For all you know, your QGIS and python calculations could both be wrong! – Spacedman Jun 1 '15 at 12:06

Browse other questions tagged or ask your own question.