rt_band_get_summary_stats()

rt_bandstats rt_band_get_summary_stats	(	rt_band	band,
		int	exclude_nodata_value,
		double	sample,
		int	inc_vals,
		uint64_t *	cK,
		double *	cM,
		double *	cQ
	)

Compute summary statistics for a band.

Parameters

band	: the band to query for summary stats
exclude_nodata_value	: if non-zero, ignore nodata values
sample	: percentage of pixels to sample
inc_vals	: flag to include values in return struct
cK	: number of pixels counted thus far in coverage
cM	: M component of 1-pass stddev for coverage
cQ	: Q component of 1-pass stddev for coverage

Returns

the summary statistics for a band or NULL

Definition at line 110 of file rt_statistics.c.

  {
        uint32_t x = 0;
        int64_t y = 0;
        uint32_t z = 0;
        uint32_t offset = 0;
        uint32_t diff = 0;
        int rtn;
        int hasnodata = FALSE;
        double nodata = 0;
        double *values = NULL;
        double value;
        int isnodata = 0;
        rt_bandstats stats = NULL;
 
        uint32_t do_sample = 0;
        uint32_t sample_size = 0;
        uint32_t sample_per = 0;
        uint32_t sample_int = 0;
        uint32_t i = 0;
        uint32_t j = 0;
        double sum = 0;
        uint32_t k = 0;
        double M = 0;
        double Q = 0;
 
#if POSTGIS_DEBUG_LEVEL > 0
        clock_t start, stop;
        double elapsed = 0;
#endif
 
        RASTER_DEBUG(3, "starting");
#if POSTGIS_DEBUG_LEVEL > 0
        start = clock();
#endif
 
        assert(NULL != band);
 
        /* band is empty (width < 1 || height < 1) */
        if (band->width < 1 || band->height < 1) {
                stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
                if (NULL == stats) {
                        rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
                        return NULL;
                }
 
                rtwarn("Band is empty as width and/or height is 0");
 
                stats->sample = 1;
                stats->sorted = 0;
                stats->values = NULL;
                stats->count = 0;
                stats->min = stats->max = 0;
                stats->sum = 0;
                stats->mean = 0;
                stats->stddev = -1;
 
                return stats;
        }
 
        hasnodata = rt_band_get_hasnodata_flag(band);
        if (hasnodata != FALSE)
                rt_band_get_nodata(band, &nodata);
        else
                exclude_nodata_value = 0;
 
        RASTER_DEBUGF(3, "nodata = %f", nodata);
        RASTER_DEBUGF(3, "hasnodata = %d", hasnodata);
        RASTER_DEBUGF(3, "exclude_nodata_value = %d", exclude_nodata_value);
 
        /* entire band is nodata */
        if (rt_band_get_isnodata_flag(band) != FALSE) {
                stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
                if (NULL == stats) {
                        rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
                        return NULL;
                }
 
                stats->sample = 1;
                stats->sorted = 0;
                stats->values = NULL;
 
                if (exclude_nodata_value) {
                        rtwarn("All pixels of band have the NODATA value");
 
                        stats->count = 0;
                        stats->min = stats->max = 0;
                        stats->sum = 0;
                        stats->mean = 0;
                        stats->stddev = -1;
                }
                else {
                        stats->count = band->width * band->height;
                        stats->min = stats->max = nodata;
                        stats->sum = stats->count * nodata;
                        stats->mean = nodata;
                        stats->stddev = 0;
                }
 
                return stats;
        }
 
        /* clamp percentage */
        if (
                (sample < 0 || FLT_EQ(sample, 0.0)) ||
                (sample > 1 || FLT_EQ(sample, 1.0))
        ) {
                do_sample = 0;
                sample = 1;
        }
        else
                do_sample = 1;
        RASTER_DEBUGF(3, "do_sample = %d", do_sample);
 
        /* sample all pixels */
        if (!do_sample) {
                sample_size = band->width * band->height;
                sample_per = band->height;
        }
        /*
         randomly sample a percentage of available pixels
         sampling method is known as
                "systematic random sample without replacement"
        */
        else {
                sample_size = round((band->width * band->height) * sample);
                sample_per = round(sample_size / band->width);
                if (sample_per < 1)
                        sample_per = 1;
                sample_int = round(band->height / sample_per);
                srand(time(NULL));
        }
 
        RASTER_DEBUGF(3, "sampling %d of %d available pixels w/ %d per set"
                , sample_size, (band->width * band->height), sample_per);
 
        if (inc_vals) {
                values = rtalloc(sizeof(double) * sample_size);
                if (NULL == values) {
                        rtwarn("Could not allocate memory for values");
                        inc_vals = 0;
                }
        }
 
        /* initialize stats */
        stats = (rt_bandstats) rtalloc(sizeof(struct rt_bandstats_t));
        if (NULL == stats) {
                rterror("rt_band_get_summary_stats: Could not allocate memory for stats");
                return NULL;
        }
        stats->sample = sample;
        stats->count = 0;
        stats->sum = 0;
        stats->mean = 0;
        stats->stddev = -1;
        stats->min = stats->max = 0;
        stats->values = NULL;
        stats->sorted = 0;
 
        for (x = 0, j = 0, k = 0; x < band->width; x++) {
                y = -1;
                diff = 0;
 
                for (i = 0, z = 0; i < sample_per; i++) {
                        if (!do_sample)
                                y = i;
                        else {
                                offset = (rand() % sample_int) + 1;
                                y += diff + offset;
                                diff = sample_int - offset;
                        }
                        RASTER_DEBUGF(5, "(x, y, z) = (%u, %lld, %u)", x, y, z);
                        if (y >= band->height || z > sample_per) break;
 
                        rtn = rt_band_get_pixel(band, x, y, &value, &isnodata);
 
                        j++;
                        if (rtn == ES_NONE && (!exclude_nodata_value || (exclude_nodata_value && !isnodata))) {
 
                                /* inc_vals set, collect pixel values */
                                if (inc_vals) values[k] = value;
 
                                /* average */
                                k++;
                                sum += value;
 
                                /*
                                        one-pass standard deviation
                                        http://www.eecs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
                                */
                                if (k == 1) {
                                        Q = 0;
                                        M = value;
                                }
                                else {
                                        Q += (((k  - 1) * pow(value - M, 2)) / k);
                                        M += ((value - M ) / k);
                                }
 
                                /* coverage one-pass standard deviation */
                                if (NULL != cK) {
                                        (*cK)++;
                                        if (*cK == 1) {
                                                *cQ = 0;
                                                *cM = value;
                                        }
                                        else {
                                                *cQ += (((*cK  - 1) * pow(value - *cM, 2)) / *cK);
                                                *cM += ((value - *cM ) / *cK);
                                        }
                                }
 
                                /* min/max */
                                if (stats->count < 1) {
                                        stats->count = 1;
                                        stats->min = stats->max = value;
                                }
                                else {
                                        if (value < stats->min)
                                                stats->min = value;
                                        if (value > stats->max)
                                                stats->max = value;
                                }
 
                        }
 
                        z++;
                }
        }
 
        RASTER_DEBUG(3, "sampling complete");
 
        stats->count = k;
        if (k > 0) {
                if (inc_vals) {
                        /* free unused memory */
                        if (sample_size != k) {
                                values = rtrealloc(values, k * sizeof(double));
                        }
 
                        stats->values = values;
                }
 
                stats->sum = sum;
                stats->mean = sum / k;
 
                /* standard deviation */
                if (!do_sample)
                        stats->stddev = sqrt(Q / k);
                /* sample deviation */
                else {
                        if (k < 2)
                                stats->stddev = -1;
                        else
                                stats->stddev = sqrt(Q / (k - 1));
                }
        }
        /* inc_vals thus values allocated but not used */
        else if (inc_vals)
                rtdealloc(values);
 
        /* if do_sample is one */
        if (do_sample && k < 1)
                rtwarn("All sampled pixels of band have the NODATA value");
 
#if POSTGIS_DEBUG_LEVEL > 0
        stop = clock();
        elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
        RASTER_DEBUGF(3, "(time, count, mean, stddev, min, max) = (%0.4f, %d, %f, %f, %f, %f)",
                elapsed, stats->count, stats->mean, stats->stddev, stats->min, stats->max);
#endif
 
        RASTER_DEBUG(3, "done");
        return stats;
}

References rt_bandstats_t::count, ES_NONE, FALSE, FLT_EQ, rt_bandstats_t::max, rt_bandstats_t::mean, rt_bandstats_t::min, RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_band_get_nodata(), rt_band_get_pixel(), rtalloc(), rtdealloc(), rterror(), rtrealloc(), rtwarn(), rt_bandstats_t::sample, rt_bandstats_t::sorted, rt_bandstats_t::stddev, rt_bandstats_t::sum, and rt_bandstats_t::values.

Referenced by RASTER_colorMap(), RASTER_histogram(), RASTER_quantile(), RASTER_summaryStats(), RASTER_summaryStats_transfn(), RASTER_summaryStatsCoverage(), test_band_stats(), and test_gdal_warp_preserves_data().

Here is the call graph for this function:

Here is the caller graph for this function: