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 5122 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;
    rt_pixtype newpixeltype;
    int ret = -1;
    Oid oid;
    FmgrInfo cbinfo;
#if POSTGIS_PGSQL_VERSION < 120
    FunctionCallInfoData cbdata;
#else
    LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS); /* Could be optimized */
#endif
    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 < 120
    InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
 
    cbdata.argnull[0] = FALSE;
    cbdata.argnull[1] = FALSE;
    cbdata.argnull[2] = FALSE;
#else
    InitFunctionCallInfoData(*cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
 
    cbdata->args[0].isnull = FALSE;
    cbdata->args[1].isnull = FALSE;
    cbdata->args[2].isnull = FALSE;
#endif
 
    /* 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();
        }
 
#if POSTGIS_PGSQL_VERSION < 120
        cbdata.arg[k] = (Datum)NULL;
        cbdata.argnull[k] = TRUE;
#else
        cbdata->args[k].value = (Datum)NULL;
        cbdata->args[k].isnull = TRUE;
#endif
    }
    else {
#if POSTGIS_PGSQL_VERSION < 120
            cbdata.arg[k] = PG_GETARG_DATUM(4);
#else
            cbdata->args[k].value = PG_GETARG_DATUM(4);
#endif
    }
 
    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;
                    }
#if POSTGIS_PGSQL_VERSION < 120
                    cbdata.argnull[0] = TRUE;
                    cbdata.arg[0] = (Datum)NULL;
#else
                    cbdata->args[0].isnull = TRUE;
                    cbdata->args[0].value = (Datum)NULL;
#endif
                }
                else {
#if POSTGIS_PGSQL_VERSION < 120
                        cbdata.argnull[0] = FALSE;
                        cbdata.arg[0] = Float8GetDatum(r);
#else
                        cbdata->args[0].isnull = FALSE;
                        cbdata->args[0].value = Float8GetDatum(r);
#endif
                }
 
                /* 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');
 
#if POSTGIS_PGSQL_VERSION < 120
                    cbdata.argnull[1] = FALSE;
                    cbdata.arg[1] = PointerGetDatum(a);
#else
                    cbdata->args[1].isnull = FALSE;
                    cbdata->args[1].value = PointerGetDatum(a);
#endif
                }
 
                POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFct: (%dx%d), r = %f",
                    x, y, r);
 
#if POSTGIS_PGSQL_VERSION < 120
                tmpnewval = FunctionCallInvoke(&cbdata);
 
                if (cbdata.isnull) {
                    newval = newnodatavalue;
                }
#else
                tmpnewval = FunctionCallInvoke(cbdata);
 
                if (cbdata->isnull)
                {
                        newval = newnodatavalue;
                }
#endif
                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);
}

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

Here is the call graph for this function: