PostGIS  3.1.6dev-r@@SVN_REVISION@@

◆ RASTER_mapAlgebraFct()

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 5143 of file rtpg_mapalgebra.c.

5144 {
5145  rt_pgraster *pgraster = NULL;
5146  rt_pgraster *pgrtn = NULL;
5147  rt_raster raster = NULL;
5148  rt_raster newrast = NULL;
5149  rt_band band = NULL;
5150  rt_band newband = NULL;
5151  int x, y, nband, width, height;
5152  double r;
5153  double newnodatavalue = 0.0;
5154  double newinitialvalue = 0.0;
5155  double newval = 0.0;
5156  rt_pixtype newpixeltype;
5157  int ret = -1;
5158  Oid oid;
5159  FmgrInfo cbinfo;
5160 #if POSTGIS_PGSQL_VERSION < 120
5161  FunctionCallInfoData cbdata;
5162 #else
5163  LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS); /* Could be optimized */
5164 #endif
5165  Datum tmpnewval;
5166  char * strFromText = NULL;
5167  int k = 0;
5168 
5169  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFct: STARTING...");
5170 
5171  /* Check raster */
5172  if (PG_ARGISNULL(0)) {
5173  elog(WARNING, "Raster is NULL. Returning NULL");
5174  PG_RETURN_NULL();
5175  }
5176 
5177 
5178  /* Deserialize raster */
5179  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5180  raster = rt_raster_deserialize(pgraster, FALSE);
5181  if (NULL == raster) {
5182  PG_FREE_IF_COPY(pgraster, 0);
5183  elog(ERROR, "RASTER_mapAlgebraFct: Could not deserialize raster");
5184  PG_RETURN_NULL();
5185  }
5186 
5187  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting arguments...");
5188 
5189  /* Get the rest of the arguments */
5190 
5191  if (PG_ARGISNULL(1))
5192  nband = 1;
5193  else
5194  nband = PG_GETARG_INT32(1);
5195 
5196  if (nband < 1)
5197  nband = 1;
5198 
5199  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Creating new empty raster...");
5200 
5205  width = rt_raster_get_width(raster);
5206  height = rt_raster_get_height(raster);
5207 
5208  newrast = rt_raster_new(width, height);
5209 
5210  if ( NULL == newrast ) {
5211 
5213  PG_FREE_IF_COPY(pgraster, 0);
5214 
5215  elog(ERROR, "RASTER_mapAlgebraFct: Could not create a new raster");
5216  PG_RETURN_NULL();
5217  }
5218 
5219  rt_raster_set_scale(newrast,
5222 
5223  rt_raster_set_offsets(newrast,
5226 
5227  rt_raster_set_skews(newrast,
5230 
5232 
5233 
5238  if (rt_raster_is_empty(newrast))
5239  {
5240  elog(NOTICE, "Raster is empty. Returning an empty raster");
5242  PG_FREE_IF_COPY(pgraster, 0);
5243 
5244  pgrtn = rt_raster_serialize(newrast);
5245  rt_raster_destroy(newrast);
5246  if (NULL == pgrtn) {
5247  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5248  PG_RETURN_NULL();
5249  }
5250 
5251  SET_VARSIZE(pgrtn, pgrtn->size);
5252  PG_RETURN_POINTER(pgrtn);
5253  }
5254 
5255  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Getting raster band %d...", nband);
5256 
5261  if (!rt_raster_has_band(raster, nband - 1)) {
5262  elog(NOTICE, "Raster does not have the required band. Returning a raster "
5263  "without a band");
5265  PG_FREE_IF_COPY(pgraster, 0);
5266 
5267  pgrtn = rt_raster_serialize(newrast);
5268  rt_raster_destroy(newrast);
5269  if (NULL == pgrtn) {
5270  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5271  PG_RETURN_NULL();
5272  }
5273 
5274  SET_VARSIZE(pgrtn, pgrtn->size);
5275  PG_RETURN_POINTER(pgrtn);
5276  }
5277 
5278  /* Get the raster band */
5280  if ( NULL == band ) {
5281  elog(NOTICE, "Could not get the required band. Returning a raster "
5282  "without a band");
5284  PG_FREE_IF_COPY(pgraster, 0);
5285 
5286  pgrtn = rt_raster_serialize(newrast);
5287  rt_raster_destroy(newrast);
5288  if (NULL == pgrtn) {
5289  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5290  PG_RETURN_NULL();
5291  }
5292 
5293  SET_VARSIZE(pgrtn, pgrtn->size);
5294  PG_RETURN_POINTER(pgrtn);
5295  }
5296 
5297  /*
5298  * Get NODATA value
5299  */
5300  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting NODATA value for band...");
5301 
5303  rt_band_get_nodata(band, &newnodatavalue);
5304  }
5305 
5306  else {
5307  newnodatavalue = rt_band_get_min_value(band);
5308  }
5309 
5310  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: NODATA value for band: %f",
5311  newnodatavalue);
5317  newinitialvalue = newnodatavalue;
5318 
5322  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Setting pixeltype...");
5323 
5324  if (PG_ARGISNULL(2)) {
5325  newpixeltype = rt_band_get_pixtype(band);
5326  }
5327 
5328  else {
5329  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5330  newpixeltype = rt_pixtype_index_from_name(strFromText);
5331  pfree(strFromText);
5332  if (newpixeltype == PT_END)
5333  newpixeltype = rt_band_get_pixtype(band);
5334  }
5335 
5336  if (newpixeltype == PT_END) {
5337 
5339  PG_FREE_IF_COPY(pgraster, 0);
5340  rt_raster_destroy(newrast);
5341 
5342  elog(ERROR, "RASTER_mapAlgebraFct: Invalid pixeltype");
5343  PG_RETURN_NULL();
5344  }
5345 
5346  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Pixeltype set to %s",
5347  rt_pixtype_name(newpixeltype));
5348 
5349  /* Get the name of the callback user function for raster values */
5350  if (PG_ARGISNULL(3)) {
5351 
5353  PG_FREE_IF_COPY(pgraster, 0);
5354  rt_raster_destroy(newrast);
5355 
5356  elog(ERROR, "RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
5357  PG_RETURN_NULL();
5358  }
5359 
5360  oid = PG_GETARG_OID(3);
5361  if (oid == InvalidOid) {
5362 
5364  PG_FREE_IF_COPY(pgraster, 0);
5365  rt_raster_destroy(newrast);
5366 
5367  elog(ERROR, "RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
5368  PG_RETURN_NULL();
5369  }
5370 
5371  fmgr_info(oid, &cbinfo);
5372 
5373  /* function cannot return set */
5374  if (cbinfo.fn_retset) {
5375 
5377  PG_FREE_IF_COPY(pgraster, 0);
5378  rt_raster_destroy(newrast);
5379 
5380  elog(ERROR, "RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
5381  PG_RETURN_NULL();
5382  }
5383  /* function should have correct # of args */
5384  else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
5385 
5387  PG_FREE_IF_COPY(pgraster, 0);
5388  rt_raster_destroy(newrast);
5389 
5390  elog(ERROR, "RASTER_mapAlgebraFct: Function does not have two or three input parameters");
5391  PG_RETURN_NULL();
5392  }
5393 
5394  if (cbinfo.fn_nargs == 2)
5395  k = 1;
5396  else
5397  k = 2;
5398 
5399  if (func_volatile(oid) == 'v') {
5400  elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5401  }
5402 
5403  /* prep function call data */
5404 #if POSTGIS_PGSQL_VERSION < 120
5405  InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5406 
5407  cbdata.argnull[0] = FALSE;
5408  cbdata.argnull[1] = FALSE;
5409  cbdata.argnull[2] = FALSE;
5410 #else
5411  InitFunctionCallInfoData(*cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5412 
5413  cbdata->args[0].isnull = FALSE;
5414  cbdata->args[1].isnull = FALSE;
5415  cbdata->args[2].isnull = FALSE;
5416 #endif
5417 
5418  /* check that the function isn't strict if the args are null. */
5419  if (PG_ARGISNULL(4)) {
5420  if (cbinfo.fn_strict) {
5421 
5423  PG_FREE_IF_COPY(pgraster, 0);
5424  rt_raster_destroy(newrast);
5425 
5426  elog(ERROR, "RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
5427  PG_RETURN_NULL();
5428  }
5429 
5430 #if POSTGIS_PGSQL_VERSION < 120
5431  cbdata.arg[k] = (Datum)NULL;
5432  cbdata.argnull[k] = TRUE;
5433 #else
5434  cbdata->args[k].value = (Datum)NULL;
5435  cbdata->args[k].isnull = TRUE;
5436 #endif
5437  }
5438  else {
5439 #if POSTGIS_PGSQL_VERSION < 120
5440  cbdata.arg[k] = PG_GETARG_DATUM(4);
5441 #else
5442  cbdata->args[k].value = PG_GETARG_DATUM(4);
5443 #endif
5444  }
5445 
5452 
5453  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Band is a nodata band, returning "
5454  "a raster filled with nodata");
5455 
5456  ret = rt_raster_generate_new_band(newrast, newpixeltype,
5457  newinitialvalue, TRUE, newnodatavalue, 0);
5458 
5460  PG_FREE_IF_COPY(pgraster, 0);
5461 
5462  /* Serialize created raster */
5463  pgrtn = rt_raster_serialize(newrast);
5464  rt_raster_destroy(newrast);
5465  if (NULL == pgrtn) {
5466  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5467  PG_RETURN_NULL();
5468  }
5469 
5470  SET_VARSIZE(pgrtn, pgrtn->size);
5471  PG_RETURN_POINTER(pgrtn);
5472  }
5473 
5474 
5479  ret = rt_raster_generate_new_band(newrast, newpixeltype,
5480  newinitialvalue, TRUE, newnodatavalue, 0);
5481 
5482  /* Get the new raster band */
5483  newband = rt_raster_get_band(newrast, 0);
5484  if ( NULL == newband ) {
5485  elog(NOTICE, "Could not modify band for new raster. Returning new "
5486  "raster with the original band");
5487 
5489  PG_FREE_IF_COPY(pgraster, 0);
5490 
5491  /* Serialize created raster */
5492  pgrtn = rt_raster_serialize(newrast);
5493  rt_raster_destroy(newrast);
5494  if (NULL == pgrtn) {
5495  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5496  PG_RETURN_NULL();
5497  }
5498 
5499  SET_VARSIZE(pgrtn, pgrtn->size);
5500  PG_RETURN_POINTER(pgrtn);
5501  }
5502 
5503  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Main computing loop (%d x %d)",
5504  width, height);
5505 
5506  for (x = 0; x < width; x++) {
5507  for(y = 0; y < height; y++) {
5508  ret = rt_band_get_pixel(band, x, y, &r, NULL);
5509 
5514  if (ret == ES_NONE) {
5515  if (FLT_EQ(r, newnodatavalue)) {
5516  if (cbinfo.fn_strict) {
5517  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
5518  continue;
5519  }
5520 #if POSTGIS_PGSQL_VERSION < 120
5521  cbdata.argnull[0] = TRUE;
5522  cbdata.arg[0] = (Datum)NULL;
5523 #else
5524  cbdata->args[0].isnull = TRUE;
5525  cbdata->args[0].value = (Datum)NULL;
5526 #endif
5527  }
5528  else {
5529 #if POSTGIS_PGSQL_VERSION < 120
5530  cbdata.argnull[0] = FALSE;
5531  cbdata.arg[0] = Float8GetDatum(r);
5532 #else
5533  cbdata->args[0].isnull = FALSE;
5534  cbdata->args[0].value = Float8GetDatum(r);
5535 #endif
5536  }
5537 
5538  /* Add pixel positions if callback has proper # of args */
5539  if (cbinfo.fn_nargs == 3) {
5540  Datum d[2];
5541  ArrayType *a;
5542 
5543  d[0] = Int32GetDatum(x+1);
5544  d[1] = Int32GetDatum(y+1);
5545 
5546  a = construct_array(d, 2, INT4OID, sizeof(int32), true, 'i');
5547 
5548 #if POSTGIS_PGSQL_VERSION < 120
5549  cbdata.argnull[1] = FALSE;
5550  cbdata.arg[1] = PointerGetDatum(a);
5551 #else
5552  cbdata->args[1].isnull = FALSE;
5553  cbdata->args[1].value = PointerGetDatum(a);
5554 #endif
5555  }
5556 
5557  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: (%dx%d), r = %f",
5558  x, y, r);
5559 
5560 #if POSTGIS_PGSQL_VERSION < 120
5561  tmpnewval = FunctionCallInvoke(&cbdata);
5562 
5563  if (cbdata.isnull) {
5564  newval = newnodatavalue;
5565  }
5566 #else
5567  tmpnewval = FunctionCallInvoke(cbdata);
5568 
5569  if (cbdata->isnull)
5570  {
5571  newval = newnodatavalue;
5572  }
5573 #endif
5574  else {
5575  newval = DatumGetFloat8(tmpnewval);
5576  }
5577 
5578  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: new value = %f",
5579  newval);
5580 
5581  rt_band_set_pixel(newband, x, y, newval, NULL);
5582  }
5583 
5584  }
5585  }
5586 
5587  /* The newrast band has been modified */
5588 
5589  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster modified, serializing it.");
5590  /* Serialize created raster */
5591 
5593  PG_FREE_IF_COPY(pgraster, 0);
5594 
5595  pgrtn = rt_raster_serialize(newrast);
5596  rt_raster_destroy(newrast);
5597  if (NULL == pgrtn)
5598  PG_RETURN_NULL();
5599 
5600  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster serialized");
5601 
5602  POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFct: returning raster");
5603 
5604  SET_VARSIZE(pgrtn, pgrtn->size);
5605  PG_RETURN_POINTER(pgrtn);
5606 }
char * r
Definition: cu_in_wkt.c:24
#define TRUE
Definition: dbfopen.c:73
#define FALSE
Definition: dbfopen.c:72
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:356
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_raster.c:181
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition: rt_raster.c:213
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_raster.c:485
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:137
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:674
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
Definition: rt_pixel.c:80
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1221
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
rt_pixtype
Definition: librtcore.h:185
@ PT_END
Definition: librtcore.h:197
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_raster.c:168
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:714
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
Definition: rt_raster.c:1342
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_serialize.c:521
#define FLT_EQ(x, y)
Definition: librtcore.h:2235
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
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:1745
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
Definition: rt_band.c:974
@ ES_NONE
Definition: librtcore.h:180
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_raster.c:363
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1730
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:631
const char * rt_pixtype_name(rt_pixtype pixtype)
Definition: rt_pixel.c:110
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition: rt_raster.c:190
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:199
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:725
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition: rt_raster.c:222
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1329
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
band
Definition: ovdump.py:58
nband
Definition: pixval.py:53
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_DEBUG(level, msg)
Definition: rtpostgis.h:61
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:65
unsigned int int32
Definition: shpopen.c:54
Struct definitions.
Definition: librtcore.h:2251

References ovdump::band, ES_NONE, FALSE, FLT_EQ, 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, pixval::x, and pixval::y.

Here is the call graph for this function: