PostGIS  2.3.7dev-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 5076 of file rtpg_mapalgebra.c.

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.

5077 {
5078  rt_pgraster *pgraster = NULL;
5079  rt_pgraster *pgrtn = NULL;
5080  rt_raster raster = NULL;
5081  rt_raster newrast = NULL;
5082  rt_band band = NULL;
5083  rt_band newband = NULL;
5084  int x, y, nband, width, height;
5085  double r;
5086  double newnodatavalue = 0.0;
5087  double newinitialvalue = 0.0;
5088  double newval = 0.0;
5089  rt_pixtype newpixeltype;
5090  int ret = -1;
5091  Oid oid;
5092  FmgrInfo cbinfo;
5093  FunctionCallInfoData cbdata;
5094  Datum tmpnewval;
5095  char * strFromText = NULL;
5096  int k = 0;
5097 
5098  POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFct: STARTING...");
5099 
5100  /* Check raster */
5101  if (PG_ARGISNULL(0)) {
5102  elog(WARNING, "Raster is NULL. Returning NULL");
5103  PG_RETURN_NULL();
5104  }
5105 
5106 
5107  /* Deserialize raster */
5108  pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5109  raster = rt_raster_deserialize(pgraster, FALSE);
5110  if (NULL == raster) {
5111  PG_FREE_IF_COPY(pgraster, 0);
5112  elog(ERROR, "RASTER_mapAlgebraFct: Could not deserialize raster");
5113  PG_RETURN_NULL();
5114  }
5115 
5116  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting arguments...");
5117 
5118  /* Get the rest of the arguments */
5119 
5120  if (PG_ARGISNULL(1))
5121  nband = 1;
5122  else
5123  nband = PG_GETARG_INT32(1);
5124 
5125  if (nband < 1)
5126  nband = 1;
5127 
5128  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Creating new empty raster...");
5129 
5134  width = rt_raster_get_width(raster);
5135  height = rt_raster_get_height(raster);
5136 
5137  newrast = rt_raster_new(width, height);
5138 
5139  if ( NULL == newrast ) {
5140 
5141  rt_raster_destroy(raster);
5142  PG_FREE_IF_COPY(pgraster, 0);
5143 
5144  elog(ERROR, "RASTER_mapAlgebraFct: Could not create a new raster");
5145  PG_RETURN_NULL();
5146  }
5147 
5148  rt_raster_set_scale(newrast,
5149  rt_raster_get_x_scale(raster),
5150  rt_raster_get_y_scale(raster));
5151 
5152  rt_raster_set_offsets(newrast,
5153  rt_raster_get_x_offset(raster),
5154  rt_raster_get_y_offset(raster));
5155 
5156  rt_raster_set_skews(newrast,
5157  rt_raster_get_x_skew(raster),
5158  rt_raster_get_y_skew(raster));
5159 
5160  rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
5161 
5162 
5167  if (rt_raster_is_empty(newrast))
5168  {
5169  elog(NOTICE, "Raster is empty. Returning an empty raster");
5170  rt_raster_destroy(raster);
5171  PG_FREE_IF_COPY(pgraster, 0);
5172 
5173  pgrtn = rt_raster_serialize(newrast);
5174  rt_raster_destroy(newrast);
5175  if (NULL == pgrtn) {
5176  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5177  PG_RETURN_NULL();
5178  }
5179 
5180  SET_VARSIZE(pgrtn, pgrtn->size);
5181  PG_RETURN_POINTER(pgrtn);
5182  }
5183 
5184  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Getting raster band %d...", nband);
5185 
5190  if (!rt_raster_has_band(raster, nband - 1)) {
5191  elog(NOTICE, "Raster does not have the required band. Returning a raster "
5192  "without a band");
5193  rt_raster_destroy(raster);
5194  PG_FREE_IF_COPY(pgraster, 0);
5195 
5196  pgrtn = rt_raster_serialize(newrast);
5197  rt_raster_destroy(newrast);
5198  if (NULL == pgrtn) {
5199  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5200  PG_RETURN_NULL();
5201  }
5202 
5203  SET_VARSIZE(pgrtn, pgrtn->size);
5204  PG_RETURN_POINTER(pgrtn);
5205  }
5206 
5207  /* Get the raster band */
5208  band = rt_raster_get_band(raster, nband - 1);
5209  if ( NULL == band ) {
5210  elog(NOTICE, "Could not get the required band. Returning a raster "
5211  "without a band");
5212  rt_raster_destroy(raster);
5213  PG_FREE_IF_COPY(pgraster, 0);
5214 
5215  pgrtn = rt_raster_serialize(newrast);
5216  rt_raster_destroy(newrast);
5217  if (NULL == pgrtn) {
5218  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5219  PG_RETURN_NULL();
5220  }
5221 
5222  SET_VARSIZE(pgrtn, pgrtn->size);
5223  PG_RETURN_POINTER(pgrtn);
5224  }
5225 
5226  /*
5227  * Get NODATA value
5228  */
5229  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting NODATA value for band...");
5230 
5231  if (rt_band_get_hasnodata_flag(band)) {
5232  rt_band_get_nodata(band, &newnodatavalue);
5233  }
5234 
5235  else {
5236  newnodatavalue = rt_band_get_min_value(band);
5237  }
5238 
5239  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: NODATA value for band: %f",
5240  newnodatavalue);
5246  newinitialvalue = newnodatavalue;
5247 
5251  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Setting pixeltype...");
5252 
5253  if (PG_ARGISNULL(2)) {
5254  newpixeltype = rt_band_get_pixtype(band);
5255  }
5256 
5257  else {
5258  strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5259  newpixeltype = rt_pixtype_index_from_name(strFromText);
5260  pfree(strFromText);
5261  if (newpixeltype == PT_END)
5262  newpixeltype = rt_band_get_pixtype(band);
5263  }
5264 
5265  if (newpixeltype == PT_END) {
5266 
5267  rt_raster_destroy(raster);
5268  PG_FREE_IF_COPY(pgraster, 0);
5269  rt_raster_destroy(newrast);
5270 
5271  elog(ERROR, "RASTER_mapAlgebraFct: Invalid pixeltype");
5272  PG_RETURN_NULL();
5273  }
5274 
5275  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Pixeltype set to %s",
5276  rt_pixtype_name(newpixeltype));
5277 
5278  /* Get the name of the callback user function for raster values */
5279  if (PG_ARGISNULL(3)) {
5280 
5281  rt_raster_destroy(raster);
5282  PG_FREE_IF_COPY(pgraster, 0);
5283  rt_raster_destroy(newrast);
5284 
5285  elog(ERROR, "RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
5286  PG_RETURN_NULL();
5287  }
5288 
5289  oid = PG_GETARG_OID(3);
5290  if (oid == InvalidOid) {
5291 
5292  rt_raster_destroy(raster);
5293  PG_FREE_IF_COPY(pgraster, 0);
5294  rt_raster_destroy(newrast);
5295 
5296  elog(ERROR, "RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
5297  PG_RETURN_NULL();
5298  }
5299 
5300  fmgr_info(oid, &cbinfo);
5301 
5302  /* function cannot return set */
5303  if (cbinfo.fn_retset) {
5304 
5305  rt_raster_destroy(raster);
5306  PG_FREE_IF_COPY(pgraster, 0);
5307  rt_raster_destroy(newrast);
5308 
5309  elog(ERROR, "RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
5310  PG_RETURN_NULL();
5311  }
5312  /* function should have correct # of args */
5313  else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
5314 
5315  rt_raster_destroy(raster);
5316  PG_FREE_IF_COPY(pgraster, 0);
5317  rt_raster_destroy(newrast);
5318 
5319  elog(ERROR, "RASTER_mapAlgebraFct: Function does not have two or three input parameters");
5320  PG_RETURN_NULL();
5321  }
5322 
5323  if (cbinfo.fn_nargs == 2)
5324  k = 1;
5325  else
5326  k = 2;
5327 
5328  if (func_volatile(oid) == 'v') {
5329  elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5330  }
5331 
5332  /* prep function call data */
5333 #if POSTGIS_PGSQL_VERSION <= 90
5334  InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL);
5335 #else
5336  InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5337 #endif
5338  memset(cbdata.argnull, FALSE, sizeof(bool) * cbinfo.fn_nargs);
5339 
5340  /* check that the function isn't strict if the args are null. */
5341  if (PG_ARGISNULL(4)) {
5342  if (cbinfo.fn_strict) {
5343 
5344  rt_raster_destroy(raster);
5345  PG_FREE_IF_COPY(pgraster, 0);
5346  rt_raster_destroy(newrast);
5347 
5348  elog(ERROR, "RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
5349  PG_RETURN_NULL();
5350  }
5351 
5352  cbdata.arg[k] = (Datum)NULL;
5353  cbdata.argnull[k] = TRUE;
5354  }
5355  else {
5356  cbdata.arg[k] = PG_GETARG_DATUM(4);
5357  }
5358 
5364  if (rt_band_get_isnodata_flag(band)) {
5365 
5366  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Band is a nodata band, returning "
5367  "a raster filled with nodata");
5368 
5369  ret = rt_raster_generate_new_band(newrast, newpixeltype,
5370  newinitialvalue, TRUE, newnodatavalue, 0);
5371 
5372  rt_raster_destroy(raster);
5373  PG_FREE_IF_COPY(pgraster, 0);
5374 
5375  /* Serialize created raster */
5376  pgrtn = rt_raster_serialize(newrast);
5377  rt_raster_destroy(newrast);
5378  if (NULL == pgrtn) {
5379  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5380  PG_RETURN_NULL();
5381  }
5382 
5383  SET_VARSIZE(pgrtn, pgrtn->size);
5384  PG_RETURN_POINTER(pgrtn);
5385  }
5386 
5387 
5392  ret = rt_raster_generate_new_band(newrast, newpixeltype,
5393  newinitialvalue, TRUE, newnodatavalue, 0);
5394 
5395  /* Get the new raster band */
5396  newband = rt_raster_get_band(newrast, 0);
5397  if ( NULL == newband ) {
5398  elog(NOTICE, "Could not modify band for new raster. Returning new "
5399  "raster with the original band");
5400 
5401  rt_raster_destroy(raster);
5402  PG_FREE_IF_COPY(pgraster, 0);
5403 
5404  /* Serialize created raster */
5405  pgrtn = rt_raster_serialize(newrast);
5406  rt_raster_destroy(newrast);
5407  if (NULL == pgrtn) {
5408  elog(ERROR, "RASTER_mapAlgebraFct: Could not serialize raster");
5409  PG_RETURN_NULL();
5410  }
5411 
5412  SET_VARSIZE(pgrtn, pgrtn->size);
5413  PG_RETURN_POINTER(pgrtn);
5414  }
5415 
5416  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Main computing loop (%d x %d)",
5417  width, height);
5418 
5419  for (x = 0; x < width; x++) {
5420  for(y = 0; y < height; y++) {
5421  ret = rt_band_get_pixel(band, x, y, &r, NULL);
5422 
5427  if (ret == ES_NONE) {
5428  if (FLT_EQ(r, newnodatavalue)) {
5429  if (cbinfo.fn_strict) {
5430  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
5431  continue;
5432  }
5433  cbdata.argnull[0] = TRUE;
5434  cbdata.arg[0] = (Datum)NULL;
5435  }
5436  else {
5437  cbdata.argnull[0] = FALSE;
5438  cbdata.arg[0] = Float8GetDatum(r);
5439  }
5440 
5441  /* Add pixel positions if callback has proper # of args */
5442  if (cbinfo.fn_nargs == 3) {
5443  Datum d[2];
5444  ArrayType *a;
5445 
5446  d[0] = Int32GetDatum(x+1);
5447  d[1] = Int32GetDatum(y+1);
5448 
5449  a = construct_array(d, 2, INT4OID, sizeof(int32), true, 'i');
5450 
5451  cbdata.argnull[1] = FALSE;
5452  cbdata.arg[1] = PointerGetDatum(a);
5453  }
5454 
5455  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: (%dx%d), r = %f",
5456  x, y, r);
5457 
5458  tmpnewval = FunctionCallInvoke(&cbdata);
5459 
5460  if (cbdata.isnull) {
5461  newval = newnodatavalue;
5462  }
5463  else {
5464  newval = DatumGetFloat8(tmpnewval);
5465  }
5466 
5467  POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: new value = %f",
5468  newval);
5469 
5470  rt_band_set_pixel(newband, x, y, newval, NULL);
5471  }
5472 
5473  }
5474  }
5475 
5476  /* The newrast band has been modified */
5477 
5478  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster modified, serializing it.");
5479  /* Serialize created raster */
5480 
5481  rt_raster_destroy(raster);
5482  PG_FREE_IF_COPY(pgraster, 0);
5483 
5484  pgrtn = rt_raster_serialize(newrast);
5485  rt_raster_destroy(newrast);
5486  if (NULL == pgrtn)
5487  PG_RETURN_NULL();
5488 
5489  POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster serialized");
5490 
5491  POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFct: returning raster");
5492 
5493  SET_VARSIZE(pgrtn, pgrtn->size);
5494  PG_RETURN_POINTER(pgrtn);
5495 }
unsigned int int32
Definition: shpopen.c:273
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_serialize.c:521
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition: rt_raster.c:213
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition: rt_raster.c:190
char * r
Definition: cu_in_wkt.c:24
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
#define FLT_EQ(x, y)
Definition: librtcore.h:2197
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
rt_pixtype
Definition: librtcore.h:197
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1338
#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
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1088
tuple nband
Definition: pixval.py:52
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:137
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
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
Definition: rt_pixel.c:80
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:199
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
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:1351
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_raster.c:363
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:356
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_raster.c:181
tuple x
Definition: pixval.py:53
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
#define FALSE
Definition: dbfopen.c:168
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_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:498
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:581
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:841
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
#define POSTGIS_RT_DEBUG(level, msg)
Definition: rtpostgis.h:53
#define TRUE
Definition: dbfopen.c:169
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:717
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition: rt_raster.c:222
tuple y
Definition: pixval.py:54

Here is the call graph for this function: