1

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.

Convert Raster

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.

enter image description here

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.

Results with and without mask step

  • Provide some test data (links are OK) and step-by-test procedure to reproduce your issue. – user30184 Feb 4 at 12:58
  • edited the question, thank you for your comments, is it clear now? – ffgg Feb 4 at 15:30
2

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.

enter image description here

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)

enter image description here

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.

enter image description here

| improve this answer | |
  • Just as context to @user30184's answer, QGIS in my experience handles nodata quite well. But sometimes nodata inconsistencies arise when combining (native) QGIS, GDAL, SAGA, Grass, etc. It can therefore be helpful to use the command line or other paramaters to make sure it is being passed correctly. I would also check in the layer's QGIS layer properties that the nodata setting is set correctly. – Houska Feb 4 at 16:00
  • Thanks for the answer, I've edited the question accordingly, I have narrowed down my problem a bit further, could you please have a look at it again? – ffgg Feb 10 at 14:04
  • Many thanks, accepted because it gives a workaround the problem: using GDAL's raster calculator tool instead of QGIS standard one for applying the mask. It works as expected with this solution. I might bring this up in the developers mailing list as suggested, depending on other users' comments and feedback – ffgg Feb 11 at 15:08

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.