PostGIS  2.1.10dev-r@@SVN_REVISION@@
Datum RASTER_summaryStatsCoverage ( PG_FUNCTION_ARGS  )

Definition at line 8101 of file rt_pg.c.

References ovdump::band, rt_bandstats_t::count, FALSE, FLT_EQ, rt_bandstats_t::max, rt_bandstats_t::mean, rt_bandstats_t::min, rtrowdump::raster, result, rt_band_destroy(), rt_band_get_summary_stats(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), rt_bandstats_t::sample, rt_bandstats_t::sorted, rtgdalraster::sql, rt_bandstats_t::stddev, rt_bandstats_t::sum, TRUE, and rt_bandstats_t::values.

8102 {
8103  text *tablenametext = NULL;
8104  char *tablename = NULL;
8105  text *colnametext = NULL;
8106  char *colname = NULL;
8107  int32_t bandindex = 1;
8108  bool exclude_nodata_value = TRUE;
8109  double sample = 0;
8110 
8111  int len = 0;
8112  char *sql = NULL;
8113  int spi_result;
8114  Portal portal;
8115  TupleDesc tupdesc;
8116  SPITupleTable *tuptable = NULL;
8117  HeapTuple tuple;
8118  Datum datum;
8119  bool isNull = FALSE;
8120 
8121  rt_pgraster *pgraster = NULL;
8122  rt_raster raster = NULL;
8123  rt_band band = NULL;
8124  int num_bands = 0;
8125  uint64_t cK = 0;
8126  double cM = 0;
8127  double cQ = 0;
8128  rt_bandstats stats = NULL;
8129  rt_bandstats rtn = NULL;
8130 
8131  int values_length = 6;
8132  Datum values[values_length];
8133  bool nulls[values_length];
8134  Datum result;
8135 
8136  /* tablename is null, return null */
8137  if (PG_ARGISNULL(0)) {
8138  elog(NOTICE, "Table name must be provided");
8139  PG_RETURN_NULL();
8140  }
8141  tablenametext = PG_GETARG_TEXT_P(0);
8142  tablename = text_to_cstring(tablenametext);
8143  if (!strlen(tablename)) {
8144  elog(NOTICE, "Table name must be provided");
8145  PG_RETURN_NULL();
8146  }
8147 
8148  /* column name is null, return null */
8149  if (PG_ARGISNULL(1)) {
8150  elog(NOTICE, "Column name must be provided");
8151  PG_RETURN_NULL();
8152  }
8153  colnametext = PG_GETARG_TEXT_P(1);
8154  colname = text_to_cstring(colnametext);
8155  if (!strlen(colname)) {
8156  elog(NOTICE, "Column name must be provided");
8157  PG_RETURN_NULL();
8158  }
8159 
8160  /* band index is 1-based */
8161  if (!PG_ARGISNULL(2))
8162  bandindex = PG_GETARG_INT32(2);
8163 
8164  /* exclude_nodata_value flag */
8165  if (!PG_ARGISNULL(3))
8166  exclude_nodata_value = PG_GETARG_BOOL(3);
8167 
8168  /* sample % */
8169  if (!PG_ARGISNULL(4)) {
8170  sample = PG_GETARG_FLOAT8(4);
8171  if (sample < 0 || sample > 1) {
8172  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
8173  rt_raster_destroy(raster);
8174  PG_RETURN_NULL();
8175  }
8176  else if (FLT_EQ(sample, 0.0))
8177  sample = 1;
8178  }
8179  else
8180  sample = 1;
8181 
8182  /* iterate through rasters of coverage */
8183  /* connect to database */
8184  spi_result = SPI_connect();
8185  if (spi_result != SPI_OK_CONNECT) {
8186  pfree(sql);
8187  elog(ERROR, "RASTER_summaryStatsCoverage: Could not connect to database using SPI");
8188  PG_RETURN_NULL();
8189  }
8190 
8191  /* create sql */
8192  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
8193  sql = (char *) palloc(len);
8194  if (NULL == sql) {
8195  if (SPI_tuptable) SPI_freetuptable(tuptable);
8196  SPI_finish();
8197  elog(ERROR, "RASTER_summaryStatsCoverage: Could not allocate memory for sql");
8198  PG_RETURN_NULL();
8199  }
8200 
8201  /* get cursor */
8202  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
8203  portal = SPI_cursor_open_with_args(
8204  "coverage",
8205  sql,
8206  0, NULL,
8207  NULL, NULL,
8208  TRUE, 0
8209  );
8210  pfree(sql);
8211 
8212  /* process resultset */
8213  SPI_cursor_fetch(portal, TRUE, 1);
8214  while (SPI_processed == 1 && SPI_tuptable != NULL) {
8215  tupdesc = SPI_tuptable->tupdesc;
8216  tuptable = SPI_tuptable;
8217  tuple = tuptable->vals[0];
8218 
8219  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
8220  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
8221 
8222  if (SPI_tuptable) SPI_freetuptable(tuptable);
8223  SPI_cursor_close(portal);
8224  SPI_finish();
8225 
8226  if (NULL != rtn) pfree(rtn);
8227  elog(ERROR, "RASTER_summaryStatsCoverage: Could not get raster of coverage");
8228  PG_RETURN_NULL();
8229  }
8230  else if (isNull) {
8231  SPI_cursor_fetch(portal, TRUE, 1);
8232  continue;
8233  }
8234 
8235  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
8236 
8237  raster = rt_raster_deserialize(pgraster, FALSE);
8238  if (!raster) {
8239 
8240  if (SPI_tuptable) SPI_freetuptable(tuptable);
8241  SPI_cursor_close(portal);
8242  SPI_finish();
8243 
8244  if (NULL != rtn) pfree(rtn);
8245  elog(ERROR, "RASTER_summaryStatsCoverage: Could not deserialize raster");
8246  PG_RETURN_NULL();
8247  }
8248 
8249  /* inspect number of bands */
8250  num_bands = rt_raster_get_num_bands(raster);
8251  if (bandindex < 1 || bandindex > num_bands) {
8252  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
8253 
8254  rt_raster_destroy(raster);
8255 
8256  if (SPI_tuptable) SPI_freetuptable(tuptable);
8257  SPI_cursor_close(portal);
8258  SPI_finish();
8259 
8260  if (NULL != rtn) pfree(rtn);
8261  PG_RETURN_NULL();
8262  }
8263 
8264  /* get band */
8265  band = rt_raster_get_band(raster, bandindex - 1);
8266  if (!band) {
8267  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
8268 
8269  rt_raster_destroy(raster);
8270 
8271  if (SPI_tuptable) SPI_freetuptable(tuptable);
8272  SPI_cursor_close(portal);
8273  SPI_finish();
8274 
8275  if (NULL != rtn) pfree(rtn);
8276  PG_RETURN_NULL();
8277  }
8278 
8279  /* we don't need the raw values, hence the zero parameter */
8280  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 0, &cK, &cM, &cQ);
8281 
8282  rt_band_destroy(band);
8283  rt_raster_destroy(raster);
8284 
8285  if (NULL == stats) {
8286  elog(NOTICE, "Could not compute summary statistics for band at index %d. Returning NULL", bandindex);
8287 
8288  if (SPI_tuptable) SPI_freetuptable(tuptable);
8289  SPI_cursor_close(portal);
8290  SPI_finish();
8291 
8292  if (NULL != rtn) pfree(rtn);
8293  PG_RETURN_NULL();
8294  }
8295 
8296  /* initialize rtn */
8297  if (stats->count > 0) {
8298  if (NULL == rtn) {
8299  rtn = (rt_bandstats) SPI_palloc(sizeof(struct rt_bandstats_t));
8300  if (NULL == rtn) {
8301 
8302  if (SPI_tuptable) SPI_freetuptable(tuptable);
8303  SPI_cursor_close(portal);
8304  SPI_finish();
8305 
8306  elog(ERROR, "RASTER_summaryStatsCoverage: Could not allocate memory for summary stats of coverage");
8307  PG_RETURN_NULL();
8308  }
8309 
8310  rtn->sample = stats->sample;
8311  rtn->count = stats->count;
8312  rtn->min = stats->min;
8313  rtn->max = stats->max;
8314  rtn->sum = stats->sum;
8315  rtn->mean = stats->mean;
8316  rtn->stddev = -1;
8317 
8318  rtn->values = NULL;
8319  rtn->sorted = 0;
8320  }
8321  else {
8322  rtn->count += stats->count;
8323  rtn->sum += stats->sum;
8324 
8325  if (stats->min < rtn->min)
8326  rtn->min = stats->min;
8327  if (stats->max > rtn->max)
8328  rtn->max = stats->max;
8329  }
8330  }
8331 
8332  pfree(stats);
8333 
8334  /* next record */
8335  SPI_cursor_fetch(portal, TRUE, 1);
8336  }
8337 
8338  if (SPI_tuptable) SPI_freetuptable(tuptable);
8339  SPI_cursor_close(portal);
8340  SPI_finish();
8341 
8342  if (NULL == rtn) {
8343  elog(ERROR, "RASTER_summaryStatsCoverage: Could not compute coverage summary stats");
8344  PG_RETURN_NULL();
8345  }
8346 
8347  /* coverage mean and deviation */
8348  rtn->mean = rtn->sum / rtn->count;
8349  /* sample deviation */
8350  if (rtn->sample > 0 && rtn->sample < 1)
8351  rtn->stddev = sqrt(cQ / (rtn->count - 1));
8352  /* standard deviation */
8353  else
8354  rtn->stddev = sqrt(cQ / rtn->count);
8355 
8356  /* Build a tuple descriptor for our result type */
8357  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
8358  ereport(ERROR, (
8359  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
8360  errmsg(
8361  "function returning record called in context "
8362  "that cannot accept type record"
8363  )
8364  ));
8365  }
8366 
8367  BlessTupleDesc(tupdesc);
8368 
8369  memset(nulls, FALSE, sizeof(bool) * values_length);
8370 
8371  values[0] = Int64GetDatum(rtn->count);
8372  if (rtn->count > 0) {
8373  values[1] = Float8GetDatum(rtn->sum);
8374  values[2] = Float8GetDatum(rtn->mean);
8375  values[3] = Float8GetDatum(rtn->stddev);
8376  values[4] = Float8GetDatum(rtn->min);
8377  values[5] = Float8GetDatum(rtn->max);
8378  }
8379  else {
8380  nulls[1] = TRUE;
8381  nulls[2] = TRUE;
8382  nulls[3] = TRUE;
8383  nulls[4] = TRUE;
8384  nulls[5] = TRUE;
8385  }
8386 
8387  /* build a tuple */
8388  tuple = heap_form_tuple(tupdesc, values, nulls);
8389 
8390  /* make the tuple into a datum */
8391  result = HeapTupleGetDatum(tuple);
8392 
8393  /* clean up */
8394  pfree(rtn);
8395 
8396  PG_RETURN_DATUM(result);
8397 }
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_api.c:5677
double mean
Definition: rt_api.h:2282
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
uint32_t count
Definition: rt_api.h:2277
struct rt_bandstats_t * rt_bandstats
Definition: rt_api.h:137
double sum
Definition: rt_api.h:2281
tuple band
Definition: ovdump.py:57
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:123
char ** result
Definition: liblwgeom.h:218
double min
Definition: rt_api.h:2279
double * values
Definition: rt_api.h:2285
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
double sample
Definition: rt_api.h:2276
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_api.c:1650
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.
Definition: rt_api.c:3269
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
#define FALSE
Definition: dbfopen.c:169
Struct definitions.
Definition: rt_api.h:2175
double max
Definition: rt_api.h:2280
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_api.c:8350
double stddev
Definition: rt_api.h:2283
#define TRUE
Definition: dbfopen.c:170

Here is the call graph for this function: