PostGIS  3.0.0dev-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 5071 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().

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