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

Definition at line 1270 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.

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