RASTER_mapAlgebraFctNgb()

Datum RASTER_mapAlgebraFctNgb

(

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 neighborhood since the nodata value has already been set by the first optimization

Definition at line 5595 of file rtpg_mapalgebra.c.

References ovdump::band, ES_NONE, FALSE, FLT_NEQ, 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(), rtpg_strtoupper(), rt_raster_serialized_t::size, TRUE, rt_raster_serialized_t::width, pixval::x, and pixval::y.

Referenced by RASTER_mapAlgebraFct().

{
    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, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
    double r, rpix;
    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;
    ArrayType * neighborDatum;
    char * strFromText = NULL;
    text * txtNodataMode = NULL;
    text * txtCallbackParam = NULL;
    int intReplace = 0;
    float fltReplace = 0;
    bool valuereplace = false, pixelreplace, nNodataOnly = true, nNullSkip = false;
    Datum * neighborData = NULL;
    bool * neighborNulls = NULL;
    int neighborDims[2];
    int neighborLbs[2];
    int16 typlen;
    bool typbyval;
    char typalign;

    POSTGIS_RT_DEBUG(2, "RASTER_mapAlgebraFctNgb: 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_mapAlgebraFctNgb: Could not deserialize raster");
                        PG_RETURN_NULL();
    }

    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: 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_mapAlgebraFctNgb: 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_mapAlgebraFctNgb: 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_mapAlgebraFctNgb: Could not serialize raster");
            PG_RETURN_NULL();
        }

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

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

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

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

    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: Setting pixeltype...");

    if (PG_ARGISNULL(2)) {
        newpixeltype = rt_band_get_pixtype(band);
    }

    else {
        strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
        POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
        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_mapAlgebraFctNgb: Invalid pixeltype");
        PG_RETURN_NULL();
    }

    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Pixeltype set to %s (%d)",
        rt_pixtype_name(newpixeltype), newpixeltype);

    /* Get the name of the callback userfunction */
    if (PG_ARGISNULL(5)) {

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

        elog(ERROR, "RASTER_mapAlgebraFctNgb: Required function is missing");
        PG_RETURN_NULL();
    }

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

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

        elog(ERROR, "RASTER_mapAlgebraFctNgb: Got invalid function object id");
        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_mapAlgebraFctNgb: Function provided must return double precision not resultset");
        PG_RETURN_NULL();
    }
    /* function should have correct # of args */
    else if (cbinfo.fn_nargs != 3) {

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

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

    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, 3, InvalidOid, NULL, NULL);
    memset(cbdata.argnull, FALSE, sizeof(bool) * 3);
#else
    InitFunctionCallInfoData(*cbdata, &cbinfo, 3, 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(7)) {
        if (cbinfo.fn_strict) {

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

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

#if POSTGIS_PGSQL_VERSION < 120
        cbdata.arg[2] = (Datum)NULL;
        cbdata.argnull[2] = TRUE;
#else
        cbdata->args[2].value = (Datum)NULL;
        cbdata->args[2].isnull = TRUE;
#endif
    }
    else {
#if POSTGIS_PGSQL_VERSION < 120
            cbdata.arg[2] = PG_GETARG_DATUM(7);
#else
            cbdata->args[2].value = PG_GETARG_DATUM(7);
#endif
    }

    if (rt_band_get_isnodata_flag(band)) {

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

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

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


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

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

    /* Get the width of the neighborhood */
    if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
        elog(NOTICE, "Neighborhood width is NULL or <= 0. 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_mapAlgebraFctNgb: Could not serialize raster");
            PG_RETURN_NULL();
        }

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

    ngbwidth = PG_GETARG_INT32(3);
    winwidth = ngbwidth * 2 + 1;

    /* Get the height of the neighborhood */
    if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
        elog(NOTICE, "Neighborhood height is NULL or <= 0. 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_mapAlgebraFctNgb: Could not serialize raster");
            PG_RETURN_NULL();
        }

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

    ngbheight = PG_GETARG_INT32(4);
    winheight = ngbheight * 2 + 1;

    /* Get the type of NODATA behavior for the neighborhoods. */
    if (PG_ARGISNULL(6)) {
        elog(NOTICE, "Neighborhood NODATA behavior defaulting to 'ignore'");
        txtNodataMode = cstring_to_text("ignore");
    }
    else {
        txtNodataMode = PG_GETARG_TEXT_P(6);
    }

    txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
    SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
    memcpy((void *)VARDATA(txtCallbackParam), (void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);

    /* pass the nodata mode into the user function */
#if POSTGIS_PGSQL_VERSION < 120
    cbdata.arg[1] = CStringGetDatum(txtCallbackParam);
#else
    cbdata->args[1].value = CStringGetDatum(txtCallbackParam);
#endif

    strFromText = text_to_cstring(txtNodataMode);
    strFromText = rtpg_strtoupper(strFromText);

    if (strcmp(strFromText, "VALUE") == 0)
        valuereplace = true;
    else if (strcmp(strFromText, "IGNORE") != 0 && strcmp(strFromText, "NULL") != 0) {
        /* if the text is not "IGNORE" or "NULL", it may be a numerical value */
        if (sscanf(strFromText, "%d", &intReplace) <= 0 && sscanf(strFromText, "%f", &fltReplace) <= 0) {
            /* the value is NOT an integer NOR a floating point */
            elog(NOTICE, "Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', "
                "'NULL', or a numeric value. Returning new raster with the original band");

            /* clean up the nodatamode string */
            pfree(txtCallbackParam);
            pfree(strFromText);

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

            SET_VARSIZE(pgrtn, pgrtn->size);
            PG_RETURN_POINTER(pgrtn);
        }
    }
    else if (strcmp(strFromText, "NULL") == 0) {
        /* this setting means that the neighborhood should be skipped if any of the values are null */
        nNullSkip = true;
    }

    POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: Main computing loop (%d x %d)",
            width, height);

    /* Allocate room for the neighborhood. */
    neighborData = (Datum *)palloc(winwidth * winheight * sizeof(Datum));
    neighborNulls = (bool *)palloc(winwidth * winheight * sizeof(bool));

    /* The dimensions of the neighborhood array, for creating a multi-dimensional array. */
    neighborDims[0] = winwidth;
    neighborDims[1] = winheight;

    /* The lower bounds for the new multi-dimensional array. */
    neighborLbs[0] = 1;
    neighborLbs[1] = 1;

    /* Get information about the type of item in the multi-dimensional array (float8). */
    get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);

    for (x = 0 + ngbwidth; x < width - ngbwidth; x++) {
        for(y = 0 + ngbheight; y < height - ngbheight; y++) {
            /* populate an array with the pixel values in the neighborhood */
            nIndex = 0;
            nNullItems = 0;
            nNodataOnly = true;
            pixelreplace = false;
            if (valuereplace) {
                ret = rt_band_get_pixel(band, x, y, &rpix, NULL);
                if (ret == ES_NONE && FLT_NEQ(rpix, newnodatavalue)) {
                    pixelreplace = true;
                }
            }
            for (u = x - ngbwidth; u <= x + ngbwidth; u++) {
                for (v = y - ngbheight; v <= y + ngbheight; v++) {
                    ret = rt_band_get_pixel(band, u, v, &r, NULL);
                    if (ret == ES_NONE) {
                        if (FLT_NEQ(r, newnodatavalue)) {
                            /* If the pixel value for this neighbor cell is not NODATA */
                            neighborData[nIndex] = Float8GetDatum((double)r);
                            neighborNulls[nIndex] = false;
                            nNodataOnly = false;
                        }
                        else {
                            /* If the pixel value for this neighbor cell is NODATA */
                            if (valuereplace && pixelreplace) {
                                /* Replace the NODATA value with the currently processing pixel. */
                                neighborData[nIndex] = Float8GetDatum((double)rpix);
                                neighborNulls[nIndex] = false;
                                /* do not increment nNullItems, since the user requested that the  */
                                /* neighborhood replace NODATA values with the central pixel value */
                            }
                            else {
                                neighborData[nIndex] = PointerGetDatum(NULL);
                                neighborNulls[nIndex] = true;
                                nNullItems++;
                            }
                        }
                    }
                    else {
                        /* Fill this will NULL if we can't read the raster pixel. */
                        neighborData[nIndex] = PointerGetDatum(NULL);
                        neighborNulls[nIndex] = true;
                        nNullItems++;
                    }
                    /* Next neighbor position */
                    nIndex++;
                }
            }

            if (!(nNodataOnly ||                     /* neighborhood only contains NODATA -- OR -- */
                (nNullSkip && nNullItems > 0) ||     /* neighborhood should skip any NODATA cells, and a NODATA cell was detected -- OR -- */
                (valuereplace && nNullItems > 0))) { /* neighborhood should replace NODATA cells with the central pixel value, and a NODATA cell was detected */
                POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: (%dx%d), %dx%d neighborhood",
                    x, y, winwidth, winheight);

                neighborDatum = construct_md_array((void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
                    FLOAT8OID, typlen, typbyval, typalign);

#if POSTGIS_PGSQL_VERSION < 120
                /* Assign the neighbor matrix as the first argument to the user function */
                cbdata.arg[0] = PointerGetDatum(neighborDatum);

                /* Invoke the user function */
                tmpnewval = FunctionCallInvoke(&cbdata);

                /* Get the return value of the user function */
                if (cbdata.isnull) {
                    newval = newnodatavalue;
                }
#else
                /* Assign the neighbor matrix as the first argument to the user function */
                cbdata->args[0].value = PointerGetDatum(neighborDatum);

                /* Invoke the user function */
                tmpnewval = FunctionCallInvoke(cbdata);

                /* Get the return value of the user function */
                if (cbdata->isnull)
                {
                        newval = newnodatavalue;
                }
#endif
                else {
                    newval = DatumGetFloat8(tmpnewval);
                }

                POSTGIS_RT_DEBUGF(3, "RASTER_mapAlgebraFctNgb: new value = %f",
                    newval);

                rt_band_set_pixel(newband, x, y, newval, NULL);
            }

            /* reset the number of null items in the neighborhood */
            nNullItems = 0;
        }
    }


    /* clean up */
    pfree(neighborNulls);
    pfree(neighborData);
    pfree(strFromText);
    pfree(txtCallbackParam);

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

    /* The newrast band has been modified */

    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
    /* Serialize created raster */

    pgrtn = rt_raster_serialize(newrast);
    rt_raster_destroy(newrast);
    if (NULL == pgrtn)
        PG_RETURN_NULL();

    POSTGIS_RT_DEBUG(3, "RASTER_mapAlgebraFctNgb: raster serialized");
    POSTGIS_RT_DEBUG(4, "RASTER_mapAlgebraFctNgb: returning raster");

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

Here is the call graph for this function:

Here is the caller graph for this function: