RASTER_nMapAlgebra()

Datum RASTER_nMapAlgebra

(

PG_FUNCTION_ARGS

)

Definition at line 526 of file rtpg_mapalgebra.c.

{
        rtpg_nmapalgebra_arg arg = NULL;
        rt_iterator itrset;
        ArrayType *maskArray;
        Oid etype;
        Datum *maskElements;
        bool *maskNulls;
        int16 typlen;
        bool typbyval;
        char typalign;
        int ndims = 0;
        int num;
        int *maskDims;
        int x,y;
 
 
        int i = 0;
        int noerr = 0;
        int allnull = 0;
        int allempty = 0;
        int noband = 0;
 
        rt_raster raster = NULL;
        rt_band band = NULL;
        rt_pgraster *pgraster = NULL;
 
        POSTGIS_RT_DEBUG(3, "Starting...");
 
        if (PG_ARGISNULL(0))
                PG_RETURN_NULL();
 
        /* init argument struct */
        arg = rtpg_nmapalgebra_arg_init();
        if (arg == NULL) {
                elog(ERROR, "RASTER_nMapAlgebra: Could not initialize argument structure");
                PG_RETURN_NULL();
        }
 
        /* let helper function process rastbandarg (0) */
        if (!rtpg_nmapalgebra_rastbandarg_process(arg, PG_GETARG_ARRAYTYPE_P(0), &allnull, &allempty, &noband)) {
                rtpg_nmapalgebra_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->numraster) {
                elog(NOTICE, "All input rasters are NULL. Returning NULL");
                rtpg_nmapalgebra_arg_destroy(arg);
                PG_RETURN_NULL();
        }
 
        /* pixel type (2) */
        if (!PG_ARGISNULL(2)) {
                char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
 
                /* Get the pixel type index */
                arg->pixtype = rt_pixtype_index_from_name(pixtypename);
                if (arg->pixtype == PT_END) {
                        rtpg_nmapalgebra_arg_destroy(arg);
                        elog(ERROR, "RASTER_nMapAlgebra: Invalid pixel type: %s", pixtypename);
                        PG_RETURN_NULL();
                }
        }
 
        /* distancex (3) */
        if (!PG_ARGISNULL(3)){
                arg->distance[0] = PG_GETARG_INT32(3);
        }else{
                arg->distance[0] = 0;
        }
        /* distancey (4) */
        if (!PG_ARGISNULL(4)){
                arg->distance[1] = PG_GETARG_INT32(4);
        }else{
                arg->distance[1] = 0;
        }
        if (arg->distance[0] < 0 || arg->distance[1] < 0) {
                rtpg_nmapalgebra_arg_destroy(arg);
                elog(ERROR, "RASTER_nMapAlgebra: Distance for X and Y axis must be greater than or equal to zero");
                PG_RETURN_NULL();
        }
 
        /* extent type (5) */
        if (!PG_ARGISNULL(5)) {
                char *extenttypename = rtpg_strtoupper(rtpg_trim(text_to_cstring(PG_GETARG_TEXT_P(5))));
                arg->extenttype = rt_util_extent_type(extenttypename);
        }
        POSTGIS_RT_DEBUGF(4, "extenttype: %d", arg->extenttype);
 
        /* custom extent (6) */
        if (arg->extenttype == ET_CUSTOM) {
                if (PG_ARGISNULL(6)) {
                        elog(NOTICE, "Custom extent is NULL. Returning NULL");
                        rtpg_nmapalgebra_arg_destroy(arg);
                        PG_RETURN_NULL();
                }
 
                arg->pgcextent = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(6));
 
                /* only need the raster header */
                arg->cextent = rt_raster_deserialize(arg->pgcextent, TRUE);
                if (arg->cextent == NULL) {
                        rtpg_nmapalgebra_arg_destroy(arg);
                        elog(ERROR, "RASTER_nMapAlgebra: Could not deserialize custom extent");
                        PG_RETURN_NULL();
                }
                else if (rt_raster_is_empty(arg->cextent)) {
                        elog(NOTICE, "Custom extent is an empty raster. Returning empty raster");
                        rtpg_nmapalgebra_arg_destroy(arg);
 
                        raster = rt_raster_new(0, 0);
                        if (raster == NULL) {
                                elog(ERROR, "RASTER_nMapAlgebra: 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);
                }
        }
 
        /* mask (7) */
 
        if( PG_ARGISNULL(7) ){
                pfree(arg->mask);
                arg->mask = NULL;
        }
        else {
                maskArray = PG_GETARG_ARRAYTYPE_P(7);
                etype = ARR_ELEMTYPE(maskArray);
                get_typlenbyvalalign(etype,&typlen,&typbyval,&typalign);
 
                switch (etype) {
                        case FLOAT4OID:
                        case FLOAT8OID:
                                break;
                        default:
                                rtpg_nmapalgebra_arg_destroy(arg);
                                elog(ERROR,"RASTER_nMapAlgebra: Mask data type must be FLOAT8 or FLOAT4");
                                PG_RETURN_NULL();
                }
 
                ndims = ARR_NDIM(maskArray);
 
                if (ndims != 2) {
                        elog(ERROR, "RASTER_nMapAlgebra: Mask Must be a 2D array");
                        rtpg_nmapalgebra_arg_destroy(arg);
                        PG_RETURN_NULL();
                }
 
                maskDims = ARR_DIMS(maskArray);
 
                if (maskDims[0] % 2 == 0 || maskDims[1] % 2 == 0) {
                        elog(ERROR,"RASTER_nMapAlgebra: Mask dimensions must be odd");
                        rtpg_nmapalgebra_arg_destroy(arg);
                        PG_RETURN_NULL();
                }
 
                deconstruct_array(
                        maskArray,
                        etype,
                        typlen, typbyval,typalign,
                        &maskElements,&maskNulls,&num
                );
 
                if (num < 1 || num != (maskDims[0] * maskDims[1])) {
                        if (num) {
                                pfree(maskElements);
                                pfree(maskNulls);
                        }
                        elog(ERROR, "RASTER_nMapAlgebra: Could not deconstruct new values array");
                        rtpg_nmapalgebra_arg_destroy(arg);
                        PG_RETURN_NULL();
                }
 
                /* allocate mem for mask array */
                arg->mask->values = palloc(sizeof(double*)* maskDims[0]);
                arg->mask->nodata =  palloc(sizeof(int*)*maskDims[0]);
                for (i = 0; i < maskDims[0]; i++) {
                        arg->mask->values[i] = (double*) palloc(sizeof(double) * maskDims[1]);
                        arg->mask->nodata[i] = (int*) palloc(sizeof(int) * maskDims[1]);
                }
 
                /* place values in to mask */
                i = 0;
                for (y = 0; y < maskDims[0]; y++) {
                        for (x = 0; x < maskDims[1]; x++) {
                                if (maskNulls[i]) {
                                        arg->mask->values[y][x] = 0;
                                        arg->mask->nodata[y][x] = 1;
                                }
                                else {
                                        switch (etype) {
                                                case FLOAT4OID:
                                                        arg->mask->values[y][x] = (double) DatumGetFloat4(maskElements[i]);
                                                        arg->mask->nodata[y][x] = 0;
                                                        break;
                                                case FLOAT8OID:
                                                        arg->mask->values[y][x] = (double) DatumGetFloat8(maskElements[i]);
                                                        arg->mask->nodata[y][x] = 0;
                                        }
                                }
                                i++;
                        }
                }
 
                /*set mask dimensions*/
                arg->mask->dimx = maskDims[0];
                arg->mask->dimy = maskDims[1];
                if (maskDims[0] == 1 && maskDims[1] == 1) {
                        arg->distance[0] = 0;
                        arg->distance[1] = 0;
                }
                else {
                        arg->distance[0] = maskDims[0] % 2;
                        arg->distance[1] = maskDims[1] % 2;
                }
        }/*end if else argisnull*/
 
        /* (8) weighted boolean */
        if (PG_ARGISNULL(8) || !PG_GETARG_BOOL(8)) {
                if (arg->mask != NULL)
                        arg->mask->weighted = 0;
        }else{
                if(arg->mask !=NULL)
                        arg->mask->weighted = 1;
        }
 
        noerr = 1;
 
        /* all rasters are empty, return empty raster */
        if (allempty == arg->numraster) {
                elog(NOTICE, "All input rasters are empty. Returning empty raster");
                noerr = 0;
        }
        /* all rasters don't have indicated band, return empty raster */
        else if (noband == arg->numraster) {
                elog(NOTICE, "All input rasters do not have bands at indicated indexes. Returning empty raster");
                noerr = 0;
        }
        if (!noerr) {
                rtpg_nmapalgebra_arg_destroy(arg);
 
                raster = rt_raster_new(0, 0);
                if (raster == NULL) {
                        elog(ERROR, "RASTER_nMapAlgebra: 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);
        }
 
        /* do regprocedure last (1) */
        if (!PG_ARGISNULL(1) || get_fn_expr_argtype(fcinfo->flinfo, 1) == REGPROCEDUREOID) {
                POSTGIS_RT_DEBUG(4, "processing callbackfunc");
                arg->callback.ufc_noid = PG_GETARG_OID(1);
 
                /* get function info */
                fmgr_info(arg->callback.ufc_noid, &(arg->callback.ufl_info));
 
                /* function cannot return set */
                noerr = 0;
                if (arg->callback.ufl_info.fn_retset) {
                        noerr = 1;
                }
                /* function should have correct # of args */
                else if (arg->callback.ufl_info.fn_nargs != 3) {
                        noerr = 2;
                }
 
                /* check that callback function return type is supported */
                if (
                        get_func_result_type(
                                arg->callback.ufc_noid,
                                &(arg->callback.ufc_rettype),
                                NULL
                        ) != TYPEFUNC_SCALAR
                ) {
                        noerr = 3;
                }
 
                if (!(
                        arg->callback.ufc_rettype == FLOAT8OID ||
                        arg->callback.ufc_rettype == FLOAT4OID ||
                        arg->callback.ufc_rettype == INT4OID ||
                        arg->callback.ufc_rettype == INT2OID
                )) {
                        noerr = 4;
                }
 
                /*
                        TODO: consider adding checks of the userfunction parameters
                                should be able to use get_fn_expr_argtype() of fmgr.c
                */
 
                if (noerr != 0) {
                        rtpg_nmapalgebra_arg_destroy(arg);
                        switch (noerr) {
                                case 4:
                                        elog(ERROR, "RASTER_nMapAlgebra: Function provided must return a double precision, float, int or smallint");
                                        break;
                                case 3:
                                        elog(ERROR, "RASTER_nMapAlgebra: Function provided must return scalar (double precision, float, int, smallint)");
                                        break;
                                case 2:
                                        elog(ERROR, "RASTER_nMapAlgebra: Function provided must have three input parameters");
                                        break;
                                case 1:
                                        elog(ERROR, "RASTER_nMapAlgebra: Function provided must return double precision, not resultset");
                                        break;
                        }
                        PG_RETURN_NULL();
                }
 
                if (func_volatile(arg->callback.ufc_noid) == 'v')
                        elog(NOTICE, "Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
 
                /* prep function call data */
                InitFunctionCallInfoData(*(arg->callback.ufc_info),
                                         &(arg->callback.ufl_info),
                                         arg->callback.ufl_info.fn_nargs,
                                         InvalidOid,
                                         NULL,
                                         NULL);
 
                arg->callback.ufc_info->args[0].isnull = FALSE;
                arg->callback.ufc_info->args[1].isnull = FALSE;
                arg->callback.ufc_info->args[2].isnull = FALSE;
                /* userargs (7) */
                if (!PG_ARGISNULL(9))
                        arg->callback.ufc_info->args[2].value = PG_GETARG_DATUM(9);
                else {
      if (arg->callback.ufl_info.fn_strict) {
                                /* build and assign an empty TEXT array */
                                /* TODO: manually free the empty array? */
                                arg->callback.ufc_info->args[2].value = PointerGetDatum(construct_empty_array(TEXTOID));
                                arg->callback.ufc_info->args[2].isnull = FALSE;
      }
                        else {
                                arg->callback.ufc_info->args[2].value = (Datum)NULL;
                                arg->callback.ufc_info->args[2].isnull = TRUE;
                        }
                }
        }
        else {
                rtpg_nmapalgebra_arg_destroy(arg);
                elog(ERROR, "RASTER_nMapAlgebra: callbackfunc must be provided");
                PG_RETURN_NULL();
        }
 
        /* determine nodataval and possibly pixtype */
        /* band to check */
        switch (arg->extenttype) {
                case ET_LAST:
                        i = arg->numraster - 1;
                        break;
                case ET_SECOND:
                        i = (arg->numraster > 1) ? 1 : 0;
                        break;
                default:
                        i = 0;
                        break;
        }
        /* find first viable band */
        if (!arg->hasband[i]) {
                for (i = 0; i < arg->numraster; i++) {
                        if (arg->hasband[i])
                                break;
                }
                if (i >= arg->numraster)
                        i = arg->numraster - 1;
        }
        band = rt_raster_get_band(arg->raster[i], arg->nband[i]);
 
        /* set pixel type if PT_END */
        if (arg->pixtype == PT_END)
                arg->pixtype = rt_band_get_pixtype(band);
 
        /* set hasnodata and nodataval */
        arg->hasnodata = 1;
        if (rt_band_get_hasnodata_flag(band))
                rt_band_get_nodata(band, &(arg->nodataval));
        else
                arg->nodataval = rt_band_get_min_value(band);
 
        POSTGIS_RT_DEBUGF(4, "pixtype, hasnodata, nodataval: %s, %d, %f", rt_pixtype_name(arg->pixtype), arg->hasnodata, arg->nodataval);
 
        /* init itrset */
        itrset = palloc(sizeof(struct rt_iterator_t) * arg->numraster);
        if (itrset == NULL) {
                rtpg_nmapalgebra_arg_destroy(arg);
                elog(ERROR, "RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
                PG_RETURN_NULL();
        }
 
        /* set itrset */
        for (i = 0; i < arg->numraster; i++) {
                itrset[i].raster = arg->raster[i];
                itrset[i].nband = arg->nband[i];
                itrset[i].nbnodata = 1;
        }
 
        /* pass everything to iterator */
        noerr = rt_raster_iterator(
                itrset, arg->numraster,
                arg->extenttype, arg->cextent,
                arg->pixtype,
                arg->hasnodata, arg->nodataval,
                arg->distance[0], arg->distance[1],
                arg->mask,
                &(arg->callback),
                rtpg_nmapalgebra_callback,
                &raster
        );
 
        /* cleanup */
        pfree(itrset);
        rtpg_nmapalgebra_arg_destroy(arg);
 
        if (noerr != ES_NONE) {
                elog(ERROR, "RASTER_nMapAlgebra: Could not run raster iterator function");
                PG_RETURN_NULL();
        }
        else if (raster == NULL)
                PG_RETURN_NULL();
 
        pgraster = rt_raster_serialize(raster);
        rt_raster_destroy(raster);
 
        POSTGIS_RT_DEBUG(3, "Finished");
 
        if (!pgraster)
                PG_RETURN_NULL();
 
        SET_VARSIZE(pgraster, pgraster->size);
        PG_RETURN_POINTER(pgraster);
}

References ovdump::band, rtpg_nmapalgebra_arg_t::callback, rtpg_nmapalgebra_arg_t::cextent, rt_mask_t::dimx, rt_mask_t::dimy, rtpg_nmapalgebra_arg_t::distance, ES_NONE, ET_CUSTOM, ET_LAST, ET_SECOND, rtpg_nmapalgebra_arg_t::extenttype, FALSE, rtpg_nmapalgebra_arg_t::hasband, rtpg_nmapalgebra_arg_t::hasnodata, rtpg_nmapalgebra_arg_t::mask, rt_iterator_t::nband, rtpg_nmapalgebra_arg_t::nband, rt_iterator_t::nbnodata, rt_mask_t::nodata, rtpg_nmapalgebra_arg_t::nodataval, rtpg_nmapalgebra_arg_t::numraster, rtpg_nmapalgebra_arg_t::pgcextent, rtpg_nmapalgebra_arg_t::pixtype, POSTGIS_RT_DEBUG, 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_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_is_empty(), rt_raster_iterator(), rt_raster_new(), rt_raster_serialize(), rt_util_extent_type(), rtpg_nmapalgebra_arg_destroy(), rtpg_nmapalgebra_arg_init(), rtpg_nmapalgebra_callback(), rtpg_nmapalgebra_rastbandarg_process(), rtpg_strtoupper(), rtpg_trim(), rt_raster_serialized_t::size, TRUE, rtpg_nmapalgebra_callback_arg::ufc_info, rtpg_nmapalgebra_callback_arg::ufc_noid, rtpg_nmapalgebra_callback_arg::ufc_rettype, rtpg_nmapalgebra_callback_arg::ufl_info, rt_mask_t::values, rt_mask_t::weighted, pixval::x, and pixval::y.

Here is the call graph for this function: