PostGIS  2.1.10dev-r@@SVN_REVISION@@
Datum RASTER_nMapAlgebraExpr ( PG_FUNCTION_ARGS  )

Definition at line 17332 of file rt_pg.c.

References ovdump::band, rtpg_nmapalgebraexpr_arg_t::bandarg, rtpg_nmapalgebraexpr_arg_t::callback, rtpg_nmapalgebra_arg_t::cextent, ES_NONE, ET_CUSTOM, ET_FIRST, ET_INTERSECTION, ET_LAST, ET_SECOND, rtpg_nmapalgebraexpr_callback_arg::expr, rtpg_nmapalgebraexpr_callback_arg::exprcount, rtpg_nmapalgebra_arg_t::extenttype, FALSE, rtpg_nmapalgebra_arg_t::hasband, rtpg_nmapalgebra_arg_t::hasnodata, rtpg_nmapalgebraexpr_callback_arg::hasval, rt_iterator_t::nband, rtpg_nmapalgebra_arg_t::nband, rt_iterator_t::nbnodata, rtpg_nmapalgebraexpr_callback_arg::nodatanodata, rtpg_nmapalgebra_arg_t::nodataval, rtpg_nmapalgebra_arg_t::numraster, rtpg_nmapalgebra_arg_t::pixtype, POSTGIS_RT_DEBUGF, PT_END, rtrowdump::raster, rt_iterator_t::raster, rtpg_nmapalgebra_arg_t::raster, rt_band_get_hasnodata_flag(), rt_band_get_min_value(), rt_band_get_nodata(), rt_band_get_pixtype(), rt_pixtype_index_from_name(), rt_pixtype_name(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_iterator(), rt_raster_new(), rt_raster_serialize(), rt_util_extent_type(), rtpg_nmapalgebra_rastbandarg_process(), rtpg_nmapalgebraexpr_arg_destroy(), rtpg_nmapalgebraexpr_arg_init(), rtpg_nmapalgebraexpr_callback(), rtpg_strreplace(), rtpg_strtoupper(), rtpg_trim(), rt_raster_serialized_t::size, rtpg_nmapalgebraexpr_callback_arg::spi_argcount, rtpg_nmapalgebraexpr_callback_arg::spi_argpos, rtpg_nmapalgebraexpr_callback_arg::spi_plan, rtgdalraster::sql, TRUE, and rtpg_nmapalgebraexpr_callback_arg::val.

17333 {
17334  MemoryContext mainMemCtx = CurrentMemoryContext;
17335  rtpg_nmapalgebraexpr_arg arg = NULL;
17336  rt_iterator itrset;
17337  uint16_t exprpos[3] = {1, 4, 5};
17338 
17339  int i = 0;
17340  int j = 0;
17341  int k = 0;
17342 
17343  int numraster = 0;
17344  int err = 0;
17345  int allnull = 0;
17346  int allempty = 0;
17347  int noband = 0;
17348  int len = 0;
17349 
17350  TupleDesc tupdesc;
17351  SPITupleTable *tuptable = NULL;
17352  HeapTuple tuple;
17353  Datum datum;
17354  bool isnull = FALSE;
17355 
17356  rt_raster raster = NULL;
17357  rt_band band = NULL;
17358  rt_pgraster *pgraster = NULL;
17359 
17360  const int argkwcount = 12;
17361  char *argkw[] = {
17362  "[rast.x]",
17363  "[rast.y]",
17364  "[rast.val]",
17365  "[rast]",
17366  "[rast1.x]",
17367  "[rast1.y]",
17368  "[rast1.val]",
17369  "[rast1]",
17370  "[rast2.x]",
17371  "[rast2.y]",
17372  "[rast2.val]",
17373  "[rast2]"
17374  };
17375 
17376  if (PG_ARGISNULL(0))
17377  PG_RETURN_NULL();
17378 
17379  /* init argument struct */
17380  arg = rtpg_nmapalgebraexpr_arg_init(argkwcount, argkw);
17381  if (arg == NULL) {
17382  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not initialize argument structure");
17383  PG_RETURN_NULL();
17384  }
17385 
17386  /* let helper function process rastbandarg (0) */
17387  if (!rtpg_nmapalgebra_rastbandarg_process(arg->bandarg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
17389  elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
17390  PG_RETURN_NULL();
17391  }
17392 
17393  POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
17394 
17395  /* all rasters are NULL, return NULL */
17396  if (allnull == arg->bandarg->numraster) {
17397  elog(NOTICE, "All input rasters are NULL. Returning NULL");
17399  PG_RETURN_NULL();
17400  }
17401 
17402  /* only work on one or two rasters */
17403  if (arg->bandarg->numraster > 1)
17404  numraster = 2;
17405  else
17406  numraster = 1;
17407 
17408  /* pixel type (2) */
17409  if (!PG_ARGISNULL(2)) {
17410  char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
17411 
17412  /* Get the pixel type index */
17413  arg->bandarg->pixtype = rt_pixtype_index_from_name(pixtypename);
17414  if (arg->bandarg->pixtype == PT_END) {
17416  elog(ERROR, "RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
17417  PG_RETURN_NULL();
17418  }
17419  }
17420  POSTGIS_RT_DEBUGF(4, "pixeltype: %d", arg->bandarg->pixtype);
17421 
17422  /* extent type (3) */
17423  if (!PG_ARGISNULL(3)) {
17424  char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(3))));
17425  arg->bandarg->extenttype = rt_util_extent_type(extenttypename);
17426  }
17427 
17428  if (arg->bandarg->extenttype == ET_CUSTOM) {
17429  if (numraster < 2) {
17430  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to FIRST");
17431  arg->bandarg->extenttype = ET_FIRST;
17432  }
17433  else {
17434  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to INTERSECTION");
17436  }
17437  }
17438  else if (numraster < 2)
17439  arg->bandarg->extenttype = ET_FIRST;
17440 
17441  POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->bandarg->extenttype);
17442 
17443  /* nodatanodataval (6) */
17444  if (!PG_ARGISNULL(6)) {
17445  arg->callback.nodatanodata.hasval = 1;
17446  arg->callback.nodatanodata.val = PG_GETARG_FLOAT8(6);
17447  }
17448 
17449  err = 0;
17450  /* all rasters are empty, return empty raster */
17451  if (allempty == arg->bandarg->numraster) {
17452  elog(NOTICE, "All input rasters are empty. Returning empty raster");
17453  err = 1;
17454  }
17455  /* all rasters don't have indicated band, return empty raster */
17456  else if (noband == arg->bandarg->numraster) {
17457  elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
17458  err = 1;
17459  }
17460  if (err) {
17462 
17463  raster = rt_raster_new(0, 0);
17464  if (raster == NULL) {
17465  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create empty raster");
17466  PG_RETURN_NULL();
17467  }
17468 
17469  pgraster = rt_raster_serialize(raster);
17470  rt_raster_destroy(raster);
17471  if (!pgraster) PG_RETURN_NULL();
17472 
17473  SET_VARSIZE(pgraster, pgraster->size);
17474  PG_RETURN_POINTER(pgraster);
17475  }
17476 
17477  /* connect SPI */
17478  if (SPI_connect() != SPI_OK_CONNECT) {
17480  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
17481  PG_RETURN_NULL();
17482  }
17483 
17484  /*
17485  process expressions
17486 
17487  exprpos elements are:
17488  1 - expression => spi_plan[0]
17489  4 - nodata1expr => spi_plan[1]
17490  5 - nodata2expr => spi_plan[2]
17491  */
17492  for (i = 0; i < arg->callback.exprcount; i++) {
17493  char *expr = NULL;
17494  char *tmp = NULL;
17495  char *sql = NULL;
17496  char place[5] = "$1";
17497 
17498  if (PG_ARGISNULL(exprpos[i]))
17499  continue;
17500 
17501  expr = text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]));
17502  POSTGIS_RT_DEBUGF(3, "raw expr of argument #%d: %s", exprpos[i], expr);
17503 
17504  for (j = 0, k = 1; j < argkwcount; j++) {
17505  /* attempt to replace keyword with placeholder */
17506  len = 0;
17507  tmp = rtpg_strreplace(expr, argkw[j], place, &len);
17508  pfree(expr);
17509  expr = tmp;
17510 
17511  if (len) {
17512  POSTGIS_RT_DEBUGF(4, "kw #%d (%s) at pos $%d", j, argkw[j], k);
17513  arg->callback.expr[i].spi_argcount++;
17514  arg->callback.expr[i].spi_argpos[j] = k++;
17515 
17516  sprintf(place, "$%d", k);
17517  }
17518  else
17519  arg->callback.expr[i].spi_argpos[j] = 0;
17520  }
17521 
17522  len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
17523  sql = (char *) palloc(len + 1);
17524  if (sql == NULL) {
17526  SPI_finish();
17527  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
17528  PG_RETURN_NULL();
17529  }
17530 
17531  strncpy(sql, "SELECT (", strlen("SELECT ("));
17532  strncpy(sql + strlen("SELECT ("), expr, strlen(expr));
17533  strncpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
17534  sql[len] = '\0';
17535 
17536  POSTGIS_RT_DEBUGF(3, "sql #%d: %s", exprpos[i], sql);
17537 
17538  /* prepared plan */
17539  if (arg->callback.expr[i].spi_argcount) {
17540  Oid *argtype = (Oid *) palloc(arg->callback.expr[i].spi_argcount * sizeof(Oid));
17541  POSTGIS_RT_DEBUGF(3, "expression parameter %d is a prepared plan", exprpos[i]);
17542  if (argtype == NULL) {
17543  pfree(sql);
17545  SPI_finish();
17546  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
17547  PG_RETURN_NULL();
17548  }
17549 
17550  /* specify datatypes of parameters */
17551  for (j = 0, k = 0; j < argkwcount; j++) {
17552  if (arg->callback.expr[i].spi_argpos[j] < 1) continue;
17553 
17554  /* positions are INT4 */
17555  if (
17556  (strstr(argkw[j], "[rast.x]") != NULL) ||
17557  (strstr(argkw[j], "[rast.y]") != NULL) ||
17558  (strstr(argkw[j], "[rast1.x]") != NULL) ||
17559  (strstr(argkw[j], "[rast1.y]") != NULL) ||
17560  (strstr(argkw[j], "[rast2.x]") != NULL) ||
17561  (strstr(argkw[j], "[rast2.y]") != NULL)
17562  )
17563  argtype[k] = INT4OID;
17564  /* everything else is FLOAT8 */
17565  else
17566  argtype[k] = FLOAT8OID;
17567 
17568  k++;
17569  }
17570 
17571  arg->callback.expr[i].spi_plan = SPI_prepare(sql, arg->callback.expr[i].spi_argcount, argtype);
17572  pfree(argtype);
17573  pfree(sql);
17574 
17575  if (arg->callback.expr[i].spi_plan == NULL) {
17577  SPI_finish();
17578  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
17579  PG_RETURN_NULL();
17580  }
17581  }
17582  /* no args, just execute query */
17583  else {
17584  POSTGIS_RT_DEBUGF(3, "expression parameter %d has no args, simply executing", exprpos[i]);
17585  err = SPI_execute(sql, TRUE, 0);
17586  pfree(sql);
17587 
17588  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
17590  SPI_finish();
17591  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
17592  PG_RETURN_NULL();
17593  }
17594 
17595  /* get output of prepared plan */
17596  tupdesc = SPI_tuptable->tupdesc;
17597  tuptable = SPI_tuptable;
17598  tuple = tuptable->vals[0];
17599 
17600  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
17601  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
17602  if (SPI_tuptable) SPI_freetuptable(tuptable);
17604  SPI_finish();
17605  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
17606  PG_RETURN_NULL();
17607  }
17608 
17609  if (!isnull) {
17610  arg->callback.expr[i].hasval = 1;
17611  arg->callback.expr[i].val = DatumGetFloat8(datum);
17612  }
17613 
17614  if (SPI_tuptable) SPI_freetuptable(tuptable);
17615  }
17616  }
17617 
17618  /* determine nodataval and possibly pixtype */
17619  /* band to check */
17620  switch (arg->bandarg->extenttype) {
17621  case ET_LAST:
17622  case ET_SECOND:
17623  if (numraster > 1)
17624  i = 1;
17625  else
17626  i = 0;
17627  break;
17628  default:
17629  i = 0;
17630  break;
17631  }
17632  /* find first viable band */
17633  if (!arg->bandarg->hasband[i]) {
17634  for (i = 0; i < numraster; i++) {
17635  if (arg->bandarg->hasband[i])
17636  break;
17637  }
17638  if (i >= numraster)
17639  i = numraster - 1;
17640  }
17641  band = rt_raster_get_band(arg->bandarg->raster[i], arg->bandarg->nband[i]);
17642 
17643  /* set pixel type if PT_END */
17644  if (arg->bandarg->pixtype == PT_END)
17645  arg->bandarg->pixtype = rt_band_get_pixtype(band);
17646 
17647  /* set hasnodata and nodataval */
17648  arg->bandarg->hasnodata = 1;
17649  if (rt_band_get_hasnodata_flag(band))
17650  rt_band_get_nodata(band, &(arg->bandarg->nodataval));
17651  else
17652  arg->bandarg->nodataval = rt_band_get_min_value(band);
17653 
17654  POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->bandarg->pixtype), arg->bandarg->hasnodata, arg->bandarg->nodataval);
17655 
17656  /* init itrset */
17657  itrset = palloc(sizeof(struct rt_iterator_t) * numraster);
17658  if (itrset == NULL) {
17660  SPI_finish();
17661  elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
17662  PG_RETURN_NULL();
17663  }
17664 
17665  /* set itrset */
17666  for (i = 0; i < numraster; i++) {
17667  itrset[i].raster = arg->bandarg->raster[i];
17668  itrset[i].nband = arg->bandarg->nband[i];
17669  itrset[i].nbnodata = 1;
17670  }
17671 
17672  /* pass everything to iterator */
17673  err = rt_raster_iterator(
17674  itrset, numraster,
17675  arg->bandarg->extenttype, arg->bandarg->cextent,
17676  arg->bandarg->pixtype,
17677  arg->bandarg->hasnodata, arg->bandarg->nodataval,
17678  0, 0,
17679  &(arg->callback),
17681  &raster
17682  );
17683 
17684  pfree(itrset);
17686 
17687  if (err != ES_NONE) {
17688  SPI_finish();
17689  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not run raster iterator function");
17690  PG_RETURN_NULL();
17691  }
17692  else if (raster == NULL) {
17693  SPI_finish();
17694  PG_RETURN_NULL();
17695  }
17696 
17697  /* switch to prior memory context to ensure memory allocated in correct context */
17698  MemoryContextSwitchTo(mainMemCtx);
17699 
17700  pgraster = rt_raster_serialize(raster);
17701  rt_raster_destroy(raster);
17702 
17703  /* finish SPI */
17704  SPI_finish();
17705 
17706  if (!pgraster)
17707  PG_RETURN_NULL();
17708 
17709  SET_VARSIZE(pgraster, pgraster->size);
17710  PG_RETURN_POINTER(pgraster);
17711 }
static char * rtpg_strtoupper(char *str)
Definition: rt_pg.c:730
static char * rtpg_strreplace(const char *str, const char *oldstr, const char *newstr, int *count)
Definition: rt_pg.c:685
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_api.c:8158
struct rtpg_nmapalgebraexpr_callback_arg::@17 expr[3]
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
static char * rtpg_trim(const char *input)
Definition: rt_pg.c:856
rtpg_nmapalgebra_arg bandarg
Definition: rt_pg.c:17049
const char * rt_pixtype_name(rt_pixtype pixtype)
Definition: rt_api.c:1168
struct rtpg_nmapalgebraexpr_callback_arg::@18 nodatanodata
Definition: rt_api.h:184
rt_extenttype rt_util_extent_type(const char *name)
Definition: rt_api.c:302
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_api.c:1900
rt_raster * raster
Definition: rt_pg.c:16283
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:123
static int rtpg_nmapalgebra_rastbandarg_process(rtpg_nmapalgebra_arg arg, ArrayType *array, int *allnull, int *allempty, int *noband)
Definition: rt_pg.c:16360
static rtpg_nmapalgebraexpr_arg rtpg_nmapalgebraexpr_arg_init(int cnt, char **kw)
Definition: rt_pg.c:17054
rt_pixtype pixtype
Definition: rt_pg.c:16289
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
Definition: rt_api.c:3073
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_api.c:3058
rt_raster raster
Definition: rt_api.h:2360
uint16_t nband
Definition: rt_api.h:2361
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
#define FALSE
Definition: dbfopen.c:169
uint8_t nbnodata
Definition: rt_api.h:2362
Struct definitions.
Definition: rt_api.h:2175
static int rtpg_nmapalgebraexpr_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
Definition: rt_pg.c:17108
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rt_pg.h:62
rt_errorstate rt_raster_iterator(rt_iterator itrset, uint16_t itrcount, rt_extenttype extenttype, rt_raster customextent, rt_pixtype pixtype, uint8_t hasnodata, double nodataval, uint16_t distancex, uint16_t distancey, void *userarg, int(*callback)(rt_iterator_arg arg, void *userarg, double *value, int *nodata), rt_raster *rtnraster)
n-raster iterator.
Definition: rt_api.c:13927
rtpg_nmapalgebraexpr_callback_arg callback
Definition: rt_pg.c:17051
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
Definition: rt_api.c:1138
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_api.c:5353
static void rtpg_nmapalgebraexpr_arg_destroy(rtpg_nmapalgebraexpr_arg arg)
Definition: rt_pg.c:17093
rt_extenttype extenttype
Definition: rt_pg.c:16295
#define TRUE
Definition: dbfopen.c:170
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_api.c:2002

Here is the call graph for this function: