PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ RASTER_valueCountCoverage()

Datum RASTER_valueCountCoverage ( PG_FUNCTION_ARGS  )

Definition at line 2620 of file rtpg_statistics.c.

References ovdump::band, genraster::count, rt_valuecount_t::count, rtpg_summarystats_arg_t::exclude_nodata_value, FALSE, FLT_EQ, free(), rt_valuecount_t::percent, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rtrowdump::raster, rt_band_destroy(), rt_band_get_value_count(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), rtgdalraster::sql, TRUE, genraster::value, rt_valuecount_t::value, and VALUES_LENGTH.

Referenced by RASTER_valueCount().

2620  {
2621  FuncCallContext *funcctx;
2622  TupleDesc tupdesc;
2623 
2624  int i;
2625  uint64_t covcount = 0;
2626  uint64_t covtotal = 0;
2627  rt_valuecount covvcnts = NULL;
2628  rt_valuecount covvcnts2;
2629  int call_cntr;
2630  int max_calls;
2631 
2632  POSTGIS_RT_DEBUG(3, "RASTER_valueCountCoverage: Starting");
2633 
2634  /* first call of function */
2635  if (SRF_IS_FIRSTCALL()) {
2636  MemoryContext oldcontext;
2637 
2638  text *tablenametext = NULL;
2639  char *tablename = NULL;
2640  text *colnametext = NULL;
2641  char *colname = NULL;
2642  int32_t bandindex = 1;
2643  bool exclude_nodata_value = TRUE;
2644  double *search_values = NULL;
2645  uint32_t search_values_count = 0;
2646  double roundto = 0;
2647 
2648  int len = 0;
2649  char *sql = NULL;
2650  int spi_result;
2651  Portal portal;
2652  SPITupleTable *tuptable = NULL;
2653  HeapTuple tuple;
2654  Datum datum;
2655  bool isNull = FALSE;
2656  rt_pgraster *pgraster = NULL;
2657  rt_raster raster = NULL;
2658  rt_band band = NULL;
2659  int num_bands = 0;
2660  uint32_t count;
2661  uint32_t total;
2662  rt_valuecount vcnts = NULL;
2663  int exists = 0;
2664 
2665  int j;
2666  int n;
2667 
2668  ArrayType *array;
2669  Oid etype;
2670  Datum *e;
2671  bool *nulls;
2672  int16 typlen;
2673  bool typbyval;
2674  char typalign;
2675 
2676  /* create a function context for cross-call persistence */
2677  funcctx = SRF_FIRSTCALL_INIT();
2678 
2679  /* switch to memory context appropriate for multiple function calls */
2680  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
2681 
2682  /* tablename is null, return null */
2683  if (PG_ARGISNULL(0)) {
2684  elog(NOTICE, "Table name must be provided");
2685  MemoryContextSwitchTo(oldcontext);
2686  SRF_RETURN_DONE(funcctx);
2687  }
2688  tablenametext = PG_GETARG_TEXT_P(0);
2689  tablename = text_to_cstring(tablenametext);
2690  if (!strlen(tablename)) {
2691  elog(NOTICE, "Table name must be provided");
2692  MemoryContextSwitchTo(oldcontext);
2693  SRF_RETURN_DONE(funcctx);
2694  }
2695  POSTGIS_RT_DEBUGF(3, "tablename = %s", tablename);
2696 
2697  /* column name is null, return null */
2698  if (PG_ARGISNULL(1)) {
2699  elog(NOTICE, "Column name must be provided");
2700  MemoryContextSwitchTo(oldcontext);
2701  SRF_RETURN_DONE(funcctx);
2702  }
2703  colnametext = PG_GETARG_TEXT_P(1);
2704  colname = text_to_cstring(colnametext);
2705  if (!strlen(colname)) {
2706  elog(NOTICE, "Column name must be provided");
2707  MemoryContextSwitchTo(oldcontext);
2708  SRF_RETURN_DONE(funcctx);
2709  }
2710  POSTGIS_RT_DEBUGF(3, "colname = %s", colname);
2711 
2712  /* band index is 1-based */
2713  if (!PG_ARGISNULL(2))
2714  bandindex = PG_GETARG_INT32(2);
2715 
2716  /* exclude_nodata_value flag */
2717  if (!PG_ARGISNULL(3))
2718  exclude_nodata_value = PG_GETARG_BOOL(3);
2719 
2720  /* search values */
2721  if (!PG_ARGISNULL(4)) {
2722  array = PG_GETARG_ARRAYTYPE_P(4);
2723  etype = ARR_ELEMTYPE(array);
2724  get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
2725 
2726  switch (etype) {
2727  case FLOAT4OID:
2728  case FLOAT8OID:
2729  break;
2730  default:
2731  MemoryContextSwitchTo(oldcontext);
2732  elog(ERROR, "RASTER_valueCountCoverage: Invalid data type for values");
2733  SRF_RETURN_DONE(funcctx);
2734  break;
2735  }
2736 
2737  deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
2738  &nulls, &n);
2739 
2740  search_values = palloc(sizeof(double) * n);
2741  for (i = 0, j = 0; i < n; i++) {
2742  if (nulls[i]) continue;
2743 
2744  switch (etype) {
2745  case FLOAT4OID:
2746  search_values[j] = (double) DatumGetFloat4(e[i]);
2747  break;
2748  case FLOAT8OID:
2749  search_values[j] = (double) DatumGetFloat8(e[i]);
2750  break;
2751  }
2752 
2753  POSTGIS_RT_DEBUGF(5, "search_values[%d] = %f", j, search_values[j]);
2754  j++;
2755  }
2756  search_values_count = j;
2757 
2758  if (j < 1) {
2759  pfree(search_values);
2760  search_values = NULL;
2761  }
2762  }
2763 
2764  /* roundto */
2765  if (!PG_ARGISNULL(5)) {
2766  roundto = PG_GETARG_FLOAT8(5);
2767  if (roundto < 0.) roundto = 0;
2768  }
2769 
2770  /* iterate through rasters of coverage */
2771  /* connect to database */
2772  spi_result = SPI_connect();
2773  if (spi_result != SPI_OK_CONNECT) {
2774 
2775  if (search_values_count) pfree(search_values);
2776 
2777  MemoryContextSwitchTo(oldcontext);
2778  elog(ERROR, "RASTER_valueCountCoverage: Cannot connect to database using SPI");
2779  SRF_RETURN_DONE(funcctx);
2780  }
2781 
2782  /* create sql */
2783  len = sizeof(char) * (strlen("SELECT \"\" FROM \"\" WHERE \"\" IS NOT NULL") + (strlen(colname) * 2) + strlen(tablename) + 1);
2784  sql = (char *) palloc(len);
2785  if (NULL == sql) {
2786 
2787  if (SPI_tuptable) SPI_freetuptable(tuptable);
2788  SPI_finish();
2789 
2790  if (search_values_count) pfree(search_values);
2791 
2792  MemoryContextSwitchTo(oldcontext);
2793  elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for sql");
2794  SRF_RETURN_DONE(funcctx);
2795  }
2796 
2797  /* get cursor */
2798  snprintf(sql, len, "SELECT \"%s\" FROM \"%s\" WHERE \"%s\" IS NOT NULL", colname, tablename, colname);
2799  POSTGIS_RT_DEBUGF(3, "RASTER_valueCountCoverage: %s", sql);
2800  portal = SPI_cursor_open_with_args(
2801  "coverage",
2802  sql,
2803  0, NULL,
2804  NULL, NULL,
2805  TRUE, 0
2806  );
2807  pfree(sql);
2808 
2809  /* process resultset */
2810  SPI_cursor_fetch(portal, TRUE, 1);
2811  while (SPI_processed == 1 && SPI_tuptable != NULL) {
2812  tupdesc = SPI_tuptable->tupdesc;
2813  tuptable = SPI_tuptable;
2814  tuple = tuptable->vals[0];
2815 
2816  datum = SPI_getbinval(tuple, tupdesc, 1, &isNull);
2817  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
2818 
2819  if (SPI_tuptable) SPI_freetuptable(tuptable);
2820  SPI_cursor_close(portal);
2821  SPI_finish();
2822 
2823  if (NULL != covvcnts) pfree(covvcnts);
2824  if (search_values_count) pfree(search_values);
2825 
2826  MemoryContextSwitchTo(oldcontext);
2827  elog(ERROR, "RASTER_valueCountCoverage: Cannot get raster of coverage");
2828  SRF_RETURN_DONE(funcctx);
2829  }
2830  else if (isNull) {
2831  SPI_cursor_fetch(portal, TRUE, 1);
2832  continue;
2833  }
2834 
2835  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(datum);
2836 
2837  raster = rt_raster_deserialize(pgraster, FALSE);
2838  if (!raster) {
2839 
2840  if (SPI_tuptable) SPI_freetuptable(tuptable);
2841  SPI_cursor_close(portal);
2842  SPI_finish();
2843 
2844  if (NULL != covvcnts) pfree(covvcnts);
2845  if (search_values_count) pfree(search_values);
2846 
2847  MemoryContextSwitchTo(oldcontext);
2848  elog(ERROR, "RASTER_valueCountCoverage: Cannot deserialize raster");
2849  SRF_RETURN_DONE(funcctx);
2850  }
2851 
2852  /* inspect number of bands*/
2853  num_bands = rt_raster_get_num_bands(raster);
2854  if (bandindex < 1 || bandindex > num_bands) {
2855  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2856 
2857  rt_raster_destroy(raster);
2858 
2859  if (SPI_tuptable) SPI_freetuptable(tuptable);
2860  SPI_cursor_close(portal);
2861  SPI_finish();
2862 
2863  if (NULL != covvcnts) pfree(covvcnts);
2864  if (search_values_count) pfree(search_values);
2865 
2866  MemoryContextSwitchTo(oldcontext);
2867  SRF_RETURN_DONE(funcctx);
2868  }
2869 
2870  /* get band */
2871  band = rt_raster_get_band(raster, bandindex - 1);
2872  if (!band) {
2873  elog(NOTICE, "Cannot find band at index %d. Returning NULL", bandindex);
2874 
2875  rt_raster_destroy(raster);
2876 
2877  if (SPI_tuptable) SPI_freetuptable(tuptable);
2878  SPI_cursor_close(portal);
2879  SPI_finish();
2880 
2881  if (NULL != covvcnts) pfree(covvcnts);
2882  if (search_values_count) pfree(search_values);
2883 
2884  MemoryContextSwitchTo(oldcontext);
2885  SRF_RETURN_DONE(funcctx);
2886  }
2887 
2888  /* get counts of values */
2889  vcnts = rt_band_get_value_count(band, (int) exclude_nodata_value, search_values, search_values_count, roundto, &total, &count);
2890  rt_band_destroy(band);
2891  rt_raster_destroy(raster);
2892  if (NULL == vcnts || !count) {
2893  elog(NOTICE, "Cannot count the values for band at index %d", bandindex);
2894 
2895  if (SPI_tuptable) SPI_freetuptable(tuptable);
2896  SPI_cursor_close(portal);
2897  SPI_finish();
2898 
2899  if (NULL != covvcnts) free(covvcnts);
2900  if (search_values_count) pfree(search_values);
2901 
2902  MemoryContextSwitchTo(oldcontext);
2903  SRF_RETURN_DONE(funcctx);
2904  }
2905 
2906  POSTGIS_RT_DEBUGF(3, "%d value counts returned", count);
2907 
2908  if (NULL == covvcnts) {
2909  covvcnts = (rt_valuecount) SPI_palloc(sizeof(struct rt_valuecount_t) * count);
2910  if (NULL == covvcnts) {
2911 
2912  if (SPI_tuptable) SPI_freetuptable(tuptable);
2913  SPI_cursor_close(portal);
2914  SPI_finish();
2915 
2916  if (search_values_count) pfree(search_values);
2917 
2918  MemoryContextSwitchTo(oldcontext);
2919  elog(ERROR, "RASTER_valueCountCoverage: Cannot allocate memory for value counts of coverage");
2920  SRF_RETURN_DONE(funcctx);
2921  }
2922 
2923  for (i = 0; i < count; i++) {
2924  covvcnts[i].value = vcnts[i].value;
2925  covvcnts[i].count = vcnts[i].count;
2926  covvcnts[i].percent = -1;
2927  }
2928 
2929  covcount = count;
2930  }
2931  else {
2932  for (i = 0; i < count; i++) {
2933  exists = 0;
2934 
2935  for (j = 0; j < covcount; j++) {
2936  if (FLT_EQ(vcnts[i].value, covvcnts[j].value)) {
2937  exists = 1;
2938  break;
2939  }
2940  }
2941 
2942  if (exists) {
2943  covvcnts[j].count += vcnts[i].count;
2944  }
2945  else {
2946  covcount++;
2947  covvcnts = SPI_repalloc(covvcnts, sizeof(struct rt_valuecount_t) * covcount);
2948  if (NULL == covvcnts) {
2949 
2950  if (SPI_tuptable) SPI_freetuptable(tuptable);
2951  SPI_cursor_close(portal);
2952  SPI_finish();
2953 
2954  if (search_values_count) pfree(search_values);
2955  if (NULL != covvcnts) free(covvcnts);
2956 
2957  MemoryContextSwitchTo(oldcontext);
2958  elog(ERROR, "RASTER_valueCountCoverage: Cannot change allocated memory for value counts of coverage");
2959  SRF_RETURN_DONE(funcctx);
2960  }
2961 
2962  covvcnts[covcount - 1].value = vcnts[i].value;
2963  covvcnts[covcount - 1].count = vcnts[i].count;
2964  covvcnts[covcount - 1].percent = -1;
2965  }
2966  }
2967  }
2968 
2969  covtotal += total;
2970 
2971  pfree(vcnts);
2972 
2973  /* next record */
2974  SPI_cursor_fetch(portal, TRUE, 1);
2975  }
2976 
2977  if (SPI_tuptable) SPI_freetuptable(tuptable);
2978  SPI_cursor_close(portal);
2979  SPI_finish();
2980 
2981  if (search_values_count) pfree(search_values);
2982 
2983  /* compute percentages */
2984  for (i = 0; i < covcount; i++) {
2985  covvcnts[i].percent = (double) covvcnts[i].count / covtotal;
2986  }
2987 
2988  /* Store needed information */
2989  funcctx->user_fctx = covvcnts;
2990 
2991  /* total number of tuples to be returned */
2992  funcctx->max_calls = covcount;
2993 
2994  /* Build a tuple descriptor for our result type */
2995  if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
2996  ereport(ERROR, (
2997  errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2998  errmsg(
2999  "function returning record called in context "
3000  "that cannot accept type record"
3001  )
3002  ));
3003  }
3004 
3005  BlessTupleDesc(tupdesc);
3006  funcctx->tuple_desc = tupdesc;
3007 
3008  MemoryContextSwitchTo(oldcontext);
3009  }
3010 
3011  /* stuff done on every call of the function */
3012  funcctx = SRF_PERCALL_SETUP();
3013 
3014  call_cntr = funcctx->call_cntr;
3015  max_calls = funcctx->max_calls;
3016  tupdesc = funcctx->tuple_desc;
3017  covvcnts2 = funcctx->user_fctx;
3018 
3019  /* do when there is more left to send */
3020  if (call_cntr < max_calls) {
3021  Datum values[VALUES_LENGTH];
3022  bool nulls[VALUES_LENGTH];
3023  HeapTuple tuple;
3024  Datum result;
3025 
3026  POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
3027 
3028  memset(nulls, FALSE, sizeof(bool) * VALUES_LENGTH);
3029 
3030  values[0] = Float8GetDatum(covvcnts2[call_cntr].value);
3031  values[1] = UInt32GetDatum(covvcnts2[call_cntr].count);
3032  values[2] = Float8GetDatum(covvcnts2[call_cntr].percent);
3033 
3034  /* build a tuple */
3035  tuple = heap_form_tuple(tupdesc, values, nulls);
3036 
3037  /* make the tuple into a datum */
3038  result = HeapTupleGetDatum(tuple);
3039 
3040  SRF_RETURN_NEXT(funcctx, result);
3041  }
3042  /* do when there is no more left */
3043  else {
3044  pfree(covvcnts2);
3045  SRF_RETURN_DONE(funcctx);
3046  }
3047 }
VALUES_LENGTH
#define VALUES_LENGTH
Definition: rtpg_statistics.c:2406
rt_band_t
Definition: librtcore.h:2264
rt_raster_get_num_bands
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
rtrowdump.raster
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
rt_valuecount_t::count
uint32_t count
Definition: librtcore.h:2376
ovdump.band
band
Definition: ovdump.py:57
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2185
rt_band_destroy
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:242
rt_valuecount_t::value
double value
Definition: librtcore.h:2375
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:65
uint32_t
unsigned int uint32_t
Definition: uthash.h:78
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
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
rt_valuecount_t::percent
double percent
Definition: librtcore.h:2377
rt_valuecount
struct rt_valuecount_t * rt_valuecount
Definition: librtcore.h:153
rt_band_get_value_count
rt_valuecount rt_band_get_value_count(rt_band band, int exclude_nodata_value, double *search_values, uint32_t search_values_count, double roundto, uint32_t *rtn_total, uint32_t *rtn_count)
Count the number of times provided value(s) occur in the band.
Definition: rt_statistics.c:1608
FALSE
#define FALSE
Definition: dbfopen.c:168
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2201
free
void free(void *)
rt_raster_t
Definition: librtcore.h:2235
genraster.value
int value
Definition: genraster.py:61
POSTGIS_RT_DEBUG
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h:61
TRUE
#define TRUE
Definition: dbfopen.c:169
rt_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:717
rt_valuecount_t
Definition: librtcore.h:2374
Here is the call graph for this function:
Here is the caller graph for this function: