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

Create a new empty raster with having the same georeference as the provided raster

If this new raster is empty (width = 0 OR height = 0) then there is nothing to compute and we return it right now

Check if the raster has the required band. Otherwise, return a raster without band

We set the initial value of the future band to nodata value. If nodata value is null, then the raster will be initialized to rt_band_get_min_value but all the values should be recomputed anyway

Set the new pixeltype

Optimization: If the raster is only filled with nodata values return right now a raster filled with the nodatavalueexpr TODO: Call rt_band_check_isnodata instead?

Create the raster receiving all the computed values. Initialize it to the new initial value

We compute a value only for the withdata value pixel since the nodata value has already been set by the first optimization

Definition at line 7427 of file rt_pg.c.

References ovdump::band, ES_NONE, FALSE, FLT_EQ, rt_raster_serialized_t::height, pixval::nband, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, PT_END, r, rtrowdump::raster, rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_band_get_min_value(), rt_band_get_nodata(), rt_band_get_pixel(), rt_band_get_pixtype(), rt_band_set_pixel(), rt_pixtype_index_from_name(), rt_pixtype_name(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_generate_new_band(), rt_raster_get_band(), rt_raster_get_height(), rt_raster_get_srid(), rt_raster_get_width(), rt_raster_get_x_offset(), rt_raster_get_x_scale(), rt_raster_get_x_skew(), rt_raster_get_y_offset(), rt_raster_get_y_scale(), rt_raster_get_y_skew(), rt_raster_has_band(), rt_raster_is_empty(), rt_raster_new(), rt_raster_serialize(), rt_raster_set_offsets(), rt_raster_set_scale(), rt_raster_set_skews(), rt_raster_set_srid(), rt_raster_serialized_t::size, TRUE, rt_raster_serialized_t::width, pixval::x, and pixval::y.

7428 {
7429  rt_pgraster *pgraster = NULL;
7430  rt_pgraster *pgrtn = NULL;
7431  rt_raster raster = NULL;
7432  rt_raster newrast = NULL;
7433  rt_band band = NULL;
7434  rt_band newband = NULL;
7435  int x, y, nband, width, height;
7436  double r;
7437  double newnodatavalue = 0.0;
7438  double newinitialvalue = 0.0;
7439  double newval = 0.0;
7440  rt_pixtype newpixeltype;
7441  int ret = -1;
7442  Oid oid;
7443  FmgrInfo cbinfo;
7444  FunctionCallInfoData cbdata;
7445  Datum tmpnewval;
7446  char * strFromText = NULL;
7447  int k = 0;
7448 
7449  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFct: STARTING...");
7450 
7451  /* Check raster */
7452  if (PG_ARGISNULL(0)) {
7453  elog(WARNING, "Raster is NULL. Returning NULL");
7454  PG_RETURN_NULL();
7455  }
7456 
7457 
7458  /* Deserialize raster */
7459  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
7460  raster = rt_raster_deserialize(pgraster, FALSE);
7461  if (NULL == raster) {
7462  PG_FREE_IF_COPY(pgraster, 0);
7463  elog(ERROR, "RASTER_mapAlgebraFct: Could not deserialize raster");
7464  PG_RETURN_NULL();
7465  }
7466 
7467  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting arguments...");
7468 
7469  /* Get the rest of the arguments */
7470 
7471  if (PG_ARGISNULL(1))
7472  nband = 1;
7473  else
7474  nband = PG_GETARG_INT32(1);
7475 
7476  if (nband < 1)
7477  nband = 1;
7478 
7479  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Creating new empty raster...");
7480 
7485  width = rt_raster_get_width(raster);
7486  height = rt_raster_get_height(raster);
7487 
7488  newrast = rt_raster_new(width, height);
7489 
7490  if ( NULL == newrast ) {
7491 
7492  rt_raster_destroy(raster);
7493  PG_FREE_IF_COPY(pgraster, 0);
7494 
7495  elog(ERROR, "RASTER_mapAlgebraFct: Could not create a new raster");
7496  PG_RETURN_NULL();
7497  }
7498 
7499  rt_raster_set_scale(newrast,
7500  rt_raster_get_x_scale(raster),
7501  rt_raster_get_y_scale(raster));
7502 
7503  rt_raster_set_offsets(newrast,
7504  rt_raster_get_x_offset(raster),
7505  rt_raster_get_y_offset(raster));
7506 
7507  rt_raster_set_skews(newrast,
7508  rt_raster_get_x_skew(raster),
7509  rt_raster_get_y_skew(raster));
7510 
7511  rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
7512 
7513 
7518  if (rt_raster_is_empty(newrast))
7519  {
7520  elog(NOTICE, "Raster is empty. Returning an empty raster");
7521  rt_raster_destroy(raster);
7522  PG_FREE_IF_COPY(pgraster, 0);
7523 
7524  pgrtn = rt_raster_serialize(newrast);
7525  rt_raster_destroy(newrast);
7526  if (NULL == pgrtn) {
7527  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
7528  PG_RETURN_NULL();
7529  }
7530 
7531  SET_VARSIZE(pgrtn, pgrtn->size);
7532  PG_RETURN_POINTER(pgrtn);
7533  }
7534 
7535  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Getting raster band %d...", nband);
7536 
7541  if (!rt_raster_has_band(raster, nband - 1)) {
7542  elog(NOTICE, "Raster does not have the required band. Returning a raster "
7543  "without a band");
7544  rt_raster_destroy(raster);
7545  PG_FREE_IF_COPY(pgraster, 0);
7546 
7547  pgrtn = rt_raster_serialize(newrast);
7548  rt_raster_destroy(newrast);
7549  if (NULL == pgrtn) {
7550  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
7551  PG_RETURN_NULL();
7552  }
7553 
7554  SET_VARSIZE(pgrtn, pgrtn->size);
7555  PG_RETURN_POINTER(pgrtn);
7556  }
7557 
7558  /* Get the raster band */
7559  band = rt_raster_get_band(raster, nband - 1);
7560  if ( NULL == band ) {
7561  elog(NOTICE, "Could not get the required band. Returning a raster "
7562  "without a band");
7563  rt_raster_destroy(raster);
7564  PG_FREE_IF_COPY(pgraster, 0);
7565 
7566  pgrtn = rt_raster_serialize(newrast);
7567  rt_raster_destroy(newrast);
7568  if (NULL == pgrtn) {
7569  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
7570  PG_RETURN_NULL();
7571  }
7572 
7573  SET_VARSIZE(pgrtn, pgrtn->size);
7574  PG_RETURN_POINTER(pgrtn);
7575  }
7576 
7577  /*
7578  * Get NODATA value
7579  */
7580  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting NODATA value for band...");
7581 
7582  if (rt_band_get_hasnodata_flag(band)) {
7583  rt_band_get_nodata(band, &newnodatavalue);
7584  }
7585 
7586  else {
7587  newnodatavalue = rt_band_get_min_value(band);
7588  }
7589 
7590  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: NODATA value for band: %f",
7591  newnodatavalue);
7597  newinitialvalue = newnodatavalue;
7598 
7602  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Setting pixeltype...");
7603 
7604  if (PG_ARGISNULL(2)) {
7605  newpixeltype = rt_band_get_pixtype(band);
7606  }
7607 
7608  else {
7609  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
7610  newpixeltype = rt_pixtype_index_from_name(strFromText);
7611  pfree(strFromText);
7612  if (newpixeltype == PT_END)
7613  newpixeltype = rt_band_get_pixtype(band);
7614  }
7615 
7616  if (newpixeltype == PT_END) {
7617 
7618  rt_raster_destroy(raster);
7619  PG_FREE_IF_COPY(pgraster, 0);
7620  rt_raster_destroy(newrast);
7621 
7622  elog(ERROR, "RASTER_mapAlgebraFct: Invalid pixeltype");
7623  PG_RETURN_NULL();
7624  }
7625 
7626  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Pixeltype set to %s",
7627  rt_pixtype_name(newpixeltype));
7628 
7629  /* Get the name of the callback user function for raster values */
7630  if (PG_ARGISNULL(3)) {
7631 
7632  rt_raster_destroy(raster);
7633  PG_FREE_IF_COPY(pgraster, 0);
7634  rt_raster_destroy(newrast);
7635 
7636  elog(ERROR, "RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
7637  PG_RETURN_NULL();
7638  }
7639 
7640  oid = PG_GETARG_OID(3);
7641  if (oid == InvalidOid) {
7642 
7643  rt_raster_destroy(raster);
7644  PG_FREE_IF_COPY(pgraster, 0);
7645  rt_raster_destroy(newrast);
7646 
7647  elog(ERROR, "RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
7648  PG_RETURN_NULL();
7649  }
7650 
7651  fmgr_info(oid, &cbinfo);
7652 
7653  /* function cannot return set */
7654  if (cbinfo.fn_retset) {
7655 
7656  rt_raster_destroy(raster);
7657  PG_FREE_IF_COPY(pgraster, 0);
7658  rt_raster_destroy(newrast);
7659 
7660  elog(ERROR, "RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
7661  PG_RETURN_NULL();
7662  }
7663  /* function should have correct # of args */
7664  else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
7665 
7666  rt_raster_destroy(raster);
7667  PG_FREE_IF_COPY(pgraster, 0);
7668  rt_raster_destroy(newrast);
7669 
7670  elog(ERROR, "RASTER_mapAlgebraFct: Function does not have two or three input parameters");
7671  PG_RETURN_NULL();
7672  }
7673 
7674  if (cbinfo.fn_nargs == 2)
7675  k = 1;
7676  else
7677  k = 2;
7678 
7679  if (func_volatile(oid) == 'v') {
7680  elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
7681  }
7682 
7683  /* prep function call data */
7684 #if POSTGIS_PGSQL_VERSION <= 90
7685  InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL);
7686 #else
7687  InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
7688 #endif
7689  memset(cbdata.argnull, FALSE, sizeof(bool) * cbinfo.fn_nargs);
7690 
7691  /* check that the function isn't strict if the args are null. */
7692  if (PG_ARGISNULL(4)) {
7693  if (cbinfo.fn_strict) {
7694 
7695  rt_raster_destroy(raster);
7696  PG_FREE_IF_COPY(pgraster, 0);
7697  rt_raster_destroy(newrast);
7698 
7699  elog(ERROR, "RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
7700  PG_RETURN_NULL();
7701  }
7702 
7703  cbdata.arg[k] = (Datum)NULL;
7704  cbdata.argnull[k] = TRUE;
7705  }
7706  else {
7707  cbdata.arg[k] = PG_GETARG_DATUM(4);
7708  }
7709 
7715  if (rt_band_get_isnodata_flag(band)) {
7716 
7717  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Band is a nodata band, returning "
7718  "a raster filled with nodata");
7719 
7720  ret = rt_raster_generate_new_band(newrast, newpixeltype,
7721  newinitialvalue, TRUE, newnodatavalue, 0);
7722 
7723  rt_raster_destroy(raster);
7724  PG_FREE_IF_COPY(pgraster, 0);
7725 
7726  /* Serialize created raster */
7727  pgrtn = rt_raster_serialize(newrast);
7728  rt_raster_destroy(newrast);
7729  if (NULL == pgrtn) {
7730  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
7731  PG_RETURN_NULL();
7732  }
7733 
7734  SET_VARSIZE(pgrtn, pgrtn->size);
7735  PG_RETURN_POINTER(pgrtn);
7736  }
7737 
7738 
7743  ret = rt_raster_generate_new_band(newrast, newpixeltype,
7744  newinitialvalue, TRUE, newnodatavalue, 0);
7745 
7746  /* Get the new raster band */
7747  newband = rt_raster_get_band(newrast, 0);
7748  if ( NULL == newband ) {
7749  elog(NOTICE, "Could not modify band for new raster. Returning new "
7750  "raster with the original band");
7751 
7752  rt_raster_destroy(raster);
7753  PG_FREE_IF_COPY(pgraster, 0);
7754 
7755  /* Serialize created raster */
7756  pgrtn = rt_raster_serialize(newrast);
7757  rt_raster_destroy(newrast);
7758  if (NULL == pgrtn) {
7759  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
7760  PG_RETURN_NULL();
7761  }
7762 
7763  SET_VARSIZE(pgrtn, pgrtn->size);
7764  PG_RETURN_POINTER(pgrtn);
7765  }
7766 
7767  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Main computing loop (%d x %d)",
7768  width, height);
7769 
7770  for (x = 0; x < width; x++) {
7771  for(y = 0; y < height; y++) {
7772  ret = rt_band_get_pixel(band, x, y, &r, NULL);
7773 
7778  if (ret == ES_NONE) {
7779  if (FLT_EQ(r, newnodatavalue)) {
7780  if (cbinfo.fn_strict) {
7781  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
7782  continue;
7783  }
7784  cbdata.argnull[0] = TRUE;
7785  cbdata.arg[0] = (Datum)NULL;
7786  }
7787  else {
7788  cbdata.argnull[0] = FALSE;
7789  cbdata.arg[0] = Float8GetDatum(r);
7790  }
7791 
7792  /* Add pixel positions if callback has proper # of args */
7793  if (cbinfo.fn_nargs == 3) {
7794  Datum d[2];
7795  ArrayType *a;
7796 
7797  d[0] = Int32GetDatum(x+1);
7798  d[1] = Int32GetDatum(y+1);
7799 
7800  a = construct_array(d, 2, INT4OID, sizeof(int32), true, 'i');
7801 
7802  cbdata.argnull[1] = FALSE;
7803  cbdata.arg[1] = PointerGetDatum(a);
7804  }
7805 
7806  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: (%dx%d), r = %f",
7807  x, y, r);
7808 
7809  tmpnewval = FunctionCallInvoke(&cbdata);
7810 
7811  if (cbdata.isnull) {
7812  newval = newnodatavalue;
7813  }
7814  else {
7815  newval = DatumGetFloat8(tmpnewval);
7816  }
7817 
7818  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: new value = %f",
7819  newval);
7820 
7821  rt_band_set_pixel(newband, x, y, newval, NULL);
7822  }
7823 
7824  }
7825  }
7826 
7827  /* The newrast band has been modified */
7828 
7829  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster modified, serializing it.");
7830  /* Serialize created raster */
7831 
7832  rt_raster_destroy(raster);
7833  PG_FREE_IF_COPY(pgraster, 0);
7834 
7835  pgrtn = rt_raster_serialize(newrast);
7836  rt_raster_destroy(newrast);
7837  if (NULL == pgrtn)
7838  PG_RETURN_NULL();
7839 
7840  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster serialized");
7841 
7842  POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFct: returning raster");
7843 
7844  SET_VARSIZE(pgrtn, pgrtn->size);
7845  PG_RETURN_POINTER(pgrtn);
7846 }
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_api.c:8158
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition: rt_api.c:5527
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
char * r
Definition: cu_in_wkt.c:25
const char * rt_pixtype_name(rt_pixtype pixtype)
Definition: rt_api.c:1168
Definition: rt_api.h:184
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_api.c:1900
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_api.c:5486
tuple band
Definition: ovdump.py:57
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
Definition: rt_api.c:8563
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition: rt_api.c:5495
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition: rt_api.c:5518
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:123
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_api.c:5668
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_api.c:5661
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_api.c:2042
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_api.c:5473
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_api.c:5442
rt_pixtype
Definition: rt_api.h:172
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rt_pg.h:58
tuple nband
Definition: pixval.py:52
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_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_api.c:2549
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_api.c:5455
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
tuple x
Definition: pixval.py:53
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_api.c:5434
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
Struct definitions.
Definition: rt_api.h:2175
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rt_pg.h:62
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_api.c:8550
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_api.c:5464
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_api.c:8350
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
#define TRUE
Definition: dbfopen.c:170
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_api.c:2002
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_api.c:5504
tuple y
Definition: pixval.py:54
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_api.c:5426
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
Definition: rt_api.c:5784
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
Definition: rt_api.c:2302

Here is the call graph for this function: