rt_raster_iterator()

rt_errorstate rt_raster_iterator	(	rt_iterator	itrset,
		uint16_t	itrcount,
		rt_extenttype	extenttype,
		rt_raster	customextent,
		rt_pixtype	pixtype,
		uint8_t	hasnodata,
		double	nodataval,
		uint16_t	distancex,
		uint16_t	distancey,
		rt_mask	mask,
		void *	userarg,
		int(*)(rt_iterator_arg arg, void *userarg, double *value, int *nodata)	callback,
		rt_raster *	rtnraster
	)

n-raster iterator.

The raster returned should be freed by the caller

Parameters

itrset	: set of rt_iterator objects.
itrcount	: number of objects in itrset.
extenttype	: type of extent for the output raster.
customextent	: raster specifying custom extent. is only used if extenttype is ET_CUSTOM.
pixtype	: the desired pixel type of the output raster's band.
hasnodata	: indicates if the band has nodata value
nodataval	: the nodata value, will be appropriately truncated to fit the pixtype size.
distancex	: the number of pixels around the specified pixel along the X axis
distancey	: the number of pixels around the specified pixel along the Y axis
mask	: the object of mask
userarg	: pointer to any argument that is passed as-is to callback.
callback	: callback function for actual processing of pixel values.
*rtnraster	: return one band raster from iterator process

The callback function must have the following signature.

int FNAME(rt_iterator_arg arg, void *userarg, double *value, int *nodata)

The callback function must return zero (error) or non-zero (success) indicating whether the function ran successfully. The parameters passed to the callback function are as follows.

• rt_iterator_arg arg: struct containing pixel values, NODATA flags and metadata
• void *userarg: NULL or calling function provides to rt_raster_iterator() for use by callback function
• double *value: value of pixel to be burned by rt_raster_iterator()
• int *nodata: flag (0 or 1) indicating that pixel to be burned is NODATA

Returns

ES_NONE on success, ES_ERROR on error

Definition at line 812 of file rt_mapalgebra.c.

  {
        /* output raster */
        rt_raster rtnrast = NULL;
        /* output raster's band */
        rt_band rtnband = NULL;
 
        /* working raster */
        rt_raster rast = NULL;
 
        _rti_iterator_arg _param = NULL;
        int allnull = 0;
        int allempty = 0;
        int aligned = 0;
        double offset[4] = {0.};
        rt_pixel npixels;
 
        int i = 0;
        int status = 0;
        int inextent = 0;
        int x = 0;
        int y = 0;
        int _x = 0;
        int _y = 0;
 
        int _width = 0;
        int _height = 0;
 
        double minval;
        double value;
        int isnodata;
        int nodata;
 
        RASTER_DEBUG(3, "Starting...");
 
        assert(itrset != NULL && itrcount > 0);
        assert(rtnraster != NULL);
 
        /* init rtnraster to NULL */
        *rtnraster = NULL;
 
        /* check that callback function is not NULL */
        if (callback == NULL) {
                rterror("rt_raster_iterator: Callback function not provided");
                return ES_ERROR;
        }
 
        /* check that custom extent is provided if extenttype = ET_CUSTOM */
        if (extenttype == ET_CUSTOM && rt_raster_is_empty(customextent)) {
                rterror("rt_raster_iterator: Custom extent cannot be empty if extent type is ET_CUSTOM");
                return ES_ERROR;
        }
 
        /* check that pixtype != PT_END */
        if (pixtype == PT_END) {
                rterror("rt_raster_iterator: Pixel type cannot be PT_END");
                return ES_ERROR;
        }
 
        /* initialize _param */
        if ((_param = _rti_iterator_arg_init()) == NULL) {
                rterror("rt_raster_iterator: Could not initialize internal variables");
                return ES_ERROR;
        }
 
        /* fill _param */
        if (!_rti_iterator_arg_populate(_param, itrset, itrcount, distancex, distancey, &allnull, &allempty)) {
                rterror("rt_raster_iterator: Could not populate for internal variables");
                _rti_iterator_arg_destroy(_param);
                return ES_ERROR;
        }
 
        /* shortcut if all null, return NULL */
        if (allnull == itrcount) {
                RASTER_DEBUG(3, "all rasters are NULL, returning NULL");
 
                _rti_iterator_arg_destroy(_param);
 
                return ES_NONE;
        }
        /* shortcut if all empty, return empty raster */
        else if (allempty == itrcount) {
                RASTER_DEBUG(3, "all rasters are empty, returning empty raster");
 
                _rti_iterator_arg_destroy(_param);
 
                rtnrast = rt_raster_new(0, 0);
                if (rtnrast == NULL) {
                        rterror("rt_raster_iterator: Could not create empty raster");
                        return ES_ERROR;
                }
                rt_raster_set_scale(rtnrast, 0, 0);
 
                *rtnraster = rtnrast;
                return ES_NONE;
        }
 
        /* check that all rasters are aligned */
        RASTER_DEBUG(3, "checking alignment of all rasters");
        rast = NULL;
 
        /* find raster to use as reference */
        /* use custom if provided */
        if (extenttype == ET_CUSTOM) {
                RASTER_DEBUG(4, "using custom extent as reference raster");
                rast = customextent;
        }
        /* use first valid one in _param->raster */
        else {
                for (i = 0; i < itrcount; i++) {
                        if (!_param->isempty[i]) {
                                RASTER_DEBUGF(4, "using raster at index %d as reference raster", i);
                                rast = _param->raster[i];
                                break;
                        }
                }
        }
 
        /* no rasters found, SHOULD NEVER BE HERE! */
        if (rast == NULL) {
                rterror("rt_raster_iterator: Could not find reference raster to use for alignment tests");
 
                _rti_iterator_arg_destroy(_param);
 
                return ES_ERROR;
        }
 
        do {
                aligned = 1;
 
                /* check custom first if set. also skip if rasters are the same */
                if (extenttype == ET_CUSTOM && rast != customextent) {
                        if (rt_raster_same_alignment(rast, customextent, &aligned, NULL) != ES_NONE) {
                                rterror("rt_raster_iterator: Could not test for alignment between reference raster and custom extent");
 
                                _rti_iterator_arg_destroy(_param);
 
                                return ES_ERROR;
                        }
 
                        RASTER_DEBUGF(5, "custom extent alignment: %d", aligned);
                        if (!aligned)
                                break;
                }
 
                for (i = 0; i < itrcount; i++) {
                        /* skip NULL rasters and if rasters are the same */
                        if (_param->isempty[i] || rast == _param->raster[i])
                                continue;
 
                        if (rt_raster_same_alignment(rast, _param->raster[i], &aligned, NULL) != ES_NONE) {
                                rterror("rt_raster_iterator: Could not test for alignment between reference raster and raster %d", i);
 
                                _rti_iterator_arg_destroy(_param);
 
                                return ES_ERROR;
                        }
                        RASTER_DEBUGF(5, "raster at index %d alignment: %d", i, aligned);
 
                        /* abort checking since a raster isn't aligned */
                        if (!aligned)
                                break;
                }
        }
        while (0);
 
        /* not aligned, error */
        if (!aligned) {
                rterror("rt_raster_iterator: The set of rasters provided (custom extent included, if appropriate) do not have the same alignment");
 
                _rti_iterator_arg_destroy(_param);
 
                return ES_ERROR;
        }
 
        /* use extenttype to build output raster (no bands though) */
        i = -1;
        switch (extenttype) {
                case ET_INTERSECTION:
                case ET_UNION:
                        /* make copy of first "real" raster */
                        rtnrast = rtalloc(sizeof(struct rt_raster_t));
                        if (rtnrast == NULL) {
                                rterror("rt_raster_iterator: Could not allocate memory for output raster");
 
                                _rti_iterator_arg_destroy(_param);
 
                                return ES_ERROR;
                        }
 
                        for (i = 0; i < itrcount; i++) {
                                if (!_param->isempty[i]) {
                                        memcpy(rtnrast, _param->raster[i], sizeof(struct rt_raster_t));
                                        break;
                                }
                        }
                        rtnrast->numBands = 0;
                        rtnrast->bands = NULL;
 
                        /* get extent of output raster */
                        rast = NULL;
                        for (i = i + 1; i < itrcount; i++) {
                                if (_param->isempty[i])
                                        continue;
 
                                status = rt_raster_from_two_rasters(rtnrast, _param->raster[i], extenttype, &rast, NULL);
                                rtdealloc(rtnrast);
 
                                if (rast == NULL || status != ES_NONE) {
                                        rterror("rt_raster_iterator: Could not compute %s extent of rasters",
                                                extenttype == ET_UNION ? "union" : "intersection"
                                        );
 
                                        _rti_iterator_arg_destroy(_param);
 
                                        return ES_ERROR;
                                }
                                else if (rt_raster_is_empty(rast)) {
                                        rtinfo("rt_raster_iterator: Computed raster for %s extent is empty",
                                                extenttype == ET_UNION ? "union" : "intersection"
                                        );
 
                                        _rti_iterator_arg_destroy(_param);
 
                                        *rtnraster = rast;
                                        return ES_NONE;
                                }
 
                                rtnrast = rast;
                                rast = NULL;
                        }
 
                        break;
                /*
                        first, second and last have similar checks
                        and continue into custom
                */
                case ET_FIRST:
                        i = 0;
                        /* FALLTHROUGH */
                case ET_SECOND:
                        if (i < 0) {
                                if (itrcount < 2)
                                        i = 0;
                                else
                                        i = 1;
                        }
                        /* FALLTHROUGH */
                case ET_LAST:
                        if (i < 0) i = itrcount - 1;
 
                        /* input raster is null, return NULL */
                        if (_param->raster[i] == NULL) {
                                RASTER_DEBUGF(3, "returning NULL as %s raster is NULL and extent type is ET_%s",
                                        (i == 0 ? "first" : (i == 1 ? "second" : "last")),
                                        (i == 0 ? "FIRST" : (i == 1 ? "SECOND" : "LAST"))
                                );
 
                                _rti_iterator_arg_destroy(_param);
 
                                return ES_NONE;
                        }
                        /* input raster is empty, return empty raster */
                        else if (_param->isempty[i]) {
                                RASTER_DEBUGF(3, "returning empty raster as %s raster is empty and extent type is ET_%s",
                                        (i == 0 ? "first" : (i == 1 ? "second" : "last")),
                                        (i == 0 ? "FIRST" : (i == 1 ? "SECOND" : "LAST"))
                                );
 
                                _rti_iterator_arg_destroy(_param);
 
                                rtnrast = rt_raster_new(0, 0);
                                if (rtnrast == NULL) {
                                        rterror("rt_raster_iterator: Could not create empty raster");
                                        return ES_ERROR;
                                }
                                rt_raster_set_scale(rtnrast, 0, 0);
 
                                *rtnraster = rtnrast;
                                return ES_NONE;
                        }
                        /* FALLTHROUGH */
                /* copy the custom extent raster */
                case ET_CUSTOM:
                        rtnrast = rtalloc(sizeof(struct rt_raster_t));
                        if (rtnrast == NULL) {
                                rterror("rt_raster_iterator: Could not allocate memory for output raster");
 
                                _rti_iterator_arg_destroy(_param);
 
                                return ES_ERROR;
                        }
 
                        switch (extenttype) {
                                case ET_CUSTOM:
                                        memcpy(rtnrast, customextent, sizeof(struct rt_raster_t));
                                        break;
                                /* first, second, last */
                                default:
                                        memcpy(rtnrast, _param->raster[i], sizeof(struct rt_raster_t));
                                        break;
                        }
                        rtnrast->numBands = 0;
                        rtnrast->bands = NULL;
                        break;
        }
 
        _width = rt_raster_get_width(rtnrast);
        _height = rt_raster_get_height(rtnrast);
 
        RASTER_DEBUGF(4, "rtnrast (width, height, ulx, uly, scalex, scaley, skewx, skewy, srid) = (%d, %d, %f, %f, %f, %f, %f, %f, %d)",
                _width,
                _height,
                rt_raster_get_x_offset(rtnrast),
                rt_raster_get_y_offset(rtnrast),
                rt_raster_get_x_scale(rtnrast),
                rt_raster_get_y_scale(rtnrast),
                rt_raster_get_x_skew(rtnrast),
                rt_raster_get_y_skew(rtnrast),
                rt_raster_get_srid(rtnrast)
        );
 
        /* init values and NODATA for use with empty rasters */
        if (!_rti_iterator_arg_empty_init(_param)) {
                rterror("rt_raster_iterator: Could not initialize empty values and NODATA");
 
                _rti_iterator_arg_destroy(_param);
                rt_raster_destroy(rtnrast);
 
                return ES_ERROR;
        }
 
        /* create output band */
        if (rt_raster_generate_new_band(
                rtnrast,
                pixtype,
                nodataval,
                hasnodata, nodataval,
                0
        ) < 0) {
                rterror("rt_raster_iterator: Could not add new band to output raster");
 
                _rti_iterator_arg_destroy(_param);
                rt_raster_destroy(rtnrast);
 
                return ES_ERROR;
        }
 
        /* get output band */
        rtnband = rt_raster_get_band(rtnrast, 0);
        if (rtnband == NULL) {
                rterror("rt_raster_iterator: Could not get new band from output raster");
 
                _rti_iterator_arg_destroy(_param);
                rt_raster_destroy(rtnrast);
 
                return ES_ERROR;
        }
 
        /* output band's minimum value */
        minval = rt_band_get_min_value(rtnband);
 
        /* initialize argument for callback function */
        if (!_rti_iterator_arg_callback_init(_param)) {
                rterror("rt_raster_iterator: Could not initialize callback function argument");
 
                _rti_iterator_arg_destroy(_param);
                rt_band_destroy(rtnband);
                rt_raster_destroy(rtnrast);
 
                return ES_ERROR;
        }
 
        /* fill _param->offset */
        for (i = 0; i < itrcount; i++) {
                if (_param->isempty[i])
                        continue;
 
                status = rt_raster_from_two_rasters(rtnrast, _param->raster[i], ET_FIRST, &rast, offset);
                rtdealloc(rast);
                if (status != ES_NONE) {
                        rterror("rt_raster_iterator: Could not compute raster offsets");
 
                        _rti_iterator_arg_destroy(_param);
                        rt_band_destroy(rtnband);
                        rt_raster_destroy(rtnrast);
 
                        return ES_ERROR;
                }
 
                _param->offset[i][0] = offset[2];
                _param->offset[i][1] = offset[3];
                RASTER_DEBUGF(4, "rast %d offset: %f %f", i, offset[2], offset[3]);
        }
 
        /* loop over each pixel (POI) of output raster */
        /* _x,_y are for output raster */
        /* x,y are for input raster */
        for (_y = 0; _y < _height; _y++) {
                for (_x = 0; _x < _width; _x++) {
                        RASTER_DEBUGF(4, "iterating output pixel (x, y) = (%d, %d)", _x, _y);
                        _param->arg->dst_pixel[0] = _x;
                        _param->arg->dst_pixel[1] = _y;
 
                        /* loop through each input raster */
                        for (i = 0; i < itrcount; i++) {
                                RASTER_DEBUGF(4, "raster %d", i);
 
                                /*
                                        empty raster
                                        OR band does not exist and flag set to use NODATA
                                        OR band is NODATA
                                */
                                if (
                                        _param->isempty[i] ||
                                        (_param->band.rtband[i] == NULL && itrset[i].nbnodata) ||
                                        _param->band.isnodata[i]
                                ) {
                                        RASTER_DEBUG(4, "empty raster, band does not exist or band is NODATA. using empty values and NODATA");
 
                                        x = _x;
                                        y = _y;
 
                                        _param->arg->values[i] = _param->empty.values;
                                        _param->arg->nodata[i] = _param->empty.nodata;
 
                                        continue;
                                }
 
                                /* input raster's X,Y */
                                x = _x - (int) _param->offset[i][0];
                                y = _y - (int) _param->offset[i][1];
                                RASTER_DEBUGF(4, "source pixel (x, y) = (%d, %d)", x, y);
 
                                _param->arg->src_pixel[i][0] = x;
                                _param->arg->src_pixel[i][1] = y;
 
                                /* neighborhood */
                                npixels = NULL;
                                status = 0;
                                if (distancex > 0 && distancey > 0) {
                                        RASTER_DEBUG(4, "getting neighborhood");
 
                                        status = rt_band_get_nearest_pixel(
                                                _param->band.rtband[i],
                                                x, y,
                                                distancex, distancey,
                                                1,
                                                &npixels
                                        );
                                        if (status < 0) {
                                                rterror("rt_raster_iterator: Could not get pixel neighborhood");
 
                                                _rti_iterator_arg_destroy(_param);
                                                rt_band_destroy(rtnband);
                                                rt_raster_destroy(rtnrast);
 
                                                return ES_ERROR;
                                        }
                                }
 
                                /* get value of POI */
                                /* get pixel's value */
                                if (
                                        (x >= 0 && x < _param->width[i]) &&
                                        (y >= 0 && y < _param->height[i])
                                ) {
                                        RASTER_DEBUG(4, "getting value of POI");
                                        if (rt_band_get_pixel(
                                                _param->band.rtband[i],
                                                x, y,
                                                &value,
                                                &isnodata
                                        ) != ES_NONE) {
                                                rterror("rt_raster_iterator: Could not get the pixel value of band");
 
                                                _rti_iterator_arg_destroy(_param);
                                                rt_band_destroy(rtnband);
                                                rt_raster_destroy(rtnrast);
 
                                                return ES_ERROR;
                                        }
                                        inextent = 1;
                                }
                                /* outside band extent, set to NODATA */
                                else {
                                        RASTER_DEBUG(4, "Outside band extent, setting value to NODATA");
                                        /* has NODATA, use NODATA */
                                        if (_param->band.hasnodata[i])
                                                value = _param->band.nodataval[i];
                                        /* no NODATA, use min possible value */
                                        else
                                                value = _param->band.minval[i];
 
                                        inextent = 0;
                                        isnodata = 1;
                                }
 
                                /* add pixel to neighborhood */
                                status++;
                                if (status > 1)
                                        npixels = (rt_pixel) rtrealloc(npixels, sizeof(struct rt_pixel_t) * status);
                                else
                                        npixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t));
 
                                if (npixels == NULL) {
                                        rterror("rt_raster_iterator: Could not reallocate memory for neighborhood");
 
                                        _rti_iterator_arg_destroy(_param);
                                        rt_band_destroy(rtnband);
                                        rt_raster_destroy(rtnrast);
 
                                        return ES_ERROR;
                                }
 
                                npixels[status - 1].x = x;
                                npixels[status - 1].y = y;
                                npixels[status - 1].nodata = 1;
                                npixels[status - 1].value = value;
 
                                /* set nodata flag */
                                if ((!_param->band.hasnodata[i] && inextent) || !isnodata) {
                                        npixels[status - 1].nodata = 0;
                                }
                                RASTER_DEBUGF(4, "value, nodata: %f, %d", value, npixels[status - 1].nodata);
 
                                /* convert set of rt_pixel to 2D array */
                                status = rt_pixel_set_to_array(
                                        npixels, status, mask,
                                        x, y,
                                        distancex, distancey,
                                        &(_param->arg->values[i]),
                                        &(_param->arg->nodata[i]),
                                        NULL, NULL
                                );
                                rtdealloc(npixels);
                                if (status != ES_NONE) {
                                        rterror("rt_raster_iterator: Could not create 2D array of neighborhood");
 
                                        _rti_iterator_arg_destroy(_param);
                                        rt_band_destroy(rtnband);
                                        rt_raster_destroy(rtnrast);
 
                                        return ES_ERROR;
                                }
                        }
 
                        /* callback */
                        RASTER_DEBUG(4, "calling callback function");
                        value = 0;
                        nodata = 0;
                        status = callback(_param->arg, userarg, &value, &nodata);
 
                        /* free memory from callback */
                        _rti_iterator_arg_callback_clean(_param);
 
                        /* handle callback status */
                        if (status == 0) {
                                rterror("rt_raster_iterator: Callback function returned an error");
 
                                _rti_iterator_arg_destroy(_param);
                                rt_band_destroy(rtnband);
                                rt_raster_destroy(rtnrast);
 
                                return ES_ERROR;
                        }
 
                        /* burn value to pixel */
                        status = 0;
                        if (!nodata) {
                                status = rt_band_set_pixel(rtnband, _x, _y, value, NULL);
                                RASTER_DEBUGF(4, "burning pixel (%d, %d) with value: %f", _x, _y, value);
                        }
                        else if (!hasnodata) {
                                status = rt_band_set_pixel(rtnband, _x, _y, minval, NULL);
                                RASTER_DEBUGF(4, "burning pixel (%d, %d) with minval: %f", _x, _y, minval);
                        }
                        else {
                                RASTER_DEBUGF(4, "NOT burning pixel (%d, %d)", _x, _y);
                        }
                        if (status != ES_NONE) {
                                rterror("rt_raster_iterator: Could not set pixel value");
 
                                _rti_iterator_arg_destroy(_param);
                                rt_band_destroy(rtnband);
                                rt_raster_destroy(rtnrast);
 
                                return ES_ERROR;
                        }
                }
        }
 
        /* lots of cleanup */
        _rti_iterator_arg_destroy(_param);
 
        *rtnraster = rtnrast;
        return ES_NONE;
}

References _rti_iterator_arg_callback_clean(), _rti_iterator_arg_callback_init(), _rti_iterator_arg_destroy(), _rti_iterator_arg_empty_init(), _rti_iterator_arg_init(), _rti_iterator_arg_populate(), _rti_iterator_arg_t::arg, _rti_iterator_arg_t::band, rt_raster_t::bands, rt_iterator_arg_t::dst_pixel, _rti_iterator_arg_t::empty, ES_ERROR, ES_NONE, ET_CUSTOM, ET_FIRST, ET_INTERSECTION, ET_LAST, ET_SECOND, ET_UNION, _rti_iterator_arg_t::hasnodata, _rti_iterator_arg_t::isempty, _rti_iterator_arg_t::isnodata, _rti_iterator_arg_t::minval, rt_iterator_t::nbnodata, rt_pixel_t::nodata, rt_iterator_arg_t::nodata, _rti_iterator_arg_t::nodata, _rti_iterator_arg_t::nodataval, rt_raster_t::numBands, _rti_iterator_arg_t::offset, PT_END, rtpixdump::rast, _rti_iterator_arg_t::raster, RASTER_DEBUG, RASTER_DEBUGF, rt_band_destroy(), rt_band_get_min_value(), rt_band_get_nearest_pixel(), rt_band_get_pixel(), rt_band_set_pixel(), rt_pixel_set_to_array(), rt_raster_destroy(), rt_raster_from_two_rasters(), 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_is_empty(), rt_raster_new(), rt_raster_same_alignment(), rt_raster_set_scale(), rtalloc(), _rti_iterator_arg_t::rtband, rtdealloc(), rterror(), rtinfo(), rtrealloc(), rt_iterator_arg_t::src_pixel, rt_pixel_t::value, genraster::value, rt_iterator_arg_t::values, _rti_iterator_arg_t::values, rt_pixel_t::x, pixval::x, rt_pixel_t::y, and pixval::y.

Referenced by RASTER_clip(), RASTER_nMapAlgebra(), RASTER_nMapAlgebraExpr(), RASTER_setPixelValuesGeomval(), RASTER_union_finalfn(), RASTER_union_transfn(), and test_raster_iterator().

Here is the call graph for this function:

Here is the caller graph for this function: