PostGIS  3.4.0dev-r@@SVN_REVISION@@

◆ RASTER_histogram()

Datum RASTER_histogram ( PG_FUNCTION_ARGS  )

Definition at line 919 of file rtpg_statistics.c.

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

References ovdump::band, rt_bandstats_t::count, genraster::count, FALSE, FLT_EQ, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rtrowdump::raster, result, 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: