PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ RASTER_histogram()

Datum RASTER_histogram ( PG_FUNCTION_ARGS  )

Definition at line 917 of file rtpg_statistics.c.

918 {
919  FuncCallContext *funcctx;
920  TupleDesc tupdesc;
921 
922  int i;
923  rt_histogram hist;
924  rt_histogram hist2;
925  int call_cntr;
926  int max_calls;
927 
928  /* first call of function */
929  if (SRF_IS_FIRSTCALL()) {
930  MemoryContext oldcontext;
931 
932  rt_pgraster *pgraster = NULL;
933  rt_raster raster = NULL;
934  rt_band band = NULL;
935  int32_t bandindex = 1;
936  int num_bands = 0;
937  bool exclude_nodata_value = TRUE;
938  double sample = 0;
939  uint32_t bin_count = 0;
940  double *bin_width = NULL;
941  uint32_t bin_width_count = 0;
942  double width = 0;
943  bool right = FALSE;
944  double min = 0;
945  double max = 0;
946  rt_bandstats stats = NULL;
947  uint32_t count;
948 
949  int j;
950  int n;
951 
952  ArrayType *array;
953  Oid etype;
954  Datum *e;
955  bool *nulls;
956  int16 typlen;
957  bool typbyval;
958  char typalign;
959 
960  POSTGIS_RT_DEBUG(3, "RASTER_histogram: Starting");
961 
962  /* create a function context for cross-call persistence */
963  funcctx = SRF_FIRSTCALL_INIT();
964 
965  /* switch to memory context appropriate for multiple function calls */
966  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
967 
968  /* pgraster is null, return nothing */
969  if (PG_ARGISNULL(0)) {
970  MemoryContextSwitchTo(oldcontext);
971  SRF_RETURN_DONE(funcctx);
972  }
973  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
974 
975  raster = rt_raster_deserialize(pgraster, FALSE);
976  if (!raster) {
977  PG_FREE_IF_COPY(pgraster, 0);
978  MemoryContextSwitchTo(oldcontext);
979  elog(ERROR, "RASTER_histogram: Cannot deserialize raster");
980  SRF_RETURN_DONE(funcctx);
981  }
982 
983  /* band index is 1-based */
984  if (!PG_ARGISNULL(1))
985  bandindex = PG_GETARG_INT32(1);
986  num_bands = rt_raster_get_num_bands(raster);
987  if (bandindex < 1 || bandindex > num_bands) {
988  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
989  rt_raster_destroy(raster);
990  PG_FREE_IF_COPY(pgraster, 0);
991  MemoryContextSwitchTo(oldcontext);
992  SRF_RETURN_DONE(funcctx);
993  }
994 
995  /* exclude_nodata_value flag */
996  if (!PG_ARGISNULL(2))
997  exclude_nodata_value = PG_GETARG_BOOL(2);
998 
999  /* sample % */
1000  if (!PG_ARGISNULL(3)) {
1001  sample = PG_GETARG_FLOAT8(3);
1002  if (sample < 0 || sample > 1) {
1003  elog(NOTICE, "Invalid sample percentage (must be between 0 and 1). Returning NULL");
1004  rt_raster_destroy(raster);
1005  PG_FREE_IF_COPY(pgraster, 0);
1006  MemoryContextSwitchTo(oldcontext);
1007  SRF_RETURN_DONE(funcctx);
1008  }
1009  else if (FLT_EQ(sample, 0.0))
1010  sample = 1;
1011  }
1012  else
1013  sample = 1;
1014 
1015  /* bin_count */
1016  if (!PG_ARGISNULL(4)) {
1017  bin_count = PG_GETARG_INT32(4);
1018  if (bin_count < 1) bin_count = 0;
1019  }
1020 
1021  /* bin_width */
1022  if (!PG_ARGISNULL(5)) {
1023  array = PG_GETARG_ARRAYTYPE_P(5);
1024  etype = ARR_ELEMTYPE(array);
1025  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1026 
1027  switch (etype) {
1028  case FLOAT4OID:
1029  case FLOAT8OID:
1030  break;
1031  default:
1032  rt_raster_destroy(raster);
1033  PG_FREE_IF_COPY(pgraster, 0);
1034  MemoryContextSwitchTo(oldcontext);
1035  elog(ERROR, "RASTER_histogram: Invalid data type for width");
1036  SRF_RETURN_DONE(funcctx);
1037  break;
1038  }
1039 
1040  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
1041  &nulls, &n);
1042 
1043  bin_width = palloc(sizeof(double) * n);
1044  for (i = 0, j = 0; i < n; i++) {
1045  if (nulls[i]) continue;
1046 
1047  switch (etype) {
1048  case FLOAT4OID:
1049  width = (double) DatumGetFloat4(e[i]);
1050  break;
1051  case FLOAT8OID:
1052  width = (double) DatumGetFloat8(e[i]);
1053  break;
1054  }
1055 
1056  if (width < 0 || FLT_EQ(width, 0.0)) {
1057  elog(NOTICE, "Invalid value for width (must be greater than 0). Returning NULL");
1058  pfree(bin_width);
1059  rt_raster_destroy(raster);
1060  PG_FREE_IF_COPY(pgraster, 0);
1061  MemoryContextSwitchTo(oldcontext);
1062  SRF_RETURN_DONE(funcctx);
1063  }
1064 
1065  bin_width[j] = width;
1066  POSTGIS_RT_DEBUGF(5, "bin_width[%d] = %f", j, bin_width[j]);
1067  j++;
1068  }
1069  bin_width_count = j;
1070 
1071  if (j < 1) {
1072  pfree(bin_width);
1073  bin_width = NULL;
1074  }
1075  }
1076 
1077  /* right */
1078  if (!PG_ARGISNULL(6))
1079  right = PG_GETARG_BOOL(6);
1080 
1081  /* min */
1082  if (!PG_ARGISNULL(7)) min = PG_GETARG_FLOAT8(7);
1083 
1084  /* max */
1085  if (!PG_ARGISNULL(8)) max = PG_GETARG_FLOAT8(8);
1086 
1087  /* get band */
1088  band = rt_raster_get_band(raster, bandindex - 1);
1089  if (!band) {
1090  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
1091  rt_raster_destroy(raster);
1092  PG_FREE_IF_COPY(pgraster, 0);
1093  MemoryContextSwitchTo(oldcontext);
1094  SRF_RETURN_DONE(funcctx);
1095  }
1096 
1097  /* get stats */
1098  stats = rt_band_get_summary_stats(band, (int) exclude_nodata_value, sample, 1, NULL, NULL, NULL);
1099  rt_band_destroy(band);
1100  rt_raster_destroy(raster);
1101  PG_FREE_IF_COPY(pgraster, 0);
1102  if (NULL == stats || NULL == stats->values) {
1103  elog(NOTICE, "Cannot compute summary statistics for band at index %d", bandindex);
1104  MemoryContextSwitchTo(oldcontext);
1105  SRF_RETURN_DONE(funcctx);
1106  }
1107  else if (stats->count < 1) {
1108  elog(NOTICE, "Cannot compute histogram for band at index %d as the band has no values", bandindex);
1109  MemoryContextSwitchTo(oldcontext);
1110  SRF_RETURN_DONE(funcctx);
1111  }
1112 
1113  /* get histogram */
1114  hist = rt_band_get_histogram(stats, (uint32_t)bin_count, bin_width, bin_width_count, right, min, max, &count);
1115  if (bin_width_count) pfree(bin_width);
1116  pfree(stats);
1117  if (NULL == hist || !count) {
1118  elog(NOTICE, "Cannot compute histogram for band at index %d", bandindex);
1119  MemoryContextSwitchTo(oldcontext);
1120  SRF_RETURN_DONE(funcctx);
1121  }
1122 
1123  POSTGIS_RT_DEBUGF(3, "%d bins returned", count);
1124 
1125  /* Store needed information */
1126  funcctx->user_fctx = hist;
1127 
1128  /* total number of tuples to be returned */
1129  funcctx->max_calls = count;
1130 
1131  /* Build a tuple descriptor for our result type */
1132  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
1133  ereport(ERROR, (
1134  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
1135  errmsg(
1136  "function returning record called in context "
1137  "that cannot accept type record"
1138  )
1139  ));
1140  }
1141 
1142  BlessTupleDesc(tupdesc);
1143  funcctx->tuple_desc = tupdesc;
1144 
1145  MemoryContextSwitchTo(oldcontext);
1146  }
1147 
1148  /* stuff done on every call of the function */
1149  funcctx = SRF_PERCALL_SETUP();
1150 
1151  call_cntr = funcctx->call_cntr;
1152  max_calls = funcctx->max_calls;
1153  tupdesc = funcctx->tuple_desc;
1154  hist2 = funcctx->user_fctx;
1155 
1156  /* do when there is more left to send */
1157  if (call_cntr < max_calls) {
1158  Datum values[VALUES_LENGTH];
1159  bool nulls[VALUES_LENGTH];
1160  HeapTuple tuple;
1161  Datum result;
1162 
1163  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
1164 
1165  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
1166 
1167  values[0] = Float8GetDatum(hist2[call_cntr].min);
1168  values[1] = Float8GetDatum(hist2[call_cntr].max);
1169  values[2] = Int64GetDatum(hist2[call_cntr].count);
1170  values[3] = Float8GetDatum(hist2[call_cntr].percent);
1171 
1172  /* build a tuple */
1173  tuple = heap_form_tuple(tupdesc, values, nulls);
1174 
1175  /* make the tuple into a datum */
1176  result = HeapTupleGetDatum(tuple);
1177 
1178  SRF_RETURN_NEXT(funcctx, result);
1179  }
1180  /* do when there is no more left */
1181  else {
1182  pfree(hist2);
1183  SRF_RETURN_DONE(funcctx);
1184  }
1185 }
TRUE
#define TRUE
Definition: dbfopen.c:169
FALSE
#define FALSE
Definition: dbfopen.c:168
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2234
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:340
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
rt_band_get_histogram
rt_histogram rt_band_get_histogram(rt_bandstats stats, uint32_t bin_count, double *bin_widths, uint32_t bin_widths_count, int right, double min, double max, uint32_t *rtn_count)
Count the distribution of data.
Definition: rt_statistics.c:412
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.
Definition: rt_statistics.c:110
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:725
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
genraster.count
int count
Definition: genraster.py:56
ovdump.band
band
Definition: ovdump.py:57
rtrowdump.raster
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
VALUES_LENGTH
#define VALUES_LENGTH
Definition: rtpg_statistics.c:2396
POSTGIS_RT_DEBUG
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h:61
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:65
rt_band_t
Definition: librtcore.h:2313
rt_bandstats_t::count
uint32_t count
Definition: librtcore.h:2360
rt_bandstats_t::values
double * values
Definition: librtcore.h:2368
rt_bandstats_t
Definition: librtcore.h:2358
rt_histogram_t
Definition: librtcore.h:2373
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2250
rt_raster_t
Definition: librtcore.h:2284
uint32_t
unsigned int uint32_t
Definition: uthash.h:78

References ovdump::band, rt_bandstats_t::count, genraster::count, FALSE, FLT_EQ, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rtrowdump::raster, rt_band_destroy(), rt_band_get_histogram(), rt_band_get_summary_stats(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), TRUE, rt_bandstats_t::values, and VALUES_LENGTH.

Here is the call graph for this function: