Verifying zonal statistics output

[chrisvwn]

Opening this issue to comprehensively address an issue from the Ishara Data blog here raised by @bidur on the blog here. Copied for convenience and to allow gathering all info together.

[chrisvwn]

Copied from https://isharadata.blogspot.com/2017/09/rnightlights-satellite-nightlight-data.html?showComment=1521461093580#c1389136912456374501. Edited to the relevant parts only.

=========================

Bidur DevkotaMarch 19, 2018 at 3:04 PM

Hi Chris, I have some more queries:

1. I manually calculated radiance statistics like mean, max, min, average for VIIRS 201204 tif image for Macau . Then I compared the same statistics by using RNightlights library:

lowStats <- getCtryNlData(ctryCode = "MAC", nlPeriods = nlRange("201204", "201204"), nlType = "VIIRS", nlStats =c("min","max","sum","mean","median"), ignoreMissing=FALSE)

I used the same tif file and same adm shapefile file of Macau for both cases. I got same values for max and min in both cases. However, I got slightly different values for mean and sum for overall Macau. The results for both are in this file https://www.dropbox.com/s/1i7wctl4qe6z4o3/macau2014.xls?dl=0

...

Bidur DevkotaMarch 19, 2018 at 7:21 PM

Code used to manually calculated radiance statistics https://www.dropbox.com/s/68u67ymt1uz03rc/ntlStatistics4Macau.r?dl=0

Hi Bidur,

1. From your code only one issue jumps out i.e. the treatment of negative values.

values(ntlCountry)[values(ntlCountry)<0] = 0 # Replace Raster Values Less than 0 to 0 ( Thresholding values by ZERO)

Since these values are unknown I recommend setting them to NA as well since zero is a valid measurement and using it will affect most aggregate calculations. Setting NA should exclude them from all calculations when you use na.rm=TRUE.

I believe was you wanted to use the threshold for high values? I would do this separately from negative values since they are on opposite ends of your valid data so set it to work on some greater-than-zero value. However, while this should resolve for calculations that are affected by NAs, this should not affect sums so there may be something else going on there. I will run your code later tonight and get back to you .

...

Hi Bidur,

Apologies for taking a while to get back. I have run your code and have seen te differences in results. So I see 2 things happening here. The first is a limitation of the aggregation of aggregated values. Rnightlights currently calculates sub-polygon aggregates at the lowest country level. Getting the mean of this at the country level will not be the same as calculating the mean directly at the country level. The same is true for other aggregate stats except the sum (which is why the example only does the sum). This is addressed in the next version of Rnightlights, out in the next week or so, where one can specify the admin level at which to perform aggregations.

The second issue is more elusive. Using your calculations I end up with a vector of 180 pixels while Rnightlights gives 168. I think something is happening here as regards pixel inclusion - maybe different strategies are being employed. I am still trying to figure out if this is a difference due to the raster::extract function or some other part of code. Will get back to you on this.

[chrisvwn]

A screenshot of the excel file mentioned above for convenience:

[chrisvwn]

Confirming issue with Rnightlights is the calculation of agg functions such as mean at the lowest level and aggregating to higher levels. As mentioned above the mean of means is not equivalent. For confirmation attaching screenshot of zonal statistics from QGIS. There is a difference in pixel counts - Rnightlights yields 168 pixels while QGIS yields 167. However, note the totals are the same. This will be investigated farther.