Saturday, 23 February 2019

qgis - gdal_calc - How to calculate where raster A is not equal to raster B?


I have two binary rasters with values NaN and 254. I want to generate a new raster that consists of the pixels that have value NaN in raster A and 254 in raster B and vice versa - those should be 1, the rest should be NaN again.


I want to use gdal_calc.py in a bash script for that. My code snippet looks like this:


gdal_calc.py -A S2_2015_AOI_NDBI_bin.tif -B S2_2015_AOI_NDSI_bin.tif --outfile=S2_2015_AOI_diff_auto.tif --calc="1*(A!=B)" --NoDataValue=0


But all it returns is a new raster with NaN everywhere.


I also tried --calc="1*(not_equal(A,B), --calc="1*(logical_or(logical_and(A==254, not_equal(B, 254)), logical_and(B==254, not_equal(A, 254))))", and several other combinations I could come up with. Nothing worked yet.


Strange thing is, it works flawlessly in the QGIS Raster Calculator with the formula ("S2_2015_AOI_NDBI_bin@1" != "S2_2015_AOI_NDSI_bin@1") = 1, but I want this to be automated. Also, this calculator throws 0 instead of NaN, which is not suitable for further processing as I plan it.


The two input rasters were previously generated by:


gdal_calc.py -A S2_2015_AOI_NDBI.tif --outfile=S2_2015_AOI_NDBI_bin.tif --calc="254*(A>-0.15)" --NoDataValue=0

and


gdal_calc.py -A S2_2015_AOI_NDSI.tif --outfile=S2_2015_AOI_NDSI_bin.tif --calc="254*(A>0.4)" --NoDataValue=0`

The rasters all have the same CRS. I'm running GDAL v. 1.11.3, it's a bit outdated, might this be the problem?





UPDATE


I solved the problem by skipping the step of creating the "binary" input rasters with 254 and NaN values (see question's second and third code blocks). The code looks like this now:


gdal_calc.py -A S2_2015_AOI_NDBI.tif -B S2_2015_AOI_NDSI.tif --outfile=out.tif --calc="1*(logical_and(A>-0.15, B<0.4))" --NoDataValue=0

This delivers the result I wanted.


While this is a suitable workaround for me personally, it still does not answer the original question.




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...