PostGIS  2.5.0dev-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 5074 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.

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