rt_raster_intersects()

rt_errorstate rt_raster_intersects	(	rt_raster	rast1,
		int	nband1,
		rt_raster	rast2,
		int	nband2,
		int *	intersects
	)

Return zero if error occurred in function.

Return ES_ERROR if error occurred in function.

Parameter intersects returns non-zero if two rasters intersect

Parameters

rast1	: the first raster whose band will be tested
nband1	: the 0-based band of raster rast1 to use if value is less than zero, bands are ignored. if nband1 gte zero, nband2 must be gte zero
rast2	: the second raster whose band will be tested
nband2	: the 0-based band of raster rast2 to use if value is less than zero, bands are ignored if nband2 gte zero, nband1 must be gte zero
intersects	: non-zero value if the two rasters' bands intersects

Returns

ES_NONE if success, ES_ERROR if error

Definition at line 988 of file rt_spatial_relationship.c.

  {
        int i;
        int within = 0;
 
        LWGEOM *hull[2] = {NULL};
        GEOSGeometry *ghull[2] = {NULL};
 
        uint16_t width1;
        uint16_t height1;
        uint16_t width2;
        uint16_t height2;
        double area1;
        double area2;
        double pixarea1;
        double pixarea2;
        rt_raster rastS = NULL;
        rt_raster rastL = NULL;
        uint16_t *widthS = NULL;
        uint16_t *heightS = NULL;
        uint16_t *widthL = NULL;
        uint16_t *heightL = NULL;
        int nbandS;
        int nbandL;
        rt_band bandS = NULL;
        rt_band bandL = NULL;
        int hasnodataS = FALSE;
        int hasnodataL = FALSE;
        double nodataS = 0;
        double nodataL = 0;
        int isnodataS = 0;
        int isnodataL = 0;
        double gtS[6] = {0};
        double igtL[6] = {0};
 
        uint32_t row;
        uint32_t rowoffset;
        uint32_t col;
        uint32_t coloffset;
 
        enum line_points {X1, Y1, X2, Y2};
        enum point {pX, pY};
        double lineS[4];
        double Qr[2];
        double valS;
        double valL;
 
        RASTER_DEBUG(3, "Starting");
 
        assert(NULL != rast1);
        assert(NULL != rast2);
        assert(NULL != intersects);
 
        if (nband1 < 0 && nband2 < 0) {
                nband1 = -1;
                nband2 = -1;
        }
        else {
                assert(nband1 >= 0 && nband1 < rt_raster_get_num_bands(rast1));
                assert(nband2 >= 0 && nband2 < rt_raster_get_num_bands(rast2));
        }
 
        /* same srid */
        if (rt_raster_get_srid(rast1) != rt_raster_get_srid(rast2)) {
                rterror("rt_raster_intersects: The two rasters provided have different SRIDs");
                *intersects = 0;
                return ES_ERROR;
        }
 
        /* raster extents need to intersect */
        do {
                int rtn;
 
                initGEOS(rtinfo, lwgeom_geos_error);
 
                rtn = 1;
 
                if ((rt_raster_get_convex_hull(rast1, &(hull[0])) != ES_NONE) || !hull[0]) {
                        break;
                }
                ghull[0] = (GEOSGeometry *) LWGEOM2GEOS(hull[0], 0);
                if (!ghull[0]) {
                        lwgeom_free(hull[0]);
                        break;
                }
 
                if ((rt_raster_get_convex_hull(rast2, &(hull[1])) != ES_NONE) || !hull[1]) {
                        GEOSGeom_destroy(ghull[0]);
                        lwgeom_free(hull[0]);
                        break;
                }
                ghull[1] = (GEOSGeometry *) LWGEOM2GEOS(hull[1], 0);
                if (!ghull[0]) {
                        GEOSGeom_destroy(ghull[0]);
                        lwgeom_free(hull[1]);
                        lwgeom_free(hull[0]);
                        break;
                }
 
                /* test to see if raster within the other */
                within = 0;
                if (GEOSWithin(ghull[0], ghull[1]) == 1)
                        within = -1;
                else if (GEOSWithin(ghull[1], ghull[0]) == 1)
                        within = 1;
 
                if (within != 0)
                        rtn = 1;
                else
                        rtn = GEOSIntersects(ghull[0], ghull[1]);
 
                for (i = 0; i < 2; i++) {
                        GEOSGeom_destroy(ghull[i]);
                        lwgeom_free(hull[i]);
                }
 
                if (rtn != 2) {
                        RASTER_DEBUGF(4, "convex hulls of rasters do %sintersect", rtn != 1 ? "NOT " : "");
                        if (rtn != 1) {
                                *intersects = 0;
                                return ES_NONE;
                        }
                        /* band isn't specified */
                        else if (nband1 < 0) {
                                *intersects = 1;
                                return ES_NONE;
                        }
                }
                else {
                        RASTER_DEBUG(4, "GEOSIntersects() returned a 2!!!!");
                }
        }
        while (0);
 
        /* smaller raster by area or width */
        width1 = rt_raster_get_width(rast1);
        height1 = rt_raster_get_height(rast1);
        width2 = rt_raster_get_width(rast2);
        height2 = rt_raster_get_height(rast2);
        pixarea1 = fabs(rt_raster_get_x_scale(rast1) * rt_raster_get_y_scale(rast1));
        pixarea2 = fabs(rt_raster_get_x_scale(rast2) * rt_raster_get_y_scale(rast2));
        area1 = fabs(width1 * height1 * pixarea1);
        area2 = fabs(width2 * height2 * pixarea2);
        RASTER_DEBUGF(4, "pixarea1, pixarea2, area1, area2 = %f, %f, %f, %f",
                pixarea1, pixarea2, area1, area2);
        if (
                (within <= 0) ||
                (area1 < area2) ||
                FLT_EQ(area1, area2) ||
                (area1 < pixarea2) || /* area of rast1 smaller than pixel area of rast2 */
                FLT_EQ(area1, pixarea2)
        ) {
                rastS = rast1;
                nbandS = nband1;
                widthS = &width1;
                heightS = &height1;
 
                rastL = rast2;
                nbandL = nband2;
                widthL = &width2;
                heightL = &height2;
        }
        else {
                rastS = rast2;
                nbandS = nband2;
                widthS = &width2;
                heightS = &height2;
 
                rastL = rast1;
                nbandL = nband1;
                widthL = &width1;
                heightL = &height1;
        }
 
        /* no band to use, set band to zero */
        if (nband1 < 0) {
                nbandS = 0;
                nbandL = 0;
        }
 
        RASTER_DEBUGF(4, "rast1 @ %p", rast1);
        RASTER_DEBUGF(4, "rast2 @ %p", rast2);
        RASTER_DEBUGF(4, "rastS @ %p", rastS);
        RASTER_DEBUGF(4, "rastL @ %p", rastL);
 
        /* load band of smaller raster */
        bandS = rt_raster_get_band(rastS, nbandS);
        if (NULL == bandS) {
                rterror("rt_raster_intersects: Could not get band %d of the first raster", nbandS);
                *intersects = 0;
                return ES_ERROR;
        }
 
        hasnodataS = rt_band_get_hasnodata_flag(bandS);
        if (hasnodataS != FALSE)
                rt_band_get_nodata(bandS, &nodataS);
 
        /* load band of larger raster */
        bandL = rt_raster_get_band(rastL, nbandL);
        if (NULL == bandL) {
                rterror("rt_raster_intersects: Could not get band %d of the first raster", nbandL);
                *intersects = 0;
                return ES_ERROR;
        }
 
        hasnodataL = rt_band_get_hasnodata_flag(bandL);
        if (hasnodataL != FALSE)
                rt_band_get_nodata(bandL, &nodataL);
 
        /* no band to use, ignore nodata */
        if (nband1 < 0) {
                hasnodataS = FALSE;
                hasnodataL = FALSE;
        }
 
        /* hasnodata(S|L) = TRUE and one of the two bands is isnodata */
        if (
                (hasnodataS && rt_band_get_isnodata_flag(bandS)) ||
                (hasnodataL && rt_band_get_isnodata_flag(bandL))
        ) {
                RASTER_DEBUG(3, "One of the two raster bands is NODATA. The two rasters do not intersect");
                *intersects = 0;
                return ES_NONE;
        }
 
        /* special case where a raster can fit inside another raster's pixel */
        if (within != 0 && ((pixarea1 > area2) || (pixarea2 > area1))) {
                RASTER_DEBUG(4, "Using special case of raster fitting into another raster's pixel");
                /* 3 x 3 search */
                for (coloffset = 0; coloffset < 3; coloffset++) {
                        for (rowoffset = 0; rowoffset < 3; rowoffset++) {
                                for (col = coloffset; col < *widthS; col += 3) {
                                        for (row = rowoffset; row < *heightS; row += 3) {
                                                if (hasnodataS == FALSE)
                                                        valS = 1;
                                                else if (rt_band_get_pixel(bandS, col, row, &valS, &isnodataS) != ES_NONE)
                                                        continue;
 
                                                if ((hasnodataS == FALSE) || !isnodataS) {
                                                        rt_raster_cell_to_geopoint(
                                                                rastS,
                                                                col, row,
                                                                &(lineS[X1]), &(lineS[Y1]),
                                                                gtS
                                                        );
 
                                                        if (rt_raster_geopoint_to_cell(
                                                                rastL,
                                                                lineS[X1], lineS[Y1],
                                                                &(Qr[pX]), &(Qr[pY]),
                                                                igtL
                                                        ) != ES_NONE) {
                                                                continue;
                                                        }
 
                                                        if ((Qr[pX] < 0 || Qr[pX] >= *widthL) ||
                                                            (Qr[pY] < 0 || Qr[pY] >= *heightL))
                                                        {
                                                                continue;
                                                        }
 
                                                        if (hasnodataS == FALSE)
                                                                valL = 1;
                                                        else if (rt_band_get_pixel(bandL, Qr[pX], Qr[pY], &valL, &isnodataL) != ES_NONE)
                                                                continue;
 
                                                        if ((hasnodataL == FALSE) || !isnodataL) {
                                                                RASTER_DEBUG(3, "The two rasters do intersect");
                                                                *intersects = 1;
                                                                return ES_NONE;
                                                        }
                                                }
                                        }
                                }
                        }
                }
                RASTER_DEBUG(4, "Smaller raster not in the other raster's pixel. Continuing");
        }
 
        RASTER_DEBUG(4, "Testing smaller raster vs larger raster");
        *intersects = rt_raster_intersects_algorithm(
                rastS, rastL,
                bandS, bandL,
                hasnodataS, hasnodataL,
                nodataS, nodataL
        );
 
        if (*intersects) return ES_NONE;
 
        RASTER_DEBUG(4, "Testing larger raster vs smaller raster");
        *intersects = rt_raster_intersects_algorithm(
                rastL, rastS,
                bandL, bandS,
                hasnodataL, hasnodataS,
                nodataL, nodataS
        );
 
        if (*intersects) return ES_NONE;
 
        RASTER_DEBUG(3, "The two rasters do not intersect");
 
        *intersects = 0;
        return ES_NONE;
}

References ES_ERROR, ES_NONE, FALSE, FLT_EQ, LWGEOM2GEOS(), lwgeom_free(), lwgeom_geos_error(), RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_band_get_nodata(), rt_band_get_pixel(), rt_raster_cell_to_geopoint(), rt_raster_geopoint_to_cell(), rt_raster_get_band(), rt_raster_get_convex_hull(), rt_raster_get_height(), rt_raster_get_num_bands(), rt_raster_get_srid(), rt_raster_get_width(), rt_raster_get_x_scale(), rt_raster_get_y_scale(), rt_raster_intersects_algorithm(), rterror(), rtinfo(), and within().

Referenced by RASTER_intersects(), and test_raster_intersects().

Here is the call graph for this function:

Here is the caller graph for this function: