PostGIS  2.5.0dev-r@@SVN_REVISION@@
Datum RASTER_valueCountCoverage ( PG_FUNCTION_ARGS  )

Definition at line 2617 of file rtpg_statistics.c.

References ovdump::band, genraster::count, rt_valuecount_t::count, 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, and rt_valuecount_t::value.

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

Here is the call graph for this function: