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

Definition at line 1273 of file rtpg_mapalgebra.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, rtpg_nmapalgebra_arg_t::nband, rt_iterator_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, rtpg_nmapalgebra_arg_t::raster, rtrowdump::raster, rt_iterator_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.

1274 {
1275  MemoryContext mainMemCtx = CurrentMemoryContext;
1276  rtpg_nmapalgebraexpr_arg arg = NULL;
1277  rt_iterator itrset;
1278  uint16_t exprpos[3] = {1, 4, 5};
1279 
1280  int i = 0;
1281  int j = 0;
1282  int k = 0;
1283 
1284  int numraster = 0;
1285  int err = 0;
1286  int allnull = 0;
1287  int allempty = 0;
1288  int noband = 0;
1289  int len = 0;
1290 
1291  TupleDesc tupdesc;
1292  SPITupleTable *tuptable = NULL;
1293  HeapTuple tuple;
1294  Datum datum;
1295  bool isnull = FALSE;
1296 
1297  rt_raster raster = NULL;
1298  rt_band band = NULL;
1299  rt_pgraster *pgraster = NULL;
1300 
1301  const int argkwcount = 12;
1302  char *argkw[] = {
1303  "[rast.x]",
1304  "[rast.y]",
1305  "[rast.val]",
1306  "[rast]",
1307  "[rast1.x]",
1308  "[rast1.y]",
1309  "[rast1.val]",
1310  "[rast1]",
1311  "[rast2.x]",
1312  "[rast2.y]",
1313  "[rast2.val]",
1314  "[rast2]"
1315  };
1316 
1317  if (PG_ARGISNULL(0))
1318  PG_RETURN_NULL();
1319 
1320  /* init argument struct */
1321  arg = rtpg_nmapalgebraexpr_arg_init(argkwcount, argkw);
1322  if (arg == NULL) {
1323  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not initialize argument structure");
1324  PG_RETURN_NULL();
1325  }
1326 
1327  /* let helper function process rastbandarg (0) */
1328  if (!rtpg_nmapalgebra_rastbandarg_process(arg->bandarg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
1330  elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
1331  PG_RETURN_NULL();
1332  }
1333 
1334  POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
1335 
1336  /* all rasters are NULL, return NULL */
1337  if (allnull == arg->bandarg->numraster) {
1338  elog(NOTICE, "All input rasters are NULL. Returning NULL");
1340  PG_RETURN_NULL();
1341  }
1342 
1343  /* only work on one or two rasters */
1344  if (arg->bandarg->numraster > 1)
1345  numraster = 2;
1346  else
1347  numraster = 1;
1348 
1349  /* pixel type (2) */
1350  if (!PG_ARGISNULL(2)) {
1351  char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
1352 
1353  /* Get the pixel type index */
1354  arg->bandarg->pixtype = rt_pixtype_index_from_name(pixtypename);
1355  if (arg->bandarg->pixtype == PT_END) {
1357  elog(ERROR, "RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
1358  PG_RETURN_NULL();
1359  }
1360  }
1361  POSTGIS_RT_DEBUGF(4, "pixeltype: %d", arg->bandarg->pixtype);
1362 
1363  /* extent type (3) */
1364  if (!PG_ARGISNULL(3)) {
1365  char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(3))));
1366  arg->bandarg->extenttype = rt_util_extent_type(extenttypename);
1367  }
1368 
1369  if (arg->bandarg->extenttype == ET_CUSTOM) {
1370  if (numraster < 2) {
1371  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to FIRST");
1372  arg->bandarg->extenttype = ET_FIRST;
1373  }
1374  else {
1375  elog(NOTICE, "CUSTOM extent type not supported. Defaulting to INTERSECTION");
1377  }
1378  }
1379  else if (numraster < 2)
1380  arg->bandarg->extenttype = ET_FIRST;
1381 
1382  POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->bandarg->extenttype);
1383 
1384  /* nodatanodataval (6) */
1385  if (!PG_ARGISNULL(6)) {
1386  arg->callback.nodatanodata.hasval = 1;
1387  arg->callback.nodatanodata.val = PG_GETARG_FLOAT8(6);
1388  }
1389 
1390  err = 0;
1391  /* all rasters are empty, return empty raster */
1392  if (allempty == arg->bandarg->numraster) {
1393  elog(NOTICE, "All input rasters are empty. Returning empty raster");
1394  err = 1;
1395  }
1396  /* all rasters don't have indicated band, return empty raster */
1397  else if (noband == arg->bandarg->numraster) {
1398  elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
1399  err = 1;
1400  }
1401  if (err) {
1403 
1404  raster = rt_raster_new(0, 0);
1405  if (raster == NULL) {
1406  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create empty raster");
1407  PG_RETURN_NULL();
1408  }
1409 
1410  pgraster = rt_raster_serialize(raster);
1411  rt_raster_destroy(raster);
1412  if (!pgraster) PG_RETURN_NULL();
1413 
1414  SET_VARSIZE(pgraster, pgraster->size);
1415  PG_RETURN_POINTER(pgraster);
1416  }
1417 
1418  /* connect SPI */
1419  if (SPI_connect() != SPI_OK_CONNECT) {
1421  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
1422  PG_RETURN_NULL();
1423  }
1424 
1425  /*
1426  process expressions
1427 
1428  exprpos elements are:
1429  1 - expression => spi_plan[0]
1430  4 - nodata1expr => spi_plan[1]
1431  5 - nodata2expr => spi_plan[2]
1432  */
1433  for (i = 0; i < arg->callback.exprcount; i++) {
1434  char *expr = NULL;
1435  char *tmp = NULL;
1436  char *sql = NULL;
1437  char place[5] = "$1";
1438 
1439  if (PG_ARGISNULL(exprpos[i]))
1440  continue;
1441 
1442  expr = text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]));
1443  POSTGIS_RT_DEBUGF(3, "raw expr of argument #%d: %s", exprpos[i], expr);
1444 
1445  for (j = 0, k = 1; j < argkwcount; j++) {
1446  /* attempt to replace keyword with placeholder */
1447  len = 0;
1448  tmp = rtpg_strreplace(expr, argkw[j], place, &len);
1449  pfree(expr);
1450  expr = tmp;
1451 
1452  if (len) {
1453  POSTGIS_RT_DEBUGF(4, "kw #%d (%s) at pos $%d", j, argkw[j], k);
1454  arg->callback.expr[i].spi_argcount++;
1455  arg->callback.expr[i].spi_argpos[j] = k++;
1456 
1457  sprintf(place, "$%d", k);
1458  }
1459  else
1460  arg->callback.expr[i].spi_argpos[j] = 0;
1461  }
1462 
1463  len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
1464  sql = (char *) palloc(len + 1);
1465  if (sql == NULL) {
1467  SPI_finish();
1468  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
1469  PG_RETURN_NULL();
1470  }
1471 
1472  strncpy(sql, "SELECT (", strlen("SELECT ("));
1473  strncpy(sql + strlen("SELECT ("), expr, strlen(expr));
1474  strncpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
1475  sql[len] = '\0';
1476 
1477  POSTGIS_RT_DEBUGF(3, "sql #%d: %s", exprpos[i], sql);
1478 
1479  /* prepared plan */
1480  if (arg->callback.expr[i].spi_argcount) {
1481  Oid *argtype = (Oid *) palloc(arg->callback.expr[i].spi_argcount * sizeof(Oid));
1482  POSTGIS_RT_DEBUGF(3, "expression parameter %d is a prepared plan", exprpos[i]);
1483  if (argtype == NULL) {
1484  pfree(sql);
1486  SPI_finish();
1487  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
1488  PG_RETURN_NULL();
1489  }
1490 
1491  /* specify datatypes of parameters */
1492  for (j = 0, k = 0; j < argkwcount; j++) {
1493  if (arg->callback.expr[i].spi_argpos[j] < 1) continue;
1494 
1495  /* positions are INT4 */
1496  if (
1497  (strstr(argkw[j], "[rast.x]") != NULL) ||
1498  (strstr(argkw[j], "[rast.y]") != NULL) ||
1499  (strstr(argkw[j], "[rast1.x]") != NULL) ||
1500  (strstr(argkw[j], "[rast1.y]") != NULL) ||
1501  (strstr(argkw[j], "[rast2.x]") != NULL) ||
1502  (strstr(argkw[j], "[rast2.y]") != NULL)
1503  )
1504  argtype[k] = INT4OID;
1505  /* everything else is FLOAT8 */
1506  else
1507  argtype[k] = FLOAT8OID;
1508 
1509  k++;
1510  }
1511 
1512  arg->callback.expr[i].spi_plan = SPI_prepare(sql, arg->callback.expr[i].spi_argcount, argtype);
1513  pfree(argtype);
1514  pfree(sql);
1515 
1516  if (arg->callback.expr[i].spi_plan == NULL) {
1518  SPI_finish();
1519  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
1520  PG_RETURN_NULL();
1521  }
1522  }
1523  /* no args, just execute query */
1524  else {
1525  POSTGIS_RT_DEBUGF(3, "expression parameter %d has no args, simply executing", exprpos[i]);
1526  err = SPI_execute(sql, TRUE, 0);
1527  pfree(sql);
1528 
1529  if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1531  SPI_finish();
1532  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
1533  PG_RETURN_NULL();
1534  }
1535 
1536  /* get output of prepared plan */
1537  tupdesc = SPI_tuptable->tupdesc;
1538  tuptable = SPI_tuptable;
1539  tuple = tuptable->vals[0];
1540 
1541  datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1542  if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1543  if (SPI_tuptable) SPI_freetuptable(tuptable);
1545  SPI_finish();
1546  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
1547  PG_RETURN_NULL();
1548  }
1549 
1550  if (!isnull) {
1551  arg->callback.expr[i].hasval = 1;
1552  arg->callback.expr[i].val = DatumGetFloat8(datum);
1553  }
1554 
1555  if (SPI_tuptable) SPI_freetuptable(tuptable);
1556  }
1557  }
1558 
1559  /* determine nodataval and possibly pixtype */
1560  /* band to check */
1561  switch (arg->bandarg->extenttype) {
1562  case ET_LAST:
1563  case ET_SECOND:
1564  if (numraster > 1)
1565  i = 1;
1566  else
1567  i = 0;
1568  break;
1569  default:
1570  i = 0;
1571  break;
1572  }
1573  /* find first viable band */
1574  if (!arg->bandarg->hasband[i]) {
1575  for (i = 0; i < numraster; i++) {
1576  if (arg->bandarg->hasband[i])
1577  break;
1578  }
1579  if (i >= numraster)
1580  i = numraster - 1;
1581  }
1582  band = rt_raster_get_band(arg->bandarg->raster[i], arg->bandarg->nband[i]);
1583 
1584  /* set pixel type if PT_END */
1585  if (arg->bandarg->pixtype == PT_END)
1586  arg->bandarg->pixtype = rt_band_get_pixtype(band);
1587 
1588  /* set hasnodata and nodataval */
1589  arg->bandarg->hasnodata = 1;
1590  if (rt_band_get_hasnodata_flag(band))
1591  rt_band_get_nodata(band, &(arg->bandarg->nodataval));
1592  else
1593  arg->bandarg->nodataval = rt_band_get_min_value(band);
1594 
1595  POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->bandarg->pixtype), arg->bandarg->hasnodata, arg->bandarg->nodataval);
1596 
1597  /* init itrset */
1598  itrset = palloc(sizeof(struct rt_iterator_t) * numraster);
1599  if (itrset == NULL) {
1601  SPI_finish();
1602  elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1603  PG_RETURN_NULL();
1604  }
1605 
1606  /* set itrset */
1607  for (i = 0; i < numraster; i++) {
1608  itrset[i].raster = arg->bandarg->raster[i];
1609  itrset[i].nband = arg->bandarg->nband[i];
1610  itrset[i].nbnodata = 1;
1611  }
1612 
1613  /* pass everything to iterator */
1614  err = rt_raster_iterator(
1615  itrset, numraster,
1616  arg->bandarg->extenttype, arg->bandarg->cextent,
1617  arg->bandarg->pixtype,
1618  arg->bandarg->hasnodata, arg->bandarg->nodataval,
1619  0, 0,
1620  NULL,
1621  &(arg->callback),
1623  &raster
1624  );
1625 
1626  pfree(itrset);
1628 
1629  if (err != ES_NONE) {
1630  SPI_finish();
1631  elog(ERROR, "RASTER_nMapAlgebraExpr: Could not run raster iterator function");
1632  PG_RETURN_NULL();
1633  }
1634  else if (raster == NULL) {
1635  SPI_finish();
1636  PG_RETURN_NULL();
1637  }
1638 
1639  /* switch to prior memory context to ensure memory allocated in correct context */
1640  MemoryContextSwitchTo(mainMemCtx);
1641 
1642  pgraster = rt_raster_serialize(raster);
1643  rt_raster_destroy(raster);
1644 
1645  /* finish SPI */
1646  SPI_finish();
1647 
1648  if (!pgraster)
1649  PG_RETURN_NULL();
1650 
1651  SET_VARSIZE(pgraster, pgraster->size);
1652  PG_RETURN_POINTER(pgraster);
1653 }
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_serialize.c:521
rtpg_nmapalgebra_arg bandarg
char * rtpg_strreplace(const char *str, const char *oldstr, const char *newstr, int *count)
Definition: rtpg_internal.c:55
static int rtpg_nmapalgebraexpr_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
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
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:57
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
struct rtpg_nmapalgebraexpr_callback_arg::@17 nodatanodata
rt_extenttype rt_util_extent_type(const char *name)
Definition: rt_util.c:191
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
Definition: rt_pixel.c:80
static rtpg_nmapalgebraexpr_arg rtpg_nmapalgebraexpr_arg_init(int cnt, char **kw)
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
static int rtpg_nmapalgebra_rastbandarg_process(rtpg_nmapalgebra_arg arg, ArrayType *array, int *allnull, int *allempty, int *noband)
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
char * rtpg_trim(const char *input)
struct rtpg_nmapalgebraexpr_callback_arg::@16 expr[3]
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
rt_raster raster
Definition: librtcore.h:2406
static void rtpg_nmapalgebraexpr_arg_destroy(rtpg_nmapalgebraexpr_arg arg)
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
uint16_t nband
Definition: librtcore.h:2407
#define FALSE
Definition: dbfopen.c:168
uint8_t nbnodata
Definition: librtcore.h:2408
Struct definitions.
Definition: librtcore.h:2213
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
Definition: rt_band.c:1612
const char * rt_pixtype_name(rt_pixtype pixtype)
Definition: rt_pixel.c:110
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, rt_mask mask, void *userarg, int(*callback)(rt_iterator_arg arg, void *userarg, double *value, int *nodata), rt_raster *rtnraster)
n-raster iterator.
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:498
rtpg_nmapalgebraexpr_callback_arg callback
#define TRUE
Definition: dbfopen.c:169
char * rtpg_strtoupper(char *str)

Here is the call graph for this function: