PostGIS  2.5.0beta2dev-r@@SVN_REVISION@@

◆ RASTER_valueCountCoverage()

Datum RASTER_valueCountCoverage ( PG_FUNCTION_ARGS  )

Definition at line 2607 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, text_to_cstring(), TRUE, genraster::value, and rt_valuecount_t::value.

Referenced by RASTER_valueCount().

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