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
Optimization: If a nodataval is provided, use it for newinitialvalue. Then, we can initialize the raster with this value and skip the computation of nodata values one by one in the main computing loop
Optimization: If the raster is only filled with nodata values return right now a raster filled with the newinitialvalue TODO: Call rt_band_check_isnodata instead?
Optimization: If expression resume to 'RAST' and hasnodataval is zero, we can just return the band from the original raster
Create the raster receiving all the computed values. Initialize it to the new initial value
Optimization: If expression is NULL, or all the pixels could be set in one step, return the initialized raster now
We compute a value only for the withdata value pixel since the nodata value has already been set by the first optimization
4488 int x,
y,
nband, width, height;
4490 double newnodatavalue = 0.0;
4491 double newinitialvalue = 0.0;
4492 double newval = 0.0;
4493 char *newexpr = NULL;
4494 char *initexpr = NULL;
4495 char *expression = NULL;
4496 int hasnodataval = 0;
4497 double nodataval = 0.;
4499 int skipcomputation = 0;
4501 const int argkwcount = 3;
4502 enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4503 char *argkw[] = {
"[rast]",
"[rast.x]",
"[rast.y]"};
4504 Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4506 Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4507 uint8_t argpos[3] = {0};
4512 SPIPlanPtr spi_plan = NULL;
4513 SPITupleTable * tuptable = NULL;
4515 char * strFromText = NULL;
4516 Datum *values = NULL;
4517 Datum datum = (Datum)NULL;
4519 bool isnull =
FALSE;
4526 if (PG_ARGISNULL(0)) {
4527 elog(NOTICE,
"Raster is NULL. Returning NULL");
4533 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4536 PG_FREE_IF_COPY(pgraster, 0);
4537 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not deserialize raster");
4543 if (PG_ARGISNULL(1))
4546 nband = PG_GETARG_INT32(1);
4552 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new empty raster...");
4563 if ( NULL == newrast ) {
4564 PG_FREE_IF_COPY(pgraster, 0);
4565 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not create a new raster");
4590 elog(NOTICE,
"Raster is empty. Returning an empty raster");
4592 PG_FREE_IF_COPY(pgraster, 0);
4596 if (NULL == pgrtn) {
4598 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4602 SET_VARSIZE(pgrtn, pgrtn->
size);
4603 PG_RETURN_POINTER(pgrtn);
4614 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
4617 PG_FREE_IF_COPY(pgraster, 0);
4621 if (NULL == pgrtn) {
4622 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4626 SET_VARSIZE(pgrtn, pgrtn->
size);
4627 PG_RETURN_POINTER(pgrtn);
4632 if ( NULL ==
band ) {
4633 elog(NOTICE,
"Could not get the required band. Returning a raster "
4636 PG_FREE_IF_COPY(pgraster, 0);
4640 if (NULL == pgrtn) {
4641 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4645 SET_VARSIZE(pgrtn, pgrtn->
size);
4646 PG_RETURN_POINTER(pgrtn);
4652 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4670 newinitialvalue = newnodatavalue;
4677 if (PG_ARGISNULL(2)) {
4682 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
4685 if (newpixeltype ==
PT_END)
4689 if (newpixeltype ==
PT_END) {
4690 PG_FREE_IF_COPY(pgraster, 0);
4691 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid pixeltype");
4700 if (!PG_ARGISNULL(3)) {
4701 expression = text_to_cstring(PG_GETARG_TEXT_P(3));
4702 len = strlen(
"SELECT (") + strlen(expression) + strlen(
")::double precision");
4703 initexpr = (
char *)palloc(len + 1);
4705 memcpy(initexpr,
"SELECT (", strlen(
"SELECT ("));
4706 memcpy(initexpr + strlen(
"SELECT ("), expression, strlen(expression));
4707 memcpy(initexpr + strlen(
"SELECT (") + strlen(expression),
")::double precision", strlen(
")::double precision"));
4708 initexpr[len] =
'\0';
4726 if (!PG_ARGISNULL(4)) {
4728 nodataval = PG_GETARG_FLOAT8(4);
4729 newinitialvalue = nodataval;
4746 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4747 "a raster filled with nodata");
4750 newinitialvalue,
TRUE, newnodatavalue, 0);
4756 PG_FREE_IF_COPY(pgraster, 0);
4761 if (NULL == pgrtn) {
4762 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4766 SET_VARSIZE(pgrtn, pgrtn->
size);
4767 PG_RETURN_POINTER(pgrtn);
4775 if (initexpr != NULL && ( !strcmp(initexpr,
"SELECT [rast]") || !strcmp(initexpr,
"SELECT [rast.val]") ) && !hasnodataval) {
4778 "Returning raster with band %d from original raster",
nband);
4790 PG_FREE_IF_COPY(pgraster, 0);
4795 if (NULL == pgrtn) {
4796 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4800 SET_VARSIZE(pgrtn, pgrtn->
size);
4801 PG_RETURN_POINTER(pgrtn);
4808 if (initexpr != NULL && strstr(initexpr,
"[rast") == NULL) {
4809 ret = SPI_connect();
4810 if (ret != SPI_OK_CONNECT) {
4811 PG_FREE_IF_COPY(pgraster, 0);
4812 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4817 ret = SPI_execute(initexpr,
FALSE, 0);
4819 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4823 SPI_freetuptable(tuptable);
4824 PG_FREE_IF_COPY(pgraster, 0);
4827 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid construction for expression");
4831 tupdesc = SPI_tuptable->tupdesc;
4832 tuptable = SPI_tuptable;
4834 tuple = tuptable->vals[0];
4835 newexpr = SPI_getvalue(tuple, tupdesc, 1);
4837 POSTGIS_RT_DEBUG(3,
"Constant expression evaluated to NULL, keeping initvalue");
4838 newval = newinitialvalue;
4840 newval = atof(newexpr);
4843 SPI_freetuptable(tuptable);
4850 skipcomputation = 1;
4856 if (!hasnodataval) {
4857 newinitialvalue = newval;
4858 skipcomputation = 2;
4862 else if (
FLT_NEQ(newval, newinitialvalue)) {
4863 skipcomputation = 2;
4872 newinitialvalue,
TRUE, newnodatavalue, 0);
4879 if (expression == NULL || skipcomputation == 2) {
4885 PG_FREE_IF_COPY(pgraster, 0);
4890 if (NULL == pgrtn) {
4891 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4895 SET_VARSIZE(pgrtn, pgrtn->
size);
4896 PG_RETURN_POINTER(pgrtn);
4899 RASTER_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new raster band...");
4903 if ( NULL == newband ) {
4904 elog(NOTICE,
"Could not modify band for new raster. Returning new "
4905 "raster with the original band");
4910 PG_FREE_IF_COPY(pgraster, 0);
4915 if (NULL == pgrtn) {
4916 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4920 SET_VARSIZE(pgrtn, pgrtn->
size);
4921 PG_RETURN_POINTER(pgrtn);
4927 if (initexpr != NULL) {
4930 pfree(initexpr); initexpr=newexpr;
4932 sprintf(place,
"$1");
4933 for (i = 0, j = 1; i < argkwcount; i++) {
4936 pfree(initexpr); initexpr=newexpr;
4938 argtype[argcount] = argkwtypes[i];
4942 sprintf(place,
"$%d", j);
4952 values = (Datum *) palloc(
sizeof(Datum) * argcount);
4953 if (values == NULL) {
4958 PG_FREE_IF_COPY(pgraster, 0);
4961 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
4966 nulls = (
char *)palloc(argcount);
4967 if (nulls == NULL) {
4972 PG_FREE_IF_COPY(pgraster, 0);
4975 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
4980 ret = SPI_connect();
4981 if (ret != SPI_OK_CONNECT) {
4986 PG_FREE_IF_COPY(pgraster, 0);
4989 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4994 spi_plan = SPI_prepare(initexpr, argcount, argtype);
4996 if (spi_plan == NULL) {
4999 PG_FREE_IF_COPY(pgraster, 0);
5006 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not prepare expression");
5011 for (
x = 0;
x < width;
x++) {
5012 for(
y = 0;
y < height;
y++) {
5020 if (skipcomputation == 0) {
5021 if (initexpr != NULL) {
5023 memset(nulls,
'n', argcount);
5025 for (i = 0; i < argkwcount; i++) {
5027 if (idx < 1)
continue;
5032 values[idx] = Int32GetDatum(
x+1);
5037 values[idx] = Int32GetDatum(
y+1);
5040 else if (i == kVAL ) {
5041 values[idx] = Float8GetDatum(
r);
5047 ret = SPI_execute_plan(spi_plan, values, nulls,
FALSE, 0);
5048 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
5049 SPI_processed != 1) {
5051 SPI_freetuptable(tuptable);
5053 SPI_freeplan(spi_plan);
5061 PG_FREE_IF_COPY(pgraster, 0);
5064 elog(ERROR,
"RASTER_mapAlgebraExpr: Error executing prepared plan");
5069 tupdesc = SPI_tuptable->tupdesc;
5070 tuptable = SPI_tuptable;
5072 tuple = tuptable->vals[0];
5073 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5074 if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5075 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) errored, skip setting",
x+1,
y+1,
r);
5076 newval = newinitialvalue;
5078 else if ( isnull ) {
5079 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting",
x+1,
y+1,
r);
5080 newval = newinitialvalue;
5082 newval = DatumGetFloat8(datum);
5085 SPI_freetuptable(tuptable);
5089 newval = newinitialvalue;
5102 if (initexpr != NULL) {
5103 SPI_freeplan(spi_plan);
5117 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: raster modified, serializing it.");
5121 PG_FREE_IF_COPY(pgraster, 0);
5128 SET_VARSIZE(pgrtn, pgrtn->
size);
5136 PG_RETURN_POINTER(pgrtn);
#define RASTER_DEBUG(level, msg)
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
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.
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
uint16_t rt_raster_get_num_bands(rt_raster raster)
uint16_t rt_raster_get_height(rt_raster raster)
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
const char * rt_pixtype_name(rt_pixtype pixtype)
uint16_t rt_raster_get_width(rt_raster raster)
int rt_raster_copy_band(rt_raster torast, rt_raster fromrast, int fromindex, int toindex)
Copy one band from one raster to another.
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
char * rtpg_strreplace(const char *str, const char *oldstr, const char *newstr, int *count)
#define POSTGIS_RT_DEBUG(level, msg)
#define POSTGIS_RT_DEBUGF(level, msg,...)