4

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.

  • What do you mean by "binary rasters"? What is their datatype? Float32? Byte (8-bit)? – Logan Byers Aug 8 '17 at 14:32
  • This means just that they only have two values, which is NA and 254. Sorry for the confusion. – bellackn Aug 8 '17 at 21:22
2

If you actually have NaN values in your rasters, the following should work (untested).

gdal_calc.py -A raster_a.tif -B raster_b.tif --outfile=output.tif \
--calc="(isnan(A) * (B==254)) + ((A==254) * isnan(B)) + nan*((isnan(A) * isnan(B)) + ((A==254) * (B==254)))" \
--NoDataValue=0

The idea is to get 1 where one raster is NaN and the other is 254 and NaN where both rasters are the same value. The operation where the rasters are the same can't be performed using equal or == because NaN doesn't evaluate that way (check np.nan == np.nan in a python console, it returns False).

EDIT: Knowing now that the rasters only contain 254 and 0 it should be as easy as doing an inequality check like you had done in QGIS. You have probably lost your image footprint mask if you had one prior to making raster_a/b

gdal_calc.py -A raster_a.tif -B raster_b.tif --outfile=output.tif \
--calc="(A != B)" \
--NoDataValue=0 --type=Byte
| improve this answer | |
  • Unfortunately, this is still returning a raster with NaNs only. – bellackn Aug 8 '17 at 21:33
  • What do you mean by "if you actually have NaN values", might these be sth. else? If I use Identify Features in QGIS and click on the input rasters where there should be NaN, it throws "no data". In the layer's transparency tab, "No data value: 0" is checked. If I uncheck this, Identify Features shows 0 at the same spots. – bellackn Aug 8 '17 at 21:40
  • @bellackn This means you don't actually have NaN values in your raster. What you think are NaN are actually 0 (interpreted as NoData because of a metadata entry). In this case, you can use the second set of calculations in my post, which is just the same inequality check you knew you need to use. – Logan Byers Aug 8 '17 at 22:31
  • This does not work, it results in a raster with NaNs only again. As you already mentioned, this is the same I already tried in the very beginning - because this would be intuitive (though it was theoretically wrong, as I know now ;)). – bellackn Aug 8 '17 at 22:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.