RASTER_mapAlgebraExpr()

Datum RASTER_mapAlgebraExpr

(

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 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

Optimization: If expression resume to a constant (it does not contain [rast)

Compute the new value, set it and we will return after creating the new 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

Definition at line 4577 of file rtpg_mapalgebra.c.

{
    rt_pgraster *pgraster = NULL;
    rt_pgraster *pgrtn = NULL;
    rt_raster raster = NULL;
    rt_raster newrast = NULL;
    rt_band band = NULL;
    rt_band newband = NULL;
    int x, y, nband, width, height;
    double r;
    double newnodatavalue = 0.0;
    double newinitialvalue = 0.0;
    double newval = 0.0;
    char *newexpr = NULL;
    char *initexpr = NULL;
    char *expression = NULL;
    int hasnodataval = 0;
    double nodataval = 0.;
    rt_pixtype newpixeltype;
    int skipcomputation = 0;
    int len = 0;
    const int argkwcount = 3;
    enum KEYWORDS { kVAL=0, kX=1, kY=2 };
    char *argkw[] = {"[rast]", "[rast.x]", "[rast.y]"};
    Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
    int argcount = 0;
    Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
    uint8_t argpos[3] = {0};
    char place[12];
    int idx = 0;
    int ret = -1;
    TupleDesc tupdesc;
    SPIPlanPtr spi_plan = NULL;
    SPITupleTable * tuptable = NULL;
    HeapTuple tuple;
    char * strFromText = NULL;
    Datum *values = NULL;
    Datum datum = (Datum)NULL;
    char *nulls = NULL;
    bool isnull = FALSE;
    int i = 0;
    int j = 0;
 
    POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraExpr: Starting...");
 
    /* Check raster */
    if (PG_ARGISNULL(0)) {
        elog(NOTICE, "Raster is NULL. Returning NULL");
        PG_RETURN_NULL();
    }
 
 
    /* Deserialize raster */
    pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
    raster = rt_raster_deserialize(pgraster, FALSE);
    if (NULL == raster) {
        PG_FREE_IF_COPY(pgraster, 0);
        elog(ERROR, "RASTER_mapAlgebraExpr: Could not deserialize raster");
        PG_RETURN_NULL();
    }
 
    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting arguments...");
 
    if (PG_ARGISNULL(1))
        nband = 1;
    else
        nband = PG_GETARG_INT32(1);
 
    if (nband < 1)
        nband = 1;
 
 
    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new empty raster...");
 
    width = rt_raster_get_width(raster);
    height = rt_raster_get_height(raster);
 
    newrast = rt_raster_new(width, height);
 
    if ( NULL == newrast ) {
        PG_FREE_IF_COPY(pgraster, 0);
        elog(ERROR, "RASTER_mapAlgebraExpr: Could not create a new raster");
        PG_RETURN_NULL();
    }
 
    rt_raster_set_scale(newrast,
            rt_raster_get_x_scale(raster),
            rt_raster_get_y_scale(raster));
 
    rt_raster_set_offsets(newrast,
            rt_raster_get_x_offset(raster),
            rt_raster_get_y_offset(raster));
 
    rt_raster_set_skews(newrast,
            rt_raster_get_x_skew(raster),
            rt_raster_get_y_skew(raster));
 
    rt_raster_set_srid(newrast, rt_raster_get_srid(raster));
 
 
    if (rt_raster_is_empty(newrast))
    {
        elog(NOTICE, "Raster is empty. Returning an empty raster");
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
 
        pgrtn = rt_raster_serialize(newrast);
        rt_raster_destroy(newrast);
        if (NULL == pgrtn) {
 
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
            PG_RETURN_NULL();
        }
 
        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
    }
 
 
    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Getting raster band %d...", nband);
 
    if (!rt_raster_has_band(raster, nband - 1)) {
        elog(NOTICE, "Raster does not have the required band. Returning a raster "
                "without a band");
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
 
        pgrtn = rt_raster_serialize(newrast);
        rt_raster_destroy(newrast);
        if (NULL == pgrtn) {
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
            PG_RETURN_NULL();
        }
 
        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
    }
 
    /* Get the raster band */
    band = rt_raster_get_band(raster, nband - 1);
    if ( NULL == band ) {
        elog(NOTICE, "Could not get the required band. Returning a raster "
                "without a band");
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
 
        pgrtn = rt_raster_serialize(newrast);
        rt_raster_destroy(newrast);
        if (NULL == pgrtn) {
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
            PG_RETURN_NULL();
        }
 
        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
    }
 
     /*
     * Get NODATA value
     */
    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Getting NODATA value for band...");
 
    if (rt_band_get_hasnodata_flag(band)) {
        rt_band_get_nodata(band, &newnodatavalue);
    }
 
    else {
        newnodatavalue = rt_band_get_min_value(band);
    }
 
    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: NODATA value for band: = %f",
        newnodatavalue);
 
    newinitialvalue = newnodatavalue;
 
    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Setting pixeltype...");
 
    if (PG_ARGISNULL(2)) {
        newpixeltype = rt_band_get_pixtype(band);
    }
 
    else {
        strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
        newpixeltype = rt_pixtype_index_from_name(strFromText);
        pfree(strFromText);
        if (newpixeltype == PT_END)
            newpixeltype = rt_band_get_pixtype(band);
    }
 
    if (newpixeltype == PT_END) {
        PG_FREE_IF_COPY(pgraster, 0);
        elog(ERROR, "RASTER_mapAlgebraExpr: Invalid pixeltype");
        PG_RETURN_NULL();
    }
 
    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Pixeltype set to %s",
        rt_pixtype_name(newpixeltype));
 
 
    /* Construct expression for raster values */
    if (!PG_ARGISNULL(3)) {
        expression = text_to_cstring(PG_GETARG_TEXT_P(3));
        len = strlen("SELECT (") + strlen(expression) + strlen(")::double precision");
        initexpr = (char *)palloc(len + 1);
 
        memcpy(initexpr, "SELECT (", strlen("SELECT ("));
        memcpy(initexpr + strlen("SELECT ("), expression, strlen(expression));
        memcpy(initexpr + strlen("SELECT (") + strlen(expression), ")::double precision", strlen(")::double precision"));
        initexpr[len] = '\0';
 
        POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression is %s", initexpr);
 
        /* We don't need this memory */
        /*
                                pfree(expression);
        expression = NULL;
                                */
    }
 
 
 
    if (!PG_ARGISNULL(4)) {
                                hasnodataval = 1;
                                nodataval = PG_GETARG_FLOAT8(4);
                                newinitialvalue = nodataval;
 
        POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new initial value = %f",
            newinitialvalue);
    }
    else
        hasnodataval = 0;
 
 
 
    if (rt_band_get_isnodata_flag(band)) {
 
        POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: Band is a nodata band, returning "
                "a raster filled with nodata");
 
        ret = rt_raster_generate_new_band(newrast, newpixeltype,
                newinitialvalue, TRUE, newnodatavalue, 0);
 
        /* Free memory */
        if (initexpr)
            pfree(initexpr);
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
 
        /* Serialize created raster */
        pgrtn = rt_raster_serialize(newrast);
        rt_raster_destroy(newrast);
        if (NULL == pgrtn) {
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
            PG_RETURN_NULL();
        }
 
        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
    }
 
 
    if (initexpr != NULL && ( !strcmp(initexpr, "SELECT [rast]") || !strcmp(initexpr, "SELECT [rast.val]") ) && !hasnodataval) {
 
        POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Expression resumes to RAST. "
                "Returning raster with band %d from original raster", nband);
 
        POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster has %d bands",
                rt_raster_get_num_bands(newrast));
 
        rt_raster_copy_band(newrast, raster, nband - 1, 0);
 
        POSTGIS_RT_DEBUGF(4, "RASTER_mapAlgebraExpr: New raster now has %d bands",
                rt_raster_get_num_bands(newrast));
 
        pfree(initexpr);
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
 
        /* Serialize created raster */
        pgrtn = rt_raster_serialize(newrast);
        rt_raster_destroy(newrast);
        if (NULL == pgrtn) {
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
            PG_RETURN_NULL();
        }
 
        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
    }
 
    if (initexpr != NULL && strstr(initexpr, "[rast") == NULL) {
        ret = SPI_connect();
        if (ret != SPI_OK_CONNECT) {
            PG_FREE_IF_COPY(pgraster, 0);
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
            PG_RETURN_NULL();
        };
 
        /* Execute the expression into newval */
        ret = SPI_execute(initexpr, FALSE, 0);
 
        if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
 
            /* Free memory allocated out of the current context */
            if (SPI_tuptable)
                SPI_freetuptable(tuptable);
            PG_FREE_IF_COPY(pgraster, 0);
 
            SPI_finish();
            elog(ERROR, "RASTER_mapAlgebraExpr: Invalid construction for expression");
            PG_RETURN_NULL();
        }
 
        tupdesc = SPI_tuptable->tupdesc;
        tuptable = SPI_tuptable;
 
        tuple = tuptable->vals[0];
        newexpr = SPI_getvalue(tuple, tupdesc, 1);
        if ( ! newexpr ) {
            POSTGIS_RT_DEBUG(3, "Constant expression evaluated to NULL, keeping initvalue");
            newval = newinitialvalue;
        } else {
            newval = atof(newexpr);
        }
 
        SPI_freetuptable(tuptable);
 
        POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: New raster value = %f",
                newval);
 
        SPI_finish();
 
        skipcomputation = 1;
 
        if (!hasnodataval) {
            newinitialvalue = newval;
            skipcomputation = 2;
        }
 
        /* Return the new raster as it will be before computing pixel by pixel */
        else if (FLT_NEQ(newval, newinitialvalue)) {
            skipcomputation = 2;
        }
    }
 
    ret = rt_raster_generate_new_band(newrast, newpixeltype,
            newinitialvalue, TRUE, newnodatavalue, 0);
 
    /*if (initexpr == NULL || skipcomputation == 2) {*/
    if (expression == NULL || skipcomputation == 2) {
 
        /* Free memory */
        if (initexpr)
            pfree(initexpr);
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
 
        /* Serialize created raster */
        pgrtn = rt_raster_serialize(newrast);
        rt_raster_destroy(newrast);
        if (NULL == pgrtn) {
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
            PG_RETURN_NULL();
        }
 
        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
    }
 
    RASTER_DEBUG(3, "RASTER_mapAlgebraExpr: Creating new raster band...");
 
    /* Get the new raster band */
    newband = rt_raster_get_band(newrast, 0);
    if ( NULL == newband ) {
        elog(NOTICE, "Could not modify band for new raster. Returning new "
                "raster with the original band");
 
        if (initexpr)
            pfree(initexpr);
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
 
        /* Serialize created raster */
        pgrtn = rt_raster_serialize(newrast);
        rt_raster_destroy(newrast);
        if (NULL == pgrtn) {
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not serialize raster");
            PG_RETURN_NULL();
        }
 
        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
    }
 
    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: Main computing loop (%d x %d)",
            width, height);
 
    if (initexpr != NULL) {
        /* Convert [rast.val] to [rast] */
        newexpr = rtpg_strreplace(initexpr, "[rast.val]", "[rast]", NULL);
        pfree(initexpr); initexpr=newexpr;
 
        sprintf(place,"$1");
        for (i = 0, j = 1; i < argkwcount; i++) {
            len = 0;
            newexpr = rtpg_strreplace(initexpr, argkw[i], place, &len);
            pfree(initexpr); initexpr=newexpr;
            if (len > 0) {
                argtype[argcount] = argkwtypes[i];
                argcount++;
                argpos[i] = j++;
 
                sprintf(place, "$%d", j);
            }
            else {
                argpos[i] = 0;
            }
        }
 
        POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: initexpr = %s", initexpr);
 
        /* define values */
        values = (Datum *) palloc(sizeof(Datum) * argcount);
        if (values == NULL) {
 
            SPI_finish();
 
            rt_raster_destroy(raster);
            PG_FREE_IF_COPY(pgraster, 0);
            rt_raster_destroy(newrast);
 
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
            PG_RETURN_NULL();
        }
 
        /* define nulls */
        nulls = (char *)palloc(argcount);
        if (nulls == NULL) {
 
            SPI_finish();
 
            rt_raster_destroy(raster);
            PG_FREE_IF_COPY(pgraster, 0);
            rt_raster_destroy(newrast);
 
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
            PG_RETURN_NULL();
        }
 
        /* Connect to SPI and prepare the expression */
        ret = SPI_connect();
        if (ret != SPI_OK_CONNECT) {
 
            if (initexpr)
                pfree(initexpr);
            rt_raster_destroy(raster);
            PG_FREE_IF_COPY(pgraster, 0);
            rt_raster_destroy(newrast);
 
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
            PG_RETURN_NULL();
        };
 
        /* Type of all arguments is FLOAT8OID */
        spi_plan = SPI_prepare(initexpr, argcount, argtype);
 
        if (spi_plan == NULL) {
 
            rt_raster_destroy(raster);
            PG_FREE_IF_COPY(pgraster, 0);
            rt_raster_destroy(newrast);
 
            SPI_finish();
 
            pfree(initexpr);
 
            elog(ERROR, "RASTER_mapAlgebraExpr: Could not prepare expression");
            PG_RETURN_NULL();
        }
    }
 
    for (x = 0; x < width; x++) {
        for(y = 0; y < height; y++) {
            ret = rt_band_get_pixel(band, x, y, &r, NULL);
 
            if (ret == ES_NONE && FLT_NEQ(r, newnodatavalue)) {
                if (skipcomputation == 0) {
                    if (initexpr != NULL) {
                        /* Reset the null arg flags. */
                        memset(nulls, 'n', argcount);
 
                        for (i = 0; i < argkwcount; i++) {
                            idx = argpos[i];
                            if (idx < 1) continue;
                            idx--;
 
                            if (i == kX ) {
                                /* x is 0 based index, but SQL expects 1 based index */
                                values[idx] = Int32GetDatum(x+1);
                                nulls[idx] = ' ';
                            }
                            else if (i == kY) {
                                /* y is 0 based index, but SQL expects 1 based index */
                                values[idx] = Int32GetDatum(y+1);
                                nulls[idx] = ' ';
                            }
                            else if (i == kVAL ) {
                                values[idx] = Float8GetDatum(r);
                                nulls[idx] = ' ';
                            }
 
                        }
 
                        ret = SPI_execute_plan(spi_plan, values, nulls, FALSE, 0);
                        if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
                                SPI_processed != 1) {
                            if (SPI_tuptable)
                                SPI_freetuptable(tuptable);
 
                            SPI_freeplan(spi_plan);
                            SPI_finish();
 
                            pfree(values);
                            pfree(nulls);
                            pfree(initexpr);
 
                            rt_raster_destroy(raster);
                            PG_FREE_IF_COPY(pgraster, 0);
                            rt_raster_destroy(newrast);
 
                            elog(ERROR, "RASTER_mapAlgebraExpr: Error executing prepared plan");
 
                            PG_RETURN_NULL();
                        }
 
                        tupdesc = SPI_tuptable->tupdesc;
                        tuptable = SPI_tuptable;
 
                        tuple = tuptable->vals[0];
                        datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
                        if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
                            POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) errored, skip setting", x+1,y+1,r);
                            newval = newinitialvalue;
                        }
                        else if ( isnull ) {
                            POSTGIS_RT_DEBUGF(3, "Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting", x+1,y+1,r);
                            newval = newinitialvalue;
                        } else {
                            newval = DatumGetFloat8(datum);
                        }
 
                        SPI_freetuptable(tuptable);
                    }
 
                    else
                        newval = newinitialvalue;
 
                    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraExpr: new value = %f",
                        newval);
                }
 
 
                rt_band_set_pixel(newband, x, y, newval, NULL);
            }
 
        }
    }
 
    if (initexpr != NULL) {
        SPI_freeplan(spi_plan);
        SPI_finish();
 
        pfree(values);
        pfree(nulls);
        pfree(initexpr);
    }
    else {
        POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: no SPI cleanup");
    }
 
 
    /* The newrast band has been modified */
 
    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster modified, serializing it.");
    /* Serialize created raster */
 
    rt_raster_destroy(raster);
    PG_FREE_IF_COPY(pgraster, 0);
 
    pgrtn = rt_raster_serialize(newrast);
    rt_raster_destroy(newrast);
    if (NULL == pgrtn)
        PG_RETURN_NULL();
 
    SET_VARSIZE(pgrtn, pgrtn->size);
 
    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraExpr: raster serialized");
 
 
    POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraExpr: returning raster");
 
 
    PG_RETURN_POINTER(pgrtn);
}

References ovdump::band, ES_NONE, FALSE, FLT_NEQ, pixval::nband, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, PT_END, r, rtrowdump::raster, RASTER_DEBUG, 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_copy_band(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_generate_new_band(), rt_raster_get_band(), rt_raster_get_height(), rt_raster_get_num_bands(), 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(), rtpg_strreplace(), rt_raster_serialized_t::size, TRUE, pixval::x, and pixval::y.

Here is the call graph for this function: