PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ 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
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}
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition cu_print.c:267
#define TRUE
Definition dbfopen.c:73
#define FALSE
Definition dbfopen.c:72
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition rt_raster.c:86
#define FLT_EQ(x, y)
Definition librtcore.h:2436
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition rt_band.c:499
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:376
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.
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.
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition rt_raster.c:385
int count
Definition genraster.py:57
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition rtrowdump.py:125
#define VALUES_LENGTH
#define POSTGIS_RT_DEBUG(level, msg)
Definition rtpostgis.h:65
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition rtpostgis.h:69
uint32_t count
Definition librtcore.h:2562
double * values
Definition librtcore.h:2570
Struct definitions.
Definition librtcore.h:2452

References rt_bandstats_t::count, FALSE, FLT_EQ, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, 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: