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 7427 of file rt_pg.c.

References ovdump::band, ES_NONE, FALSE, FLT_EQ, rt_raster_serialized_t::height, 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, rt_raster_serialized_t::width, pixval::x, and pixval::y.

{
    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;
    rt_pixtype newpixeltype;
    int ret = -1;
    Oid oid;
    FmgrInfo cbinfo;
    FunctionCallInfoData cbdata;
    Datum tmpnewval;
    char * strFromText = NULL;
    int k = 0;

    POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFct: STARTING...");

    /* Check raster */
    if (PG_ARGISNULL(0)) {
        elog(WARNING, "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_mapAlgebraFct: Could not deserialize raster");
                        PG_RETURN_NULL();    
    }

    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Getting arguments...");

    /* Get the rest of the arguments */

    if (PG_ARGISNULL(1))
        nband = 1;
    else
        nband = PG_GETARG_INT32(1);

    if (nband < 1)
        nband = 1;
    
    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: 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 ) {

        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);

        elog(ERROR, "RASTER_mapAlgebraFct: 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_mapAlgebraFct: Could not serialize raster");
            PG_RETURN_NULL();
        }

        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
    }

    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: 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_mapAlgebraFct: 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_mapAlgebraFct: Could not serialize raster");
            PG_RETURN_NULL();
        }

        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
    }

    /*
    * Get NODATA value
    */
    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: 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_mapAlgebraFct: NODATA value for band: %f",
            newnodatavalue);
    newinitialvalue = newnodatavalue;

    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: 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) {

        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
        rt_raster_destroy(newrast);

        elog(ERROR, "RASTER_mapAlgebraFct: Invalid pixeltype");
        PG_RETURN_NULL();
    }    
    
    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Pixeltype set to %s",
        rt_pixtype_name(newpixeltype));

    /* Get the name of the callback user function for raster values */
    if (PG_ARGISNULL(3)) {

        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
        rt_raster_destroy(newrast);

        elog(ERROR, "RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
        PG_RETURN_NULL();
    }

    oid = PG_GETARG_OID(3);
    if (oid == InvalidOid) {

        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
        rt_raster_destroy(newrast);

        elog(ERROR, "RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
        PG_RETURN_NULL();
    }

    fmgr_info(oid, &cbinfo);

    /* function cannot return set */
    if (cbinfo.fn_retset) {

        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
        rt_raster_destroy(newrast);

        elog(ERROR, "RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
        PG_RETURN_NULL();
    }
    /* function should have correct # of args */
    else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {

        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
        rt_raster_destroy(newrast);

        elog(ERROR, "RASTER_mapAlgebraFct: Function does not have two or three input parameters");
        PG_RETURN_NULL();
    }

    if (cbinfo.fn_nargs == 2)
        k = 1;
    else 
        k = 2;

    if (func_volatile(oid) == 'v') {
        elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
    }

    /* prep function call data */
#if POSTGIS_PGSQL_VERSION <= 90
    InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL);
#else
    InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
#endif
    memset(cbdata.argnull, FALSE, sizeof(bool) * cbinfo.fn_nargs);
    
    /* check that the function isn't strict if the args are null. */
    if (PG_ARGISNULL(4)) {
        if (cbinfo.fn_strict) {

            rt_raster_destroy(raster);
            PG_FREE_IF_COPY(pgraster, 0);
            rt_raster_destroy(newrast);

            elog(ERROR, "RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
            PG_RETURN_NULL();
        }

        cbdata.arg[k] = (Datum)NULL;
        cbdata.argnull[k] = TRUE;
    }
    else {
        cbdata.arg[k] = PG_GETARG_DATUM(4);
    }

    if (rt_band_get_isnodata_flag(band)) {

        POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Band is a nodata band, returning "
                "a raster filled with nodata");

        ret = rt_raster_generate_new_band(newrast, newpixeltype,
                newinitialvalue, TRUE, newnodatavalue, 0);

        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_mapAlgebraFct: Could not serialize raster");
            PG_RETURN_NULL();
        }

        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);               
    }


    ret = rt_raster_generate_new_band(newrast, newpixeltype,
            newinitialvalue, TRUE, newnodatavalue, 0);

    /* 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");

        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_mapAlgebraFct: Could not serialize raster");
            PG_RETURN_NULL();
        }

        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);      
    }
    
    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: Main computing loop (%d x %d)",
            width, height);

    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) {
                if (FLT_EQ(r, newnodatavalue)) {
                    if (cbinfo.fn_strict) {
                        POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
                        continue;
                    }
                    cbdata.argnull[0] = TRUE;
                    cbdata.arg[0] = (Datum)NULL;
                }
                else {
                    cbdata.argnull[0] = FALSE;
                    cbdata.arg[0] = Float8GetDatum(r);
                }

                /* Add pixel positions if callback has proper # of args */
                if (cbinfo.fn_nargs == 3) {
                    Datum d[2];
                    ArrayType *a;

                    d[0] = Int32GetDatum(x+1);
                    d[1] = Int32GetDatum(y+1);

                    a = construct_array(d, 2, INT4OID, sizeof(int32), true, 'i');

                    cbdata.argnull[1] = FALSE;
                    cbdata.arg[1] = PointerGetDatum(a);
                }

                POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: (%dx%d), r = %f",
                    x, y, r);
                   
                tmpnewval = FunctionCallInvoke(&cbdata);

                if (cbdata.isnull) {
                    newval = newnodatavalue;
                }
                else {
                    newval = DatumGetFloat8(tmpnewval);
                }

                POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: new value = %f", 
                    newval);
                
                rt_band_set_pixel(newband, x, y, newval, NULL);
            }

        }
    }
    
    /* The newrast band has been modified */

    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: 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();

    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFct: raster serialized");

    POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFct: returning raster");
    
    SET_VARSIZE(pgrtn, pgrtn->size);    
    PG_RETURN_POINTER(pgrtn);
}

Here is the call graph for this function: