RASTER_nMapAlgebraExpr()

Datum RASTER_nMapAlgebraExpr

(

PG_FUNCTION_ARGS

)

Definition at line 1335 of file rtpg_mapalgebra.c.

{
        MemoryContext mainMemCtx = CurrentMemoryContext;
        rtpg_nmapalgebraexpr_arg arg = NULL;
        rt_iterator itrset;
        uint16_t exprpos[3] = {1, 4, 5};
 
        int i = 0;
        int j = 0;
        int k = 0;
 
        int numraster = 0;
        int err = 0;
        int allnull = 0;
        int allempty = 0;
        int noband = 0;
        int len = 0;
 
        TupleDesc tupdesc;
        SPITupleTable *tuptable = NULL;
        HeapTuple tuple;
        Datum datum;
        bool isnull = FALSE;
 
        rt_raster raster = NULL;
        rt_band band = NULL;
        rt_pgraster *pgraster = NULL;
 
        const int argkwcount = 12;
        char *argkw[] = {
                "[rast.x]",
                "[rast.y]",
                "[rast.val]",
                "[rast]",
                "[rast1.x]",
                "[rast1.y]",
                "[rast1.val]",
                "[rast1]",
                "[rast2.x]",
                "[rast2.y]",
                "[rast2.val]",
                "[rast2]"
        };
 
        if (PG_ARGISNULL(0))
                PG_RETURN_NULL();
 
        /* init argument struct */
        arg = rtpg_nmapalgebraexpr_arg_init(argkwcount, argkw);
        if (arg == NULL) {
                elog(ERROR, "RASTER_nMapAlgebraExpr: Could not initialize argument structure");
                PG_RETURN_NULL();
        }
 
        /* let helper function process rastbandarg (0) */
        if (!rtpg_nmapalgebra_rastbandarg_process(arg->bandarg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
                rtpg_nmapalgebraexpr_arg_destroy(arg);
                elog(ERROR, "RASTER_nMapAlgebra: Could not process rastbandarg");
                PG_RETURN_NULL();
        }
 
        POSTGIS_RT_DEBUGF(4, "allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
 
        /* all rasters are NULL, return NULL */
        if (allnull == arg->bandarg->numraster) {
                elog(NOTICE, "All input rasters are NULL. Returning NULL");
                rtpg_nmapalgebraexpr_arg_destroy(arg);
                PG_RETURN_NULL();
        }
 
        /* only work on one or two rasters */
        if (arg->bandarg->numraster > 1)
                numraster = 2;
        else
                numraster = 1;
 
        /* pixel type (2) */
        if (!PG_ARGISNULL(2)) {
                char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
 
                /* Get the pixel type index */
                arg->bandarg->pixtype = rt_pixtype_index_from_name(pixtypename);
                if (arg->bandarg->pixtype == PT_END) {
                        rtpg_nmapalgebraexpr_arg_destroy(arg);
                        elog(ERROR, "RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
                        PG_RETURN_NULL();
                }
        }
        POSTGIS_RT_DEBUGF(4, "pixeltype: %d", arg->bandarg->pixtype);
 
        /* extent type (3) */
        if (!PG_ARGISNULL(3)) {
                char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(3))));
                arg->bandarg->extenttype = rt_util_extent_type(extenttypename);
        }
 
        if (arg->bandarg->extenttype == ET_CUSTOM) {
                if (numraster < 2) {
                        elog(NOTICE, "CUSTOM extent type not supported. Defaulting to FIRST");
                        arg->bandarg->extenttype = ET_FIRST;
                }
                else {
                        elog(NOTICE, "CUSTOM extent type not supported. Defaulting to INTERSECTION");
                        arg->bandarg->extenttype = ET_INTERSECTION;
                }
        }
        else if (numraster < 2)
                arg->bandarg->extenttype = ET_FIRST;
 
        POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->bandarg->extenttype);
 
        /* nodatanodataval (6) */
        if (!PG_ARGISNULL(6)) {
                arg->callback.nodatanodata.hasval = 1;
                arg->callback.nodatanodata.val = PG_GETARG_FLOAT8(6);
        }
 
        err = 0;
        /* all rasters are empty, return empty raster */
        if (allempty == arg->bandarg->numraster) {
                elog(NOTICE, "All input rasters are empty. Returning empty raster");
                err = 1;
        }
        /* all rasters don't have indicated band, return empty raster */
        else if (noband == arg->bandarg->numraster) {
                elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
                err = 1;
        }
        if (err) {
                rtpg_nmapalgebraexpr_arg_destroy(arg);
 
                raster = rt_raster_new(0, 0);
                if (raster == NULL) {
                        elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create empty raster");
                        PG_RETURN_NULL();
                }
 
                pgraster = rt_raster_serialize(raster);
                rt_raster_destroy(raster);
                if (!pgraster) PG_RETURN_NULL();
 
                SET_VARSIZE(pgraster, pgraster->size);
                PG_RETURN_POINTER(pgraster);
        }
 
        /* connect SPI */
        if (SPI_connect() != SPI_OK_CONNECT) {
                rtpg_nmapalgebraexpr_arg_destroy(arg);
                elog(ERROR, "RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
                PG_RETURN_NULL();
        }
 
        /*
                process expressions
 
                exprpos elements are:
                        1 - expression => spi_plan[0]
                        4 - nodata1expr => spi_plan[1]
                        5 - nodata2expr => spi_plan[2]
        */
        for (i = 0; i < arg->callback.exprcount; i++) {
                char *expr = NULL;
                char *tmp = NULL;
                char *sql = NULL;
                char place[12] = "$1";
 
                if (PG_ARGISNULL(exprpos[i]))
                        continue;
 
                expr = text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]));
                POSTGIS_RT_DEBUGF(3, "raw expr of argument #%d: %s", exprpos[i], expr);
 
                for (j = 0, k = 1; j < argkwcount; j++) {
                        /* attempt to replace keyword with placeholder */
                        len = 0;
                        tmp = rtpg_strreplace(expr, argkw[j], place, &len);
                        pfree(expr);
                        expr = tmp;
 
                        if (len) {
                                POSTGIS_RT_DEBUGF(4, "kw #%d (%s) at pos $%d", j, argkw[j], k);
                                arg->callback.expr[i].spi_argcount++;
                                arg->callback.expr[i].spi_argpos[j] = k++;
 
                                sprintf(place, "$%d", k);
                        }
                        else
                                arg->callback.expr[i].spi_argpos[j] = 0;
                }
 
                len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision");
                sql = (char *) palloc(len + 1);
                if (sql == NULL) {
                        rtpg_nmapalgebraexpr_arg_destroy(arg);
                        SPI_finish();
                        elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
                        PG_RETURN_NULL();
                }
 
                memcpy(sql, "SELECT (", strlen("SELECT ("));
                memcpy(sql + strlen("SELECT ("), expr, strlen(expr));
                memcpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision"));
                sql[len] = '\0';
 
                POSTGIS_RT_DEBUGF(3, "sql #%d: %s", exprpos[i], sql);
 
                /* prepared plan */
                if (arg->callback.expr[i].spi_argcount) {
                        Oid *argtype = (Oid *) palloc(arg->callback.expr[i].spi_argcount * sizeof(Oid));
                        POSTGIS_RT_DEBUGF(3, "expression parameter %d is a prepared plan", exprpos[i]);
                        if (argtype == NULL) {
                                pfree(sql);
                                rtpg_nmapalgebraexpr_arg_destroy(arg);
                                SPI_finish();
                                elog(ERROR, "RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
                                PG_RETURN_NULL();
                        }
 
                        /* specify datatypes of parameters */
                        for (j = 0, k = 0; j < argkwcount; j++) {
                                if (arg->callback.expr[i].spi_argpos[j] < 1) continue;
 
                                /* positions are INT4 */
                                if (
                                        (strstr(argkw[j], "[rast.x]") != NULL) ||
                                        (strstr(argkw[j], "[rast.y]") != NULL) ||
                                        (strstr(argkw[j], "[rast1.x]") != NULL) ||
                                        (strstr(argkw[j], "[rast1.y]") != NULL) ||
                                        (strstr(argkw[j], "[rast2.x]") != NULL) ||
                                        (strstr(argkw[j], "[rast2.y]") != NULL)
                                )
                                        argtype[k] = INT4OID;
                                /* everything else is FLOAT8 */
                                else
                                        argtype[k] = FLOAT8OID;
 
                                k++;
                        }
 
                        arg->callback.expr[i].spi_plan = SPI_prepare(sql, arg->callback.expr[i].spi_argcount, argtype);
                        pfree(argtype);
                        pfree(sql);
 
                        if (arg->callback.expr[i].spi_plan == NULL) {
                                rtpg_nmapalgebraexpr_arg_destroy(arg);
                                SPI_finish();
                                elog(ERROR, "RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
                                PG_RETURN_NULL();
                        }
                }
                /* no args, just execute query */
                else {
                        POSTGIS_RT_DEBUGF(3, "expression parameter %d has no args, simply executing", exprpos[i]);
                        err = SPI_execute(sql, TRUE, 0);
                        pfree(sql);
 
                        if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
                                rtpg_nmapalgebraexpr_arg_destroy(arg);
                                SPI_finish();
                                elog(ERROR, "RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
                                PG_RETURN_NULL();
                        }
 
                        /* get output of prepared plan */
                        tupdesc = SPI_tuptable->tupdesc;
                        tuptable = SPI_tuptable;
                        tuple = tuptable->vals[0];
 
                        datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
                        if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
                                if (SPI_tuptable) SPI_freetuptable(tuptable);
                                rtpg_nmapalgebraexpr_arg_destroy(arg);
                                SPI_finish();
                                elog(ERROR, "RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
                                PG_RETURN_NULL();
                        }
 
                        if (!isnull) {
                                arg->callback.expr[i].hasval = 1;
                                arg->callback.expr[i].val = DatumGetFloat8(datum);
                        }
 
                        if (SPI_tuptable) SPI_freetuptable(tuptable);
                }
        }
 
        /* determine nodataval and possibly pixtype */
        /* band to check */
        switch (arg->bandarg->extenttype) {
                case ET_LAST:
                case ET_SECOND:
                        if (numraster > 1)
                                i = 1;
                        else
                                i = 0;
                        break;
                default:
                        i = 0;
                        break;
        }
        /* find first viable band */
        if (!arg->bandarg->hasband[i]) {
                for (i = 0; i < numraster; i++) {
                        if (arg->bandarg->hasband[i])
                                break;
                }
                if (i >= numraster)
                        i = numraster - 1;
        }
        band = rt_raster_get_band(arg->bandarg->raster[i], arg->bandarg->nband[i]);
 
        /* set pixel type if PT_END */
        if (arg->bandarg->pixtype == PT_END)
                arg->bandarg->pixtype = rt_band_get_pixtype(band);
 
        /* set hasnodata and nodataval */
        arg->bandarg->hasnodata = 1;
        if (rt_band_get_hasnodata_flag(band))
                rt_band_get_nodata(band, &(arg->bandarg->nodataval));
        else
                arg->bandarg->nodataval = rt_band_get_min_value(band);
 
        POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->bandarg->pixtype), arg->bandarg->hasnodata, arg->bandarg->nodataval);
 
        /* init itrset */
        itrset = palloc(sizeof(struct rt_iterator_t) * numraster);
        if (itrset == NULL) {
                rtpg_nmapalgebraexpr_arg_destroy(arg);
                SPI_finish();
                elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
                PG_RETURN_NULL();
        }
 
        /* set itrset */
        for (i = 0; i < numraster; i++) {
                itrset[i].raster = arg->bandarg->raster[i];
                itrset[i].nband = arg->bandarg->nband[i];
                itrset[i].nbnodata = 1;
        }
 
        /* pass everything to iterator */
        err = rt_raster_iterator(
                itrset, numraster,
                arg->bandarg->extenttype, arg->bandarg->cextent,
                arg->bandarg->pixtype,
                arg->bandarg->hasnodata, arg->bandarg->nodataval,
                0, 0,
                NULL,
                &(arg->callback),
                rtpg_nmapalgebraexpr_callback,
                &raster
        );
 
        pfree(itrset);
        rtpg_nmapalgebraexpr_arg_destroy(arg);
 
        if (err != ES_NONE) {
                SPI_finish();
                elog(ERROR, "RASTER_nMapAlgebraExpr: Could not run raster iterator function");
                PG_RETURN_NULL();
        }
        else if (raster == NULL) {
                SPI_finish();
                PG_RETURN_NULL();
        }
 
        /* switch to prior memory context to ensure memory allocated in correct context */
        MemoryContextSwitchTo(mainMemCtx);
 
        pgraster = rt_raster_serialize(raster);
        rt_raster_destroy(raster);
 
        /* finish SPI */
        SPI_finish();
 
        if (!pgraster)
                PG_RETURN_NULL();
 
        SET_VARSIZE(pgraster, pgraster->size);
        PG_RETURN_POINTER(pgraster);
}

References ovdump::band, rtpg_nmapalgebraexpr_arg_t::bandarg, rtpg_nmapalgebraexpr_arg_t::callback, rtpg_nmapalgebra_arg_t::cextent, ES_NONE, ET_CUSTOM, ET_FIRST, ET_INTERSECTION, ET_LAST, ET_SECOND, rtpg_nmapalgebraexpr_callback_arg::expr, rtpg_nmapalgebraexpr_callback_arg::exprcount, rtpg_nmapalgebra_arg_t::extenttype, FALSE, rtpg_nmapalgebra_arg_t::hasband, rtpg_nmapalgebra_arg_t::hasnodata, rtpg_nmapalgebraexpr_callback_arg::hasval, rt_iterator_t::nband, rtpg_nmapalgebra_arg_t::nband, rt_iterator_t::nbnodata, rtpg_nmapalgebraexpr_callback_arg::nodatanodata, rtpg_nmapalgebra_arg_t::nodataval, rtpg_nmapalgebra_arg_t::numraster, rtpg_nmapalgebra_arg_t::pixtype, POSTGIS_RT_DEBUGF, PT_END, rt_iterator_t::raster, rtpg_nmapalgebra_arg_t::raster, rtrowdump::raster, rt_band_get_hasnodata_flag(), rt_band_get_min_value(), rt_band_get_nodata(), rt_band_get_pixtype(), rt_pixtype_index_from_name(), rt_pixtype_name(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_iterator(), rt_raster_new(), rt_raster_serialize(), rt_util_extent_type(), rtpg_nmapalgebra_rastbandarg_process(), rtpg_nmapalgebraexpr_arg_destroy(), rtpg_nmapalgebraexpr_arg_init(), rtpg_nmapalgebraexpr_callback(), rtpg_strreplace(), rtpg_strtoupper(), rtpg_trim(), rt_raster_serialized_t::size, rtpg_nmapalgebraexpr_callback_arg::spi_argcount, rtpg_nmapalgebraexpr_callback_arg::spi_argpos, rtpg_nmapalgebraexpr_callback_arg::spi_plan, rtgdalraster::sql, TRUE, and rtpg_nmapalgebraexpr_callback_arg::val.

Here is the call graph for this function: