PostGIS  2.5.2dev-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 5073 of file rtpg_mapalgebra.c.

References ovdump::band, ES_NONE, FALSE, FLT_EQ, rt_raster_serialized_t::height, pixval::nband, PG_FUNCTION_INFO_V1(), POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, PT_END, r, rtrowdump::raster, RASTER_mapAlgebraFctNgb(), 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, text_to_cstring(), TRUE, rt_raster_serialized_t::width, pixval::x, and pixval::y.

Referenced by RASTER_mapAlgebraExpr().

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