QGIS Zonal Statistics and No_data

Question

I think QGIS Zonal Statistics algorithm does not handle No_data values correctly.

The expected behaviour is that if I set some raster cells to No_data (for example with the method suggested here, using gdal_translate), then Zonal Statistics will ignore those values.

But in fact it treats them as -3.4e+38, so all the statistics are wrong. See for example Antartica in this pictures, all 0 values in layer 1 have been set to No_data in the Converted layer, so Mean and Sum stats should be 0. But they are not.

I tried with SAGA's "raster statistics for polygons", and there is the same problem.

Is it worth filing a bug report? Or am I missing something?

Just for reference, layer 1 is N balance on landscape, taken from earthstat.

EDIT

I tested it all again after @user30184's answer. The problem actually arises if I apply a mask to the raster layer first, to change all values <= 0 to 0.

Using the raster calculator:

(layer1 > 0) * layer1

Then I convert 0 values to no_data with the method described above, and finally I run Zonal Statistics.

You can see the results in the picture below. Statistics with "_" as prefix are calculated without the mask layer step. Statistics with "t" as prefix are calculated with the mask layer. As you can see for Antartica, I get funny values in the "t" results.

Answer 1

The description of the -a_nodata parameter from the documentation of gdal_translate https://gdal.org/programs/gdal_translate.html

-a_nodata [value]

Assign a specified nodata value to output bands. It can be set to none to avoid setting a nodata value to the output file if one exists for the source file. Note that, if the input dataset has a nodata value, this does not cause pixel values that are equal to that nodata value to be changed to the value specified with this option.

I believe that in your original image value -3.4e+38 is used to mean nodata. Then you fired a gdal_translate command

gdal_translate -a_nodata 0.0 ...

What happens is that all pixel values of 0.0 has been marked as nodata in the target image. If -3.4e+38 used to mean nodata before it does not mean it anymore but all those pixels contain now real data for QGIS and SAGA. The -a_nodata parameter does not change the pixel values, it is just touching the metadata.

Without having test data it is hard to know but it feels like you should have -3.4e+38 as nodata value. If it is not so in the original image (check it with gdalinfo) you can do set if by running

gdal_translate -a_nodata -3.4e+38  ...

EDIT

I downloaded the image from your link and I am correcting myself. The original image has no nodata set.

gdalinfo NitrogenBalanceOnLandscape_140Crops.tif
…
Band 1 Block=432x224 Type=Float32, ColorInterp=Gray
  Overviews: 2160x1080, 1080x540, 540x270, 270x135

I set the nodata to 0 with the similar command that you used and now I have

gdalinfo with_nodata.tif
…
Band 1 Block=4320x1 Type=Float32, ColorInterp=Gray
  NoData Value=0

I digitized two polygons with QGIS, one over the ocean and one partly over Africa. Then I used the zonal statistics tool from the processing toolbox. The results in the table show zero and null for the all nodata area, and reasonable values for the African triangle. For me it seems that nodata pixels are correctly excluded from the zonal analysis.

EDIT 2

By following the workflow that is described in the question I am getting the same results. For some reason that I can't explain the workflow original image -> QGIS raster calculator (layer1 > 0) * layer1 -> assign nodata=0 with gdal_tranlate results an image that has 4320 pixels with value -3.4028234663853e+38

  Min=-340282346638529993179660072199368212480.000 Max=2752376.250
  Minimum=-340282346638529993179660072199368212480.000, Maximum=2752376.250, Mean=-2090745287322099897642461654474031104.000, StdDev=26590834825873999547824535346249465856.000
  256 buckets from -3.4095e+38 to 6.6722e+35:
  4320 0 0 0 0 0 0 0 0 0 0 0 0 0

I can get correct result with gdal_calc https://gdal.org/programs/gdal_calc.html but even here I can't explain why.

The command is

gdal_calc -A nitro.tif --outfile nitro_calc.tif --calc="(A>0)*A" --NoDataValue=0

Compare the resuls from the original workflow and the gdal_calc workflow (columns _ccount, _csum, and _cmean)

Understanding what happens requires probably very deep knowledge on how QGIS, GDAL, and numpy work. Perhaps you can reach best experts from qgis-users or qgis-developers mailing list.

EDIT 3

The origin of the issue is probably in the QGIS raster calculator. The lowest pixel row of the original image contains just 0 values and raster calculator with (layer1 > 0) * layer1 should not do anything for the lowest row. But actually it considers that the lowest row is nodata and changes the values which used to be 0 into -3.4e+38. Later GDAL is told to label 0 as nodata and then -3.4e+38 turns into real data and pixels are used in the zonal statistics.