/*- rt_band ----------------------------------------------------------*/
-struct rt_extband_t {
- uint8_t bandNum;
- char* path; /* externally owned ? */
-};
-
-struct rt_band_t {
- rt_pixtype pixtype;
- int32_t offline;
- uint16_t width;
- uint16_t height;
- int32_t hasnodata; /* a flag indicating if this band contains nodata values */
- int32_t isnodata; /* a flag indicating if this band is filled only with
- nodata values */
- double nodataval; /* int will be converted ... */
- int32_t ownsData; /* XXX mloskot: its behaviour needs to be documented */
-
- union {
- void* mem; /* actual data, externally owned */
- struct rt_extband_t offline;
- } data;
-
-};
-
rt_band
rt_band_new_inline(uint16_t width, uint16_t height,
rt_pixtype pixtype, uint32_t hasnodata, double nodataval,
int
rt_band_set_nodata(rt_band band, double val) {
rt_pixtype pixtype = PT_END;
- //double oldnodataval = band->nodataval;
+ /*
+ double oldnodataval = band->nodataval;
+ */
int32_t checkvalint = 0;
uint32_t checkvaluint = 0;
RASTER_DEBUGF(3, "rt_band_set_nodata: band->nodataval = %f", band->nodataval);
- // the nodata value was just set, so this band has NODATA
+ /* the nodata value was just set, so this band has NODATA */
rt_band_set_hasnodata_flag(band, 1);
#ifdef POSTGIS_RASTER_WARN_ON_TRUNCATION
return index;
}
-
+/**
+ * Convert an xr, yr raster point to an xw, yw point on map
+ *
+ * @param raster : the raster to get info from
+ * @param xr : the pixel's column
+ * @param yr : the pixel's row
+ * @param xw : output parameter, X ordinate of the geographical point
+ * @param yw : output parameter, Y ordinate of the geographical point
+ */
void
rt_raster_cell_to_geopoint(rt_raster raster,
- double x, double y,
- double* x1, double* y1) {
+ double xr, double yr,
+ double* xw, double* yw
+) {
+ assert(NULL != raster);
+ assert(NULL != xw);
+ assert(NULL != yw);
+ /* Six parameters affine transformation */
+ *xw = raster->scaleX * xr + raster->skewX * yr + raster->ipX;
+ *yw = raster->scaleY * yr + raster->skewY * xr + raster->ipY;
- assert(NULL != raster);
- assert(NULL != x1);
- assert(NULL != y1);
+ RASTER_DEBUGF(3, "rt_raster_cell_to_geopoint(%g, %g) = %g, %g",
+ xr, yr, *xw, *yw);
- /* Six parameters affine transformation */
- *x1 = raster->scaleX * x + raster->skewX * y + raster->ipX;
- *y1 = raster->scaleY * y + raster->skewY * x + raster->ipY;
+ /*
+ {
+ double gt[] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+ double x1;
+ double y1;
+
+ gt[0] = rt_raster_get_x_offset(raster);
+ gt[1] = rt_raster_get_x_scale(raster);
+ gt[2] = rt_raster_get_x_skew(raster);
+ gt[3] = rt_raster_get_y_offset(raster);
+ gt[4] = rt_raster_get_y_skew(raster);
+ gt[5] = rt_raster_get_y_scale(raster);
+
+ GDALApplyGeoTransform(gt, xr, yr, &x1, &y1);
+ RASTER_DEBUGF(4, "GDALApplyGeoTransform for (%f, %f) = (%f, %f)",
+ xr, yr, x1, y1);
+ }
+ */
+}
- RASTER_DEBUGF(3, "rt_raster_cell_to_geopoint(%g,%g)", x, y);
- RASTER_DEBUGF(3, " ipx/y:%g/%g", raster->ipX, raster->ipY);
- RASTER_DEBUGF(3, "cell_to_geopoint: ipX:%g, ipY:%g, %g,%g -> %g,%g",
- raster->ipX, raster->ipY, x, y, *x1, *y1);
+/**
+ * Convert an xw,yw map point to a xr,yr raster point
+ *
+ * @param raster : the raster to get info from
+ * @param xw : X ordinate of the geographical point
+ * @param yw : Y ordinate of the geographical point
+ * @param xr : output parameter, the pixel's column
+ * @param yr : output parameter, the pixel's row
+ *
+ * @return if zero, error occurred in function
+ */
+int
+rt_raster_geopoint_to_cell(rt_raster raster,
+ double xw, double yw,
+ double *xr, double *yr
+) {
+ double d;
+
+ assert(NULL != raster);
+ assert(NULL != xr);
+ assert(NULL != yr);
+
+ /* Six parameters affine transformation */
+ d = raster->skewX * raster->skewY - raster->scaleX * raster->scaleY;
+ if (FLT_EQ(d, 0)) {
+ rterror("rt_raster_geopoint_to_cell: Attempting to compute coordinates on a raster with scale equal to 0");
+ return 0;
+ }
+
+ *xr = (raster->skewX * (yw - raster->ipY) - raster->scaleY * (xw - raster->ipX)) / d;
+ *yr = (raster->scaleX * (yw - raster->ipY) + raster->skewY * (raster->ipX - xw)) / (-1 * d);
+
+ *xr = floor(*xr);
+ *yr = floor(*yr);
+
+ RASTER_DEBUGF(3, "rt_raster_geopoint_to_cell(%g, %g) = %g, %g",
+ xw, yw, *xr, *yr);
+
+ /*
+ {
+ double gt[] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+ double igt[6];
+ double x1;
+ double y1;
+
+ gt[0] = rt_raster_get_x_offset(raster);
+ gt[1] = rt_raster_get_x_scale(raster);
+ gt[2] = rt_raster_get_x_skew(raster);
+ gt[3] = rt_raster_get_y_offset(raster);
+ gt[4] = rt_raster_get_y_skew(raster);
+ gt[5] = rt_raster_get_y_scale(raster);
+ GDALInvGeoTransform(gt, igt);
+
+ GDALApplyGeoTransform(igt, xw, yw, &x1, &y1);
+ RASTER_DEBUGF(4, "GDALApplyGeoTransform for (%f, %f) = (%f, %f)",
+ xw, yw, x1, y1);
+ }
+ */
+
+ return 1;
}
rt_geomval
* postgresql memory context to work with, and the memory created
* for pszSrcText is created outside this context.
**/
- //rtdealloc(pszSrcText);
+ /*rtdealloc(pszSrcText);*/
free(pszSrcText);
pszSrcText = NULL;
/* And we should check if the band is a nodata band */
/* TODO: No!! This is too slow */
- //rt_band_check_is_nodata(band);
+ /*rt_band_check_is_nodata(band);*/
return band;
}
if (band->isnodata) *ptr |= BANDTYPE_FLAG_ISNODATA;
ptr += 1;
-#if 0 // no padding required for WKB
+#if 0 /* no padding required for WKB */
/* Add padding (if needed) */
if (pixbytes > 1) {
memset(ptr, '\0', pixbytes - 1);
return 0;
}
-#if 0 // no padding for WKB
+#if 0 /* no padding for WKB */
/* Consistency checking (ptr is pixbytes-aligned) */
assert(!((uint64_t) ptr % pixbytes));
#endif
}
}
-#if 0 // no padding for WKB
+#if 0 /* no padding for WKB */
/* Pad up to 8-bytes boundary */
while ((uint64_t) ptr % 8) {
*ptr = 0;
/* Align size to 8-bytes boundary (trailing padding) */
/* XXX jorgearevalo: bug here. If the size is actually 8-bytes aligned,
this line will add 8 bytes trailing padding, and it's not necessary */
- //size += 8 - (size % 8);
+ /*size += 8 - (size % 8);*/
if (size % 8)
size += 8 - (size % 8);
* @param grid_yw : the Y value of point on grid to align raster to
* @param skew_x : the X skew of the raster
* @param skew_y : the Y skew of the raster
+ * @param options : array of options. only option is "ALL_TOUCHED"
*
* @return the raster of the provided geometry
*/
double *scale_x, double *scale_y,
double *ul_xw, double *ul_yw,
double *grid_xw, double *grid_yw,
- double *skew_x, double *skew_y
+ double *skew_x, double *skew_y,
+ char **options
) {
rt_raster rast;
int i = 0;
OGRSpatialReferenceH src_sr = NULL;
OGRGeometryH src_geom;
OGREnvelope src_env;
+ OGRwkbGeometryType wkbtype = wkbUnknown;
int ul_user = 0;
double djunk = 0;
RASTER_DEBUGF(3, "scale (x, y) = %f, %f", _scale_x, _scale_y);
RASTER_DEBUGF(3, "dim (x, y) = %d, %d", _width, _height);
- if (
- (wkbFlatten(OGR_G_GetGeometryType(src_geom)) == wkbPoint) &&
- FLT_NEQ(_scale_x, 0) &&
- FLT_NEQ(_scale_y, 0)
+ /*
+ if geometry is a point or point-like and bounds not set, increase
+ extent by half-pixel to avoid missing points on border
+ */
+ wkbtype = wkbFlatten(OGR_G_GetGeometryType(src_geom));
+ if ((
+ (wkbtype == wkbPoint) ||
+ (wkbtype == wkbMultiPoint)
+ ) &&
+ FLT_EQ(_width, 0) &&
+ FLT_EQ(_height, 0)
) {
src_env.MinX -= (_scale_x / 2.);
src_env.MaxX += (_scale_x / 2.);
1, &src_geom,
NULL, NULL,
_value,
- NULL,
+ options,
NULL, NULL
);
rtdealloc(band_list);
return rast;
}
+
+static
+int rt_raster_intersects_algorithm(
+ rt_raster rast1, rt_raster rast2,
+ rt_band band1, rt_band band2,
+ int hasnodata1, int hasnodata2,
+ double nodata1, double nodata2
+) {
+ int i;
+ int byHeight = 1;
+ uint32_t dimValue;
+
+ uint32_t row;
+ uint32_t rowoffset;
+ uint32_t col;
+ uint32_t coloffset;
+
+ enum line_points {X1, Y1, X2, Y2};
+ enum point {pX, pY};
+ double line1[4];
+ double line2[4];
+ double P[2];
+ double Qw[2];
+ double Qr[2];
+ double d;
+ double val1;
+ int noval1;
+ double val2;
+ int noval2;
+ uint32_t adjacent[8] = {0};
+
+ double xscale;
+ double yscale;
+
+ uint16_t width1;
+ uint16_t height1;
+ uint16_t width2;
+ uint16_t height2;
+
+ width1 = rt_raster_get_width(rast1);
+ height1 = rt_raster_get_height(rast1);
+ width2 = rt_raster_get_width(rast2);
+ height2 = rt_raster_get_height(rast2);
+
+ /* sampling scale */
+ xscale = fmin(rt_raster_get_x_scale(rast1), rt_raster_get_x_scale(rast2)) / 10.;
+ yscale = fmin(rt_raster_get_y_scale(rast1), rt_raster_get_y_scale(rast2)) / 10.;
+
+ /* see if skew made rast2's rows are parallel to rast1's cols */
+ rt_raster_cell_to_geopoint(
+ rast1,
+ 0, 0,
+ &(line1[X1]), &(line1[Y1])
+ );
+
+ rt_raster_cell_to_geopoint(
+ rast1,
+ 0, height1,
+ &(line1[X2]), &(line1[Y2])
+ );
+
+ rt_raster_cell_to_geopoint(
+ rast2,
+ 0, 0,
+ &(line2[X1]), &(line2[Y1])
+ );
+
+ rt_raster_cell_to_geopoint(
+ rast2,
+ width2, 0,
+ &(line2[X2]), &(line2[Y2])
+ );
+
+ /* parallel vertically */
+ if (FLT_EQ(line1[X2] - line1[X1], 0.) && FLT_EQ(line2[X2] - line2[X1], 0.))
+ byHeight = 0;
+ /* parallel */
+ else if (FLT_EQ(((line1[Y2] - line1[Y1]) / (line1[X2] - line1[X1])), ((line2[Y2] - line2[Y1]) / (line2[X2] - line2[X1]))))
+ byHeight = 0;
+
+ if (byHeight)
+ dimValue = height2;
+ else
+ dimValue = width2;
+ RASTER_DEBUGF(4, "byHeight: %d, dimValue: %d", byHeight, dimValue);
+
+ /* 3 x 3 search */
+ for (coloffset = 0; coloffset < 3; coloffset++) {
+ for (rowoffset = 0; rowoffset < 3; rowoffset++) {
+ /* smaller raster */
+ for (col = coloffset; col <= width1; col += 3) {
+
+ rt_raster_cell_to_geopoint(
+ rast1,
+ col, 0,
+ &(line1[X1]), &(line1[Y1])
+ );
+
+ rt_raster_cell_to_geopoint(
+ rast1,
+ col, height1,
+ &(line1[X2]), &(line1[Y2])
+ );
+
+ /* larger raster */
+ for (row = rowoffset; row <= dimValue; row += 3) {
+
+ if (byHeight) {
+ rt_raster_cell_to_geopoint(
+ rast2,
+ 0, row,
+ &(line2[X1]), &(line2[Y1])
+ );
+
+ rt_raster_cell_to_geopoint(
+ rast2,
+ width2, row,
+ &(line2[X2]), &(line2[Y2])
+ );
+ }
+ else {
+ rt_raster_cell_to_geopoint(
+ rast2,
+ row, 0,
+ &(line2[X1]), &(line2[Y1])
+ );
+
+ rt_raster_cell_to_geopoint(
+ rast2,
+ row, height2,
+ &(line2[X2]), &(line2[Y2])
+ );
+ }
+
+ RASTER_DEBUGF(4, "(col, row) = (%d, %d)", col, row);
+ RASTER_DEBUGF(4, "line1(x1, y1, x2, y2) = (%f, %f, %f, %f)",
+ line1[X1], line1[Y1], line1[X2], line1[Y2]);
+ RASTER_DEBUGF(4, "line2(x1, y1, x2, y2) = (%f, %f, %f, %f)",
+ line2[X1], line2[Y1], line2[X2], line2[Y2]);
+
+ /* intersection */
+ /* http://en.wikipedia.org/wiki/Line-line_intersection */
+ d = ((line1[X1] - line1[X2]) * (line2[Y1] - line2[Y2])) - ((line1[Y1] - line1[Y2]) * (line2[X1] - line2[X2]));
+ /* no intersection */
+ if (FLT_EQ(d, 0.)) {
+ continue;
+ }
+
+ P[pX] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[X1] - line2[X2])) -
+ ((line1[X1] - line1[X2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
+ P[pX] = P[pX] / d;
+
+ P[pY] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[Y1] - line2[Y2])) -
+ ((line1[Y1] - line1[Y2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
+ P[pY] = P[pY] / d;
+
+ RASTER_DEBUGF(4, "P(x, y) = (%f, %f)", P[pX], P[pY]);
+
+ /* intersection within bounds */
+ if ((
+ (FLT_EQ(P[pX], line1[X1]) || FLT_EQ(P[pX], line1[X2])) ||
+ (P[pX] > fmin(line1[X1], line1[X2]) && (P[pX] < fmax(line1[X1], line1[X2])))
+ ) && (
+ (FLT_EQ(P[pY], line1[Y1]) || FLT_EQ(P[pY], line1[Y2])) ||
+ (P[pY] > fmin(line1[Y1], line1[Y2]) && (P[pY] < fmax(line1[Y1], line1[Y2])))
+ ) && (
+ (FLT_EQ(P[pX], line2[X1]) || FLT_EQ(P[pX], line2[X2])) ||
+ (P[pX] > fmin(line2[X1], line2[X2]) && (P[pX] < fmax(line2[X1], line2[X2])))
+ ) && (
+ (FLT_EQ(P[pY], line2[Y1]) || FLT_EQ(P[pY], line2[Y2])) ||
+ (P[pY] > fmin(line2[Y1], line2[Y2]) && (P[pY] < fmax(line2[Y1], line2[Y2])))
+ )) {
+ RASTER_DEBUG(4, "within bounds");
+
+ for (i = 0; i < 8; i++) adjacent[i] = 0;
+
+ /* test points around intersection */
+ for (i = 0; i < 8; i++) {
+ switch (i) {
+ case 7:
+ Qw[pX] = P[pX] - xscale;
+ Qw[pY] = P[pY] + yscale;
+ break;
+ /* 270 degrees = 09:00 */
+ case 6:
+ Qw[pX] = P[pX] - xscale;
+ Qw[pY] = P[pY];
+ break;
+ case 5:
+ Qw[pX] = P[pX] - xscale;
+ Qw[pY] = P[pY] - yscale;
+ break;
+ /* 180 degrees = 06:00 */
+ case 4:
+ Qw[pX] = P[pX];
+ Qw[pY] = P[pY] - yscale;
+ break;
+ case 3:
+ Qw[pX] = P[pX] + xscale;
+ Qw[pY] = P[pY] - yscale;
+ break;
+ /* 90 degrees = 03:00 */
+ case 2:
+ Qw[pX] = P[pX] + xscale;
+ Qw[pY] = P[pY];
+ break;
+ /* 45 degrees */
+ case 1:
+ Qw[pX] = P[pX] + xscale;
+ Qw[pY] = P[pY] + yscale;
+ break;
+ /* 0 degrees = 00:00 */
+ case 0:
+ Qw[pX] = P[pX];
+ Qw[pY] = P[pY] + yscale;
+ break;
+ }
+
+ /* unable to convert point to cell */
+ noval1 = 0;
+ if (!rt_raster_geopoint_to_cell(
+ rast1,
+ Qw[pX], Qw[pY],
+ &(Qr[pX]), &(Qr[pY])
+ )) {
+ noval1 = 1;
+ }
+ /* cell is outside bounds of grid */
+ else if (
+ (Qr[pX] < 0 || Qr[pX] > width1 || FLT_EQ(Qr[pX], width1)) ||
+ (Qr[pY] < 0 || Qr[pY] > height1 || FLT_EQ(Qr[pY], height1))
+ ) {
+ noval1 = 1;
+ }
+ else if (hasnodata1 == FALSE)
+ val1 = 1;
+ /* unable to get value at cell */
+ else if (rt_band_get_pixel(band1, Qr[pX], Qr[pY], &val1) < 0)
+ noval1 = 1;
+
+ /* unable to convert point to cell */
+ noval2 = 0;
+ if (!rt_raster_geopoint_to_cell(
+ rast2,
+ Qw[pX], Qw[pY],
+ &(Qr[pX]), &(Qr[pY])
+ )) {
+ noval2 = 1;
+ }
+ /* cell is outside bounds of grid */
+ else if (
+ (Qr[pX] < 0 || Qr[pX] > width2 || FLT_EQ(Qr[pX], width2)) ||
+ (Qr[pY] < 0 || Qr[pY] > height2 || FLT_EQ(Qr[pY], height2))
+ ) {
+ noval2 = 1;
+ }
+ else if (hasnodata2 == FALSE)
+ val2 = 1;
+ /* unable to get value at cell */
+ else if (rt_band_get_pixel(band2, Qr[pX], Qr[pY], &val2) < 0)
+ noval2 = 1;
+
+ if (!noval1) {
+ RASTER_DEBUGF(4, "val1 = %f", val1);
+ }
+ if (!noval2) {
+ RASTER_DEBUGF(4, "val2 = %f", val2);
+ }
+
+ /* pixels touch */
+ if (!noval1 && (
+ (hasnodata1 == FALSE) || (
+ (hasnodata1 != FALSE) &&
+ FLT_NEQ(val1, nodata1)
+ )
+ )) {
+ adjacent[i]++;
+ }
+ if (!noval2 && (
+ (hasnodata2 == FALSE) || (
+ (hasnodata2 != FALSE) &&
+ FLT_NEQ(val2, nodata2)
+ )
+ )) {
+ adjacent[i] += 3;
+ }
+
+ /* two pixel values not present */
+ if (noval1 || noval2) {
+ RASTER_DEBUGF(4, "noval1 = %d, noval2 = %d", noval1, noval2);
+ continue;
+ }
+
+ /* pixels valid, so intersect */
+ if ((
+ (hasnodata1 == FALSE) || (
+ (hasnodata1 != FALSE) &&
+ FLT_NEQ(val1, nodata1)
+ )
+ ) && (
+ (hasnodata2 == FALSE) || (
+ (hasnodata2 != FALSE) &&
+ FLT_NEQ(val2, nodata2)
+ )
+ )) {
+ RASTER_DEBUG(3, "The two rasters do intersect");
+
+ return 1;
+ }
+ }
+
+ /* pixels touch */
+ for (i = 0; i < 4; i++) {
+ RASTER_DEBUGF(4, "adjacent[%d] = %d, adjacent[%d] = %d"
+ , i, adjacent[i], i + 4, adjacent[i + 4]);
+ if (adjacent[i] == 0) continue;
+
+ if (adjacent[i] + adjacent[i + 4] == 4) {
+ RASTER_DEBUG(3, "The two rasters touch");
+
+ return 1;
+ }
+ }
+ }
+ else {
+ RASTER_DEBUG(4, "outside of bounds");
+ }
+ }
+ }
+ }
+ }
+
+ return 0;
+}
+
+/**
+ * Return zero if error occurred in function.
+ * Parameter intersects returns non-zero if two rasters intersect
+ *
+ * @param rast1 : the first raster whose band will be tested
+ * @param 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
+ * @param rast2 : the second raster whose band will be tested
+ * @param 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
+ * @param intersects : non-zero value if the two rasters' bands intersects
+ *
+ * @return if zero, an error occurred in function
+ */
+int
+rt_raster_intersects(
+ rt_raster rast1, int nband1,
+ rt_raster rast2, int nband2,
+ int *intersects
+) {
+ int i;
+ int j;
+ int within = 0;
+
+ LWPOLY *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;
+
+ 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);
+
+ 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 0;
+ }
+
+ /* raster extents need to intersect */
+ do {
+ int rtn;
+
+ initGEOS(lwnotice, lwgeom_geos_error);
+
+ rtn = 1;
+ for (i = 0; i < 2; i++) {
+ hull[i] = rt_raster_get_convex_hull(i < 1 ? rast1 : rast2);
+ if (NULL == hull[i]) {
+ for (j = 0; j < i; j++) {
+ GEOSGeom_destroy(ghull[j]);
+ lwpoly_free(hull[j]);
+ }
+ rtn = 0;
+ break;
+ }
+ ghull[i] = (GEOSGeometry *) LWGEOM2GEOS((LWGEOM *) hull[i]);
+ if (NULL == ghull[i]) {
+ for (j = 0; j < i; j++) {
+ GEOSGeom_destroy(ghull[j]);
+ lwpoly_free(hull[j]);
+ }
+ lwpoly_free(hull[i]);
+ rtn = 0;
+ break;
+ }
+ }
+ if (!rtn) 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]);
+ lwpoly_free(hull[i]);
+ }
+
+ if (rtn != 2) {
+ RASTER_DEBUGF(4, "convex hulls of rasters do %sintersect", rtn != 1 ? "NOT " : "");
+ if (rtn != 1) {
+ *intersects = 0;
+ return 1;
+ }
+ /* band isn't specified */
+ else if (nband1 < 0) {
+ *intersects = 1;
+ return 1;
+ }
+ }
+ }
+ 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: Unable to get band %d of the first raster\n", nbandS);
+ *intersects = 0;
+ return 0;
+ }
+ if (bandS->offline) {
+ rterror("rt_raster_intersects not implemented yet for OFFDB bands");
+ *intersects = 0;
+ return 0;
+ }
+
+ hasnodataS = rt_band_get_hasnodata_flag(bandS);
+ if (hasnodataS != FALSE)
+ nodataS = rt_band_get_nodata(bandS);
+
+ /* load band of larger raster */
+ bandL = rt_raster_get_band(rastL, nbandL);
+ if (NULL == bandL) {
+ rterror("rt_raster_intersects: Unable to get band %d of the first raster\n", nbandL);
+ *intersects = 0;
+ return 0;
+ }
+ if (bandL->offline) {
+ rterror("rt_raster_intersects not implemented yet for OFFDB bands");
+ *intersects = 0;
+ return 0;
+ }
+
+ hasnodataL = rt_band_get_hasnodata_flag(bandL);
+ if (hasnodataL != FALSE)
+ nodataL = rt_band_get_nodata(bandL);
+
+ /* no band to use, ignore nodata */
+ if (nband1 < 0) {
+ hasnodataS = FALSE;
+ hasnodataL = FALSE;
+ }
+
+ /* 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) < 0)
+ continue;
+
+ if ((hasnodataS == FALSE) || (
+ (hasnodataS != FALSE) &&
+ FLT_NEQ(valS, nodataS)
+ )) {
+ rt_raster_cell_to_geopoint(
+ rastS,
+ col, row,
+ &(lineS[X1]), &(lineS[Y1])
+ );
+
+ if (!rt_raster_geopoint_to_cell(
+ rastL,
+ lineS[X1], lineS[Y1],
+ &(Qr[pX]), &(Qr[pY])
+ )) {
+ continue;
+ }
+
+ if (
+ (Qr[pX] < 0 || Qr[pX] > *widthL || FLT_EQ(Qr[pX], *widthL)) ||
+ (Qr[pY] < 0 || Qr[pY] > *heightL || FLT_EQ(Qr[pY], *heightL))
+ ) {
+ continue;
+ }
+
+ if (hasnodataS == FALSE)
+ valL = 1;
+ else if (rt_band_get_pixel(bandL, Qr[pX], Qr[pY], &valL) < 0)
+ continue;
+
+ if ((hasnodataL == FALSE) || (
+ (hasnodataL != FALSE) &&
+ FLT_NEQ(valL, nodataL)
+ )) {
+ RASTER_DEBUG(3, "The two rasters do intersect");
+ *intersects = 1;
+ return 1;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ *intersects = rt_raster_intersects_algorithm(
+ rastS, rastL,
+ bandS, bandL,
+ hasnodataS, hasnodataL,
+ nodataS, nodataL
+ );
+
+ if (*intersects) return 1;
+
+ *intersects = rt_raster_intersects_algorithm(
+ rastL, rastS,
+ bandL, bandS,
+ hasnodataL, hasnodataS,
+ nodataL, nodataS
+ );
+
+ if (*intersects) return 1;
+
+ RASTER_DEBUG(3, "The two rasters do not intersect");
+
+ *intersects = 0;
+ return 1;
+}
#endif
#endif
+/* Microsoft does not support fmax or fmin in math.h */
+#if defined(WIN32)
+#if !defined(HAVE_MAX)
+#define fmax max
+#endif
+#if !defined(HAVE_MIN)
+#define fmin min
+#endif
+#endif
+
#include <stdlib.h> /* For size_t, srand and rand */
#include <stdint.h> /* For C99 int types */
#include "liblwgeom.h"
+#include "lwgeom_geos.h"
#include "gdal_alg.h"
#include "gdal_frmts.h"
typedef struct rt_gdaldriver_t* rt_gdaldriver;
typedef struct rt_reclassexpr_t* rt_reclassexpr;
-
-/**
- * Struct definitions
- *
- * These structs are defined here as they are needed elsewhere
- * including rt_pg/rt_pg.c and reduce duplicative declarations
- *
- */
-struct rt_raster_serialized_t {
- /*---[ 8 byte boundary ]---{ */
- uint32_t size; /* required by postgresql: 4 bytes */
- uint16_t version; /* format version (this is version 0): 2 bytes */
- uint16_t numBands; /* Number of bands: 2 bytes */
-
- /* }---[ 8 byte boundary ]---{ */
- double scaleX; /* pixel width: 8 bytes */
-
- /* }---[ 8 byte boundary ]---{ */
- double scaleY; /* pixel height: 8 bytes */
-
- /* }---[ 8 byte boundary ]---{ */
- double ipX; /* insertion point X: 8 bytes */
-
- /* }---[ 8 byte boundary ]---{ */
- double ipY; /* insertion point Y: 8 bytes */
-
- /* }---[ 8 byte boundary ]---{ */
- double skewX; /* skew about the X axis: 8 bytes */
-
- /* }---[ 8 byte boundary ]---{ */
- double skewY; /* skew about the Y axis: 8 bytes */
-
- /* }---[ 8 byte boundary ]--- */
- int32_t srid; /* Spatial reference id: 4 bytes */
- uint16_t width; /* pixel columns: 2 bytes */
- uint16_t height; /* pixel rows: 2 bytes */
-};
-
-/* NOTE: the initial part of this structure matches the layout
- * of data in the serialized form version 0, starting
- * from the numBands element
- */
-struct rt_raster_t {
- uint32_t size;
- uint16_t version;
-
- /* Number of bands, all share the same dimension
- * and georeference */
- uint16_t numBands;
-
- /* Georeference (in projection units) */
- double scaleX; /* pixel width */
- double scaleY; /* pixel height */
- double ipX; /* geo x ordinate of the corner of upper-left pixel */
- double ipY; /* geo y ordinate of the corner of bottom-right pixel */
- double skewX; /* skew about the X axis*/
- double skewY; /* skew about the Y axis */
-
- int32_t srid; /* spatial reference id */
- uint16_t width; /* pixel columns - max 65535 */
- uint16_t height; /* pixel rows - max 65535 */
- rt_band *bands; /* actual bands */
-
-};
-
-/* WKT string representing each polygon in WKT format acompagned by its
-correspoding value */
-struct rt_geomval_t {
- int srid;
- double val;
- char * geom;
-};
-
-/* summary stats of specified band */
-struct rt_bandstats_t {
- double sample;
- uint32_t count;
-
- double min;
- double max;
- double sum;
- double mean;
- double stddev;
-
- double *values;
- int sorted; /* flag indicating that values is sorted ascending by value */
-};
-
-/* histogram bin(s) of specified band */
-struct rt_histogram_t {
- uint32_t count;
- double percent;
-
- double min;
- double max;
-
- int inc_min;
- int inc_max;
-};
-
-/* quantile(s) of the specified band */
-struct rt_quantile_t {
- double quantile;
- double value;
-};
-
-/* listed-list structures for rt_band_get_quantiles_stream */
-struct quantile_llist {
- uint8_t algeq; /* AL-GEQ (1) or AL-GT (0) */
- double quantile;
- uint64_t tau; /* position in sequence */
-
- struct quantile_llist_element *head; /* H index 0 */
- struct quantile_llist_element *tail; /* H index last */
- uint32_t count; /* # of elements in H */
-
- /* faster access to elements at specific intervals */
- struct quantile_llist_index *index;
- uint32_t index_max; /* max # of elements in index */
-
- uint64_t sum1; /* N1H */
- uint64_t sum2; /* N2H */
-};
-
-struct quantile_llist_element {
- double value;
- uint32_t count;
-
- struct quantile_llist_element *prev;
- struct quantile_llist_element *next;
-};
-
-struct quantile_llist_index {
- struct quantile_llist_element *element;
- uint32_t index;
-};
-
-/* number of times a value occurs */
-struct rt_valuecount_t {
- double value;
- uint32_t count;
- double percent;
-};
-
-/* reclassification expression */
-struct rt_reclassexpr_t {
- struct rt_reclassrange {
- double min;
- double max;
- int inc_min; /* include min */
- int inc_max; /* include max */
- int exc_min; /* exceed min */
- int exc_max; /* exceed max */
- } src, dst;
-};
-
-/* gdal driver information */
-struct rt_gdaldriver_t {
- int idx;
- char *short_name;
- char *long_name;
- char *create_options;
-};
-
/**
* Global functions for memory/logging handlers.
*/
rt_quantile rt_band_get_quantiles(rt_bandstats stats,
double *quantiles, int quantiles_count, uint32_t *rtn_count);
+struct quantile_llist;
int quantile_llist_destroy(struct quantile_llist **list,
uint32_t list_count);
int32_t rt_raster_get_srid(rt_raster raster);
/**
- * Convert an x,y raster point to an x1,y1 point on map
+ * Convert an xr, yr raster point to an xw, yw point on map
*
* @param raster : the raster to get info from
- * @param x : the pixel's column
- * @param y : the pixel's row
- * @param x1 : output parameter, X ordinate of the geographical point
- * @param y1 : output parameter, Y ordinate of the geographical point
+ * @param xr : the pixel's column
+ * @param yr : the pixel's row
+ * @param xw : output parameter, X ordinate of the geographical point
+ * @param yw : output parameter, Y ordinate of the geographical point
*/
void rt_raster_cell_to_geopoint(rt_raster raster,
- double x, double y,
- double* x1, double* y1);
+ double xr, double yr,
+ double* xw, double* yw);
+
+/**
+ * Convert an xw, yw map point to a xr, yr raster point
+ *
+ * @param raster : the raster to get info from
+ * @param xw : X ordinate of the geographical point
+ * @param yw : Y ordinate of the geographical point
+ * @param xr : output parameter, the pixel's column
+ * @param yr : output parameter, the pixel's row
+ *
+ * @return if zero, error occurred in function
+ */
+int rt_raster_geopoint_to_cell(rt_raster raster,
+ double xw, double yw,
+ double *xr, double *yr);
/**
* Get raster's polygon convex hull.
* @param grid_yw : the Y value of point on grid to align raster to
* @param skew_x : the X skew of the raster
* @param skew_y : the Y skew of the raster
+ * @param options : array of options. only option is "ALL_TOUCHED"
*
* @return the raster of the provided geometry
*/
double *scale_x, double *scale_y,
double *ul_xw, double *ul_yw,
double *grid_xw, double *grid_yw,
- double *skew_x, double *skew_y);
+ double *skew_x, double *skew_y,
+ char **options
+);
+
+/**
+ * Return zero if error occurred in function.
+ * Parameter intersects returns non-zero if two rasters intersect
+ *
+ * @param rast1 : the first raster whose band will be tested
+ * @param 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
+ * @param rast2 : the second raster whose band will be tested
+ * @param 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
+ * @param intersects : non-zero value if the two rasters' bands intersects
+ *
+ * @return if zero, an error occurred in function
+ */
+int rt_raster_intersects(
+ rt_raster rast1, int nband1,
+ rt_raster rast2, int nband2,
+ int *intersects
+);
/*- utilities -------------------------------------------------------*/
*/
#define ROUND(x, y) (((x > 0.0) ? floor((x * pow(10, y) + 0.5)) : ceil((x * pow(10, y) - 0.5))) / pow(10, y));
+/**
+ * Struct definitions
+ *
+ * These structs are defined here as they are needed elsewhere
+ * including rt_pg/rt_pg.c and reduce duplicative declarations
+ *
+ */
+struct rt_raster_serialized_t {
+ /*---[ 8 byte boundary ]---{ */
+ uint32_t size; /* required by postgresql: 4 bytes */
+ uint16_t version; /* format version (this is version 0): 2 bytes */
+ uint16_t numBands; /* Number of bands: 2 bytes */
+
+ /* }---[ 8 byte boundary ]---{ */
+ double scaleX; /* pixel width: 8 bytes */
+
+ /* }---[ 8 byte boundary ]---{ */
+ double scaleY; /* pixel height: 8 bytes */
+
+ /* }---[ 8 byte boundary ]---{ */
+ double ipX; /* insertion point X: 8 bytes */
+
+ /* }---[ 8 byte boundary ]---{ */
+ double ipY; /* insertion point Y: 8 bytes */
+
+ /* }---[ 8 byte boundary ]---{ */
+ double skewX; /* skew about the X axis: 8 bytes */
+
+ /* }---[ 8 byte boundary ]---{ */
+ double skewY; /* skew about the Y axis: 8 bytes */
+
+ /* }---[ 8 byte boundary ]--- */
+ int32_t srid; /* Spatial reference id: 4 bytes */
+ uint16_t width; /* pixel columns: 2 bytes */
+ uint16_t height; /* pixel rows: 2 bytes */
+};
+
+/* NOTE: the initial part of this structure matches the layout
+ * of data in the serialized form version 0, starting
+ * from the numBands element
+ */
+struct rt_raster_t {
+ uint32_t size;
+ uint16_t version;
+
+ /* Number of bands, all share the same dimension
+ * and georeference */
+ uint16_t numBands;
+
+ /* Georeference (in projection units) */
+ double scaleX; /* pixel width */
+ double scaleY; /* pixel height */
+ double ipX; /* geo x ordinate of the corner of upper-left pixel */
+ double ipY; /* geo y ordinate of the corner of bottom-right pixel */
+ double skewX; /* skew about the X axis*/
+ double skewY; /* skew about the Y axis */
+
+ int32_t srid; /* spatial reference id */
+ uint16_t width; /* pixel columns - max 65535 */
+ uint16_t height; /* pixel rows - max 65535 */
+ rt_band *bands; /* actual bands */
+
+};
+
+struct rt_extband_t {
+ uint8_t bandNum;
+ char* path; /* externally owned ? */
+};
+
+struct rt_band_t {
+ rt_pixtype pixtype;
+ int32_t offline;
+ uint16_t width;
+ uint16_t height;
+ int32_t hasnodata; /* a flag indicating if this band contains nodata values */
+ int32_t isnodata; /* a flag indicating if this band is filled only with
+ nodata values */
+ double nodataval; /* int will be converted ... */
+ int32_t ownsData; /* XXX mloskot: its behaviour needs to be documented */
+
+ union {
+ void* mem; /* actual data, externally owned */
+ struct rt_extband_t offline;
+ } data;
+
+};
+
+/* WKT string representing each polygon in WKT format accompanied by its
+corresponding value */
+struct rt_geomval_t {
+ int srid;
+ double val;
+ char * geom;
+};
+
+/* summary stats of specified band */
+struct rt_bandstats_t {
+ double sample;
+ uint32_t count;
+
+ double min;
+ double max;
+ double sum;
+ double mean;
+ double stddev;
+
+ double *values;
+ int sorted; /* flag indicating that values is sorted ascending by value */
+};
+
+/* histogram bin(s) of specified band */
+struct rt_histogram_t {
+ uint32_t count;
+ double percent;
+
+ double min;
+ double max;
+
+ int inc_min;
+ int inc_max;
+};
+
+/* quantile(s) of the specified band */
+struct rt_quantile_t {
+ double quantile;
+ double value;
+};
+
+/* listed-list structures for rt_band_get_quantiles_stream */
+struct quantile_llist {
+ uint8_t algeq; /* AL-GEQ (1) or AL-GT (0) */
+ double quantile;
+ uint64_t tau; /* position in sequence */
+
+ struct quantile_llist_element *head; /* H index 0 */
+ struct quantile_llist_element *tail; /* H index last */
+ uint32_t count; /* # of elements in H */
+
+ /* faster access to elements at specific intervals */
+ struct quantile_llist_index *index;
+ uint32_t index_max; /* max # of elements in index */
+
+ uint64_t sum1; /* N1H */
+ uint64_t sum2; /* N2H */
+};
+
+struct quantile_llist_element {
+ double value;
+ uint32_t count;
+
+ struct quantile_llist_element *prev;
+ struct quantile_llist_element *next;
+};
+
+struct quantile_llist_index {
+ struct quantile_llist_element *element;
+ uint32_t index;
+};
+
+/* number of times a value occurs */
+struct rt_valuecount_t {
+ double value;
+ uint32_t count;
+ double percent;
+};
+
+/* reclassification expression */
+struct rt_reclassexpr_t {
+ struct rt_reclassrange {
+ double min;
+ double max;
+ int inc_min; /* include min */
+ int inc_max; /* include max */
+ int exc_min; /* exceed min */
+ int exc_max; /* exceed max */
+ } src, dst;
+};
+
+/* gdal driver information */
+struct rt_gdaldriver_t {
+ int idx;
+ char *short_name;
+ char *long_name;
+ char *create_options;
+};
+
#endif /* RT_API_H_INCLUDED */
/* get raster band's meta data */
Datum RASTER_bandmetadata(PG_FUNCTION_ARGS);
+/* determine if two rasters intersect */
+Datum RASTER_intersects(PG_FUNCTION_ARGS);
+
/* Replace function taken from
* http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3
*/
double *skew_x = NULL;
double *skew_y = NULL;
+ char **options = NULL;
+ int options_len = 0;
+
uint32_t num_bands = 0;
int srid = SRID_UNKNOWN;
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: skew (x, y) = %f, %f", skew[0], skew[1]);
+ /* all touched */
+ if (!PG_ARGISNULL(14) && PG_GETARG_BOOL(14) == TRUE) {
+ if (options_len < 1) {
+ options_len = 1;
+ options = (char **) palloc(sizeof(char *) * options_len);
+ }
+ else {
+ options_len++;
+ options = (char **) repalloc(options, sizeof(char *) * options_len);
+ }
+
+ options[options_len - 1] = palloc(sizeof(char*) * (strlen("ALL_TOUCHED=TRUE") + 1));
+ options[options_len - 1] = "ALL_TOUCHED=TRUE";
+ }
+
+ if (options_len) {
+ options_len++;
+ options = (char **) repalloc(options, sizeof(char *) * options_len);
+ options[options_len - 1] = NULL;
+ }
+
/* get geometry's srid */
#ifdef GSERIALIZED_ON
srid = gserialized_get_srid(pggeom);
pfree(hasnodatas);
pfree(nodatavals);
}
+ if (options_len) pfree(options);
lwgeom_free(geom);
PG_FREE_IF_COPY(pggeom, 0);
scale_x, scale_y,
ul_xw, ul_yw,
grid_xw, grid_yw,
- skew_x, skew_y
+ skew_x, skew_y,
+ options
);
if (pixtypes_len) pfree(pixtypes);
pfree(hasnodatas);
pfree(nodatavals);
}
+ if (options_len) pfree(options);
if (!rast) {
- elog(ERROR, "RASTER_asRaster: Could not create rasterize geometry");
+ elog(ERROR, "RASTER_asRaster: Could not rasterize geometry");
PG_RETURN_NULL();
}
PG_RETURN_DATUM(result);
}
+/**
+ * See if two rasters intersect
+ */
+PG_FUNCTION_INFO_V1(RASTER_intersects);
+Datum RASTER_intersects(PG_FUNCTION_ARGS)
+{
+ const int set_count = 2;
+ rt_pgraster *pgrast;
+ rt_raster rast[2] = {NULL};
+ uint32_t bandindex[2] = {0};
+ uint32_t hasbandindex[2] = {0};
+
+ uint32_t i;
+ uint32_t j;
+ uint32_t k;
+ uint32_t numBands;
+ int rtn;
+ int intersects;
+
+ LWPOLY *hull[2] = {NULL};
+ GEOSGeometry *ghull[2] = {NULL};
+
+ for (i = 0, j = 0; i < set_count; i++) {
+ /* pgrast is null, return null */
+ if (PG_ARGISNULL(j)) {
+ for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+ PG_RETURN_NULL();
+ }
+ pgrast = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
+ j++;
+
+ /* raster */
+ rast[i] = rt_raster_deserialize(pgrast, FALSE);
+ if (!rast[i]) {
+ elog(ERROR, "RASTER_intersects: Could not deserialize the %s raster", i < 1 ? "first" : "second");
+ for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+ PG_RETURN_NULL();
+ }
+
+ /* numbands */
+ numBands = rt_raster_get_num_bands(rast[i]);
+ if (numBands < 1) {
+ elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
+ for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+ PG_RETURN_NULL();
+ }
+
+ /* band index */
+ if (!PG_ARGISNULL(j)) {
+ bandindex[i] = PG_GETARG_INT32(j);
+ if (bandindex[i] < 1 || bandindex[i] > numBands) {
+ elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
+ for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+ PG_RETURN_NULL();
+ }
+ hasbandindex[i] = 1;
+ }
+ else
+ hasbandindex[i] = 0;
+ POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
+ POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
+ j++;
+ }
+
+ /* hasbandindex must be balanced */
+ if (
+ (hasbandindex[0] && !hasbandindex[1]) ||
+ (!hasbandindex[0] && hasbandindex[1])
+ ) {
+ elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
+ for (k = 0; k < i; k++) rt_raster_destroy(rast[k]);
+ PG_RETURN_NULL();
+ }
+
+ /* SRID must match */
+ if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
+ elog(NOTICE, "The two rasters provided have different SRIDs");
+ for (k = 0; k < set_count; k++) rt_raster_destroy(rast[k]);
+ PG_RETURN_NULL();
+ }
+
+ /* raster extents need to intersect */
+ do {
+ initGEOS(lwnotice, lwgeom_geos_error);
+
+ rtn = 1;
+ for (i = 0; i < 2; i++) {
+ hull[i] = rt_raster_get_convex_hull(rast[i]);
+ if (NULL == hull[i]) {
+ for (j = 0; j < i; j++) {
+ GEOSGeom_destroy(ghull[j]);
+ lwpoly_free(hull[j]);
+ }
+ rtn = 0;
+ break;
+ }
+ ghull[i] = (GEOSGeometry *) LWGEOM2GEOS((LWGEOM *) hull[i]);
+ if (NULL == ghull[i]) {
+ for (j = 0; j < i; j++) {
+ GEOSGeom_destroy(ghull[j]);
+ lwpoly_free(hull[j]);
+ }
+ lwpoly_free(hull[i]);
+ rtn = 0;
+ break;
+ }
+ }
+ if (!rtn) break;
+
+ rtn = GEOSIntersects(ghull[0], ghull[1]);
+
+ for (i = 0; i < 2; i++) {
+ GEOSGeom_destroy(ghull[i]);
+ lwpoly_free(hull[i]);
+ }
+
+ if (rtn != 2) {
+ if (rtn != 1) {
+ for (k = 0; k < set_count; k++) rt_raster_destroy(rast[k]);
+ PG_RETURN_BOOL(0);
+ }
+ /* band isn't specified */
+ else if (!hasbandindex[0]) {
+ for (k = 0; k < set_count; k++) rt_raster_destroy(rast[k]);
+ PG_RETURN_BOOL(1);
+ }
+ }
+ }
+ while (0);
+
+ rtn = rt_raster_intersects(
+ rast[0], (hasbandindex[0] ? bandindex[0] - 1 : -1),
+ rast[1], (hasbandindex[1] ? bandindex[1] - 1 : -1),
+ &intersects
+ );
+ for (k = 0; k < set_count; k++) rt_raster_destroy(rast[k]);
+
+ if (!rtn) {
+ elog(ERROR, "RASTER_intersects: Unable to test for intersection on the two rasters");
+ PG_RETURN_NULL();
+ }
+
+ PG_RETURN_BOOL(intersects);
+}
+
/* ---------------------------------------------------------------- */
/* Memory allocation / error reporting hooks */
/* ---------------------------------------------------------------- */
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
gridx double precision DEFAULT NULL, gridy double precision DEFAULT NULL,
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
AS 'MODULE_PATHNAME', 'RASTER_asRaster'
pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
- AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $6, $7, $8, NULL, NULL, $4, $5, $9, $10) $$
+ AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $6, $7, $8, NULL, NULL, $4, $5, $9, $10, $11) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
- AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $4, $5, $6, $7, $8, NULL, NULL, $9, $10) $$
+ AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $4, $5, $6, $7, $8, NULL, NULL, $9, $10, $11) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
- AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $6, $7, $8, NULL, NULL, $4, $5, $9, $10) $$
+ AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $6, $7, $8, NULL, NULL, $4, $5, $9, $10, $11) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
- AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $4, $5, $6, $7, $8, NULL, NULL, $9, $10) $$
+ AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $4, $5, $6, $7, $8, NULL, NULL, $9, $10, $11) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
pixeltype text,
value double precision DEFAULT 1,
nodataval double precision DEFAULT 0,
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
- AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10) $$
+ AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10, $11) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
value double precision DEFAULT 1,
nodataval double precision DEFAULT 0,
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
- AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL, $9, $10) $$
+ AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL, $9, $10, $11) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
pixeltype text,
value double precision DEFAULT 1,
nodataval double precision DEFAULT 0,
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
- AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10) $$
+ AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10, $11) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
value double precision DEFAULT 1,
nodataval double precision DEFAULT 0,
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
- AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL,$9, $10) $$
+ AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL,$9, $10, $11) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
ref raster,
pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
value double precision[] DEFAULT ARRAY[1]::double precision[],
- nodataval double precision[] DEFAULT ARRAY[0]::double precision[]
+ nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
+ touched boolean DEFAULT FALSE
)
RETURNS raster
AS $$
g := geom;
END IF;
- RETURN _st_asraster(g, scale_x, scale_y, NULL, NULL, $3, $4, $5, NULL, NULL, ul_x, ul_y, skew_x, skew_y);
+ RETURN _st_asraster(g, scale_x, scale_y, NULL, NULL, $3, $4, $5, NULL, NULL, ul_x, ul_y, skew_x, skew_y, $6);
END;
$$ LANGUAGE 'plpgsql' STABLE;
ref raster,
pixeltype text,
value double precision DEFAULT 1,
- nodataval double precision DEFAULT 0
+ nodataval double precision DEFAULT 0,
+ touched boolean DEFAULT FALSE
)
RETURNS raster
- AS $$ SELECT st_asraster($1, $2, ARRAY[$3]::text[], ARRAY[$4]::double precision[], ARRAY[$5]::double precision[]) $$
+ AS $$ SELECT st_asraster($1, $2, ARRAY[$3]::text[], ARRAY[$4]::double precision[], ARRAY[$5]::double precision[], $6) $$
LANGUAGE 'sql' STABLE;
-----------------------------------------------------------------------
LANGUAGE 'plpgsql' IMMUTABLE STRICT;
---------------------------------------------------------------------------------
--- ST_World2RasteCoordX(rast raster, xw float8)
+-- ST_World2RasterCoordX(rast raster, xw float8)
-- Returns the column number of the pixels covering the provided world X coordinate
-- for a non-rotated raster.
-- This function works even if the world coordinate is outside the raster extent.
RESTRICT = contsel, JOIN = contjoinsel
);
+-----------------------------------------------------------------------
+-- Raster/Raster Spatial Relationship
+-----------------------------------------------------------------------
+CREATE OR REPLACE FUNCTION _st_intersects(rastA raster, nbandA integer, rastB raster, nbandB integer)
+ RETURNS boolean
+ AS 'MODULE_PATHNAME', 'RASTER_intersects'
+ LANGUAGE 'C' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_intersects(rastA raster, nbandA integer, rastB raster, nbandB integer)
+ RETURNS boolean
+ AS $$ SELECT $1 && $3 AND _st_intersects($1, $2, $3, $4) $$
+ LANGUAGE 'SQL' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_intersects(rastA raster, rastB raster, nbandA integer DEFAULT NULL, nbandB integer DEFAULT NULL)
+ RETURNS boolean
+ AS $$ SELECT $1 && $2 AND _st_intersects($1, $3, $2, $4) $$
+ LANGUAGE 'SQL' IMMUTABLE;
+
-----------------------------------------------------------------------
-- Raster/Geometry Spatial Relationship
-----------------------------------------------------------------------
+CREATE OR REPLACE FUNCTION _st_intersects(rast raster, geom geometry, nband integer DEFAULT NULL)
+ RETURNS boolean
+ AS $$
+ DECLARE
+ gr raster;
+ scale double precision;
+ BEGIN
+ IF ST_Intersects(geom, ST_ConvexHull(rast)) IS NOT TRUE THEN
+ RETURN FALSE;
+ ELSEIF nband IS NULL THEN
+ RETURN TRUE;
+ END IF;
+
+ -- scale is set to 1/100th of raster for granularity
+ SELECT CASE WHEN scaley < scalex THEN scaley ELSE scalex END / 100 INTO scale FROM ST_Metadata(rast);
+ gr := ST_AsRaster(geom, scale, scale);
+ IF gr IS NULL THEN
+ RAISE EXCEPTION 'Unable to convert geometry to a raster';
+ RETURN FALSE;
+ END IF;
+
+ RETURN ST_Intersects(rast, 1, gr, 1);
+ END;
+ $$ LANGUAGE 'plpgsql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_intersects(rast raster, geom geometry, nband integer DEFAULT NULL)
+ RETURNS boolean
+ AS $$ SELECT $1 && $2 AND _st_intersects($1, $2, $3) $$
+ LANGUAGE 'SQL' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_intersects(rast raster, nband integer, geom geometry)
+ RETURNS boolean
+ AS $$ SELECT $1 && $3 AND _st_intersects($1, $3, $2) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
-----------------------------------------------------------------------
--- _st_intersects(geomin geometry, rast raster, band integer, hasnodata boolean)
--- If hasnodata is true, check for the presence of withvalue pixels in the area
+-- _st_intersects(geom geometry, rast raster, nband integer)
+-- If band is provided, check for the presence of withvalue pixels in the area
-- shared by the raster and the geometry. If only nodata value pixels are found, the
-- geometry does not intersect with the raster.
-----------------------------------------------------------------------
-CREATE OR REPLACE FUNCTION _st_intersects(geomin geometry, rast raster, band integer, hasnodata boolean)
- RETURNS boolean AS
- $$
- DECLARE
- nodata float8 := 0.0;
- geomintersect geometry;
- x1w double precision := 0.0;
- x2w double precision := 0.0;
- y1w double precision := 0.0;
- y2w double precision := 0.0;
- x1 integer := 0;
- x2 integer := 0;
- x3 integer := 0;
- x4 integer := 0;
- y1 integer := 0;
- y2 integer := 0;
- y3 integer := 0;
- y4 integer := 0;
- x integer := 0;
- y integer := 0;
- xinc integer := 0;
- yinc integer := 0;
- pixelval double precision;
- bintersect boolean := FALSE;
- gtype text;
- scale float8;
- w int;
- h int;
- BEGIN
+-- This function can not be STRICT
+CREATE OR REPLACE FUNCTION _st_intersects(geom geometry, rast raster, nband integer DEFAULT NULL)
+ RETURNS boolean AS $$
+ DECLARE
+ hasnodata boolean := TRUE;
+ nodata float8 := 0.0;
+ convexhull geometry;
+ geomintersect geometry;
+ x1w double precision := 0.0;
+ x2w double precision := 0.0;
+ y1w double precision := 0.0;
+ y2w double precision := 0.0;
+ x1 integer := 0;
+ x2 integer := 0;
+ x3 integer := 0;
+ x4 integer := 0;
+ y1 integer := 0;
+ y2 integer := 0;
+ y3 integer := 0;
+ y4 integer := 0;
+ x integer := 0;
+ y integer := 0;
+ xinc integer := 0;
+ yinc integer := 0;
+ pixelval double precision;
+ bintersect boolean := FALSE;
+ gtype text;
+ scale float8;
+ w int;
+ h int;
+ BEGIN
+ convexhull := ST_ConvexHull(rast);
+ IF nband IS NOT NULL THEN
+ SELECT hasnodatavalue INTO hasnodata FROM ST_BandMetaData(rast, nband);
+ END IF;
- -- Get the intersection between with the geometry.
- -- We will search for withvalue pixel only in this area.
- geomintersect := st_intersection(geomin, st_convexhull(rast));
+ IF ST_Intersects(geom, convexhull) IS NOT TRUE THEN
+ RETURN FALSE;
+ ELSEIF nband IS NULL OR hasnodata IS FALSE THEN
+ RETURN TRUE;
+ END IF;
---RAISE NOTICE 'geomintersect1=%', astext(geomintersect);
+ -- Get the intersection between with the geometry.
+ -- We will search for withvalue pixel only in this area.
+ geomintersect := st_intersection(geom, convexhull);
- -- If the intersection is empty, return false
- IF st_isempty(geomintersect) THEN
- RETURN FALSE;
- END IF;
+--RAISE NOTICE 'geomintersect=%', st_astext(geomintersect);
- -- If the band does not have a nodatavalue, there is no need to search for with value pixels
- IF NOT hasnodata OR st_bandnodatavalue(rast, band) IS NULL THEN
- RETURN TRUE;
- END IF;
+ -- If the intersection is empty, return false
+ IF st_isempty(geomintersect) THEN
+ RETURN FALSE;
+ END IF;
- -- We create a minimalistic buffer around the intersection in order to scan every pixels
- -- that would touch the edge or intersect with the geometry
- SELECT (scalex * skewy), width, height INTO scale, w, h FROM ST_Metadata(rast);
- geomintersect := st_buffer(geomintersect, scale / 1000000);
+ -- We create a minimalistic buffer around the intersection in order to scan every pixels
+ -- that would touch the edge or intersect with the geometry
+ SELECT (scalex * skewy), width, height INTO scale, w, h FROM ST_Metadata(rast);
+ geomintersect := st_buffer(geomintersect, scale / 1000000);
--RAISE NOTICE 'geomintersect2=%', astext(geomintersect);
- -- Find the world coordinates of the bounding box of the intersecting area
- x1w := st_xmin(geomintersect);
- y1w := st_ymin(geomintersect);
- x2w := st_xmax(geomintersect);
- y2w := st_ymax(geomintersect);
- nodata := st_bandnodatavalue(rast, band);
+ -- Find the world coordinates of the bounding box of the intersecting area
+ x1w := st_xmin(geomintersect);
+ y1w := st_ymin(geomintersect);
+ x2w := st_xmax(geomintersect);
+ y2w := st_ymax(geomintersect);
+ nodata := st_bandnodatavalue(rast, nband);
--RAISE NOTICE 'x1w=%, y1w=%, x2w=%, y2w=%', x1w, y1w, x2w, y2w;
- -- Convert world coordinates to raster coordinates
- x1 := st_world2rastercoordx(rast, x1w, y1w);
- y1 := st_world2rastercoordy(rast, x1w, y1w);
- x2 := st_world2rastercoordx(rast, x2w, y1w);
- y2 := st_world2rastercoordy(rast, x2w, y1w);
- x3 := st_world2rastercoordx(rast, x1w, y2w);
- y3 := st_world2rastercoordy(rast, x1w, y2w);
- x4 := st_world2rastercoordx(rast, x2w, y2w);
- y4 := st_world2rastercoordy(rast, x2w, y2w);
+ -- Convert world coordinates to raster coordinates
+ x1 := st_world2rastercoordx(rast, x1w, y1w);
+ y1 := st_world2rastercoordy(rast, x1w, y1w);
+ x2 := st_world2rastercoordx(rast, x2w, y1w);
+ y2 := st_world2rastercoordy(rast, x2w, y1w);
+ x3 := st_world2rastercoordx(rast, x1w, y2w);
+ y3 := st_world2rastercoordy(rast, x1w, y2w);
+ x4 := st_world2rastercoordx(rast, x2w, y2w);
+ y4 := st_world2rastercoordy(rast, x2w, y2w);
--RAISE NOTICE 'x1=%, y1=%, x2=%, y2=%, x3=%, y3=%, x4=%, y4=%', x1, y1, x2, y2, x3, y3, x4, y4;
- -- Order the raster coordinates for the upcoming FOR loop.
- x1 := int4smaller(int4smaller(int4smaller(x1, x2), x3), x4);
- y1 := int4smaller(int4smaller(int4smaller(y1, y2), y3), y4);
- x2 := int4larger(int4larger(int4larger(x1, x2), x3), x4);
- y2 := int4larger(int4larger(int4larger(y1, y2), y3), y4);
+ -- Order the raster coordinates for the upcoming FOR loop.
+ x1 := int4smaller(int4smaller(int4smaller(x1, x2), x3), x4);
+ y1 := int4smaller(int4smaller(int4smaller(y1, y2), y3), y4);
+ x2 := int4larger(int4larger(int4larger(x1, x2), x3), x4);
+ y2 := int4larger(int4larger(int4larger(y1, y2), y3), y4);
- -- Make sure the range is not lower than 1.
- -- This can happen when world coordinate are exactly on the left border
- -- of the raster and that they do not span on more than one pixel.
- x1 := int4smaller(int4larger(x1, 1), w);
- y1 := int4smaller(int4larger(y1, 1), h);
+ -- Make sure the range is not lower than 1.
+ -- This can happen when world coordinate are exactly on the left border
+ -- of the raster and that they do not span on more than one pixel.
+ x1 := int4smaller(int4larger(x1, 1), w);
+ y1 := int4smaller(int4larger(y1, 1), h);
- -- Also make sure the range does not exceed the width and height of the raster.
- -- This can happen when world coordinate are exactly on the lower right border
- -- of the raster.
- x2 := int4smaller(x2, w);
- y2 := int4smaller(y2, h);
+ -- Also make sure the range does not exceed the width and height of the raster.
+ -- This can happen when world coordinate are exactly on the lower right border
+ -- of the raster.
+ x2 := int4smaller(x2, w);
+ y2 := int4smaller(y2, h);
--RAISE NOTICE 'x1=%, y1=%, x2=%, y2=%', x1, y1, x2, y2;
- -- Search exhaustively for withvalue pixel on a moving 3x3 grid
- -- (very often more efficient than searching on a mere 1x1 grid)
- FOR xinc in 0..2 LOOP
- FOR yinc in 0..2 LOOP
- FOR x IN x1+xinc..x2 BY 3 LOOP
- FOR y IN y1+yinc..y2 BY 3 LOOP
- -- Check first if the pixel intersects with the geometry. Often many won't.
- bintersect := NOT st_isempty(st_intersection(st_pixelaspolygon(rast, band, x, y), geomin));
-
- IF bintersect THEN
- -- If the pixel really intersects, check its value. Return TRUE if with value.
- pixelval := st_value(rast, band, x, y);
- IF pixelval != nodata THEN
- RETURN TRUE;
- END IF;
- END IF;
- END LOOP;
- END LOOP;
- END LOOP;
- END LOOP;
- RETURN FALSE;
- END;
- $$
- LANGUAGE 'plpgsql' IMMUTABLE STRICT;
-
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(geometry, raster, integer)
- RETURNS boolean AS
- $$ SELECT $1 && $2 AND _st_intersects($1, $2, $3, TRUE);
- $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(raster, integer, geometry)
- RETURNS boolean AS
- $$ SELECT st_intersects($3, $1, $2);
- $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(geometry, raster)
- RETURNS boolean AS
- $$ SELECT st_intersects($1, $2, 1);
- $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(raster, geometry)
- RETURNS boolean AS
- $$ SELECT st_intersects($2, $1, 1);
- $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(geometry, raster, integer, boolean)
- RETURNS boolean AS
- $$ SELECT $1 && $2 AND _st_intersects($1, $2, $3, $4);
- $$ LANGUAGE 'SQL' IMMUTABLE;
-
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(raster, integer, boolean, geometry)
- RETURNS boolean AS
- $$ SELECT st_intersects($4, $1, $2, $3);
- $$ LANGUAGE 'SQL' IMMUTABLE;
+ -- Search exhaustively for withvalue pixel on a moving 3x3 grid
+ -- (very often more efficient than searching on a mere 1x1 grid)
+ FOR xinc in 0..2 LOOP
+ FOR yinc in 0..2 LOOP
+ FOR x IN x1+xinc..x2 BY 3 LOOP
+ FOR y IN y1+yinc..y2 BY 3 LOOP
+ -- Check first if the pixel intersects with the geometry. Often many won't.
+ bintersect := NOT st_isempty(st_intersection(st_pixelaspolygon(rast, nband, x, y), geom));
+
+ IF bintersect THEN
+ -- If the pixel really intersects, check its value. Return TRUE if with value.
+ pixelval := st_value(rast, nband, x, y);
+ IF pixelval != nodata THEN
+ RETURN TRUE;
+ END IF;
+ END IF;
+ END LOOP;
+ END LOOP;
+ END LOOP;
+ END LOOP;
--- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(geometry, raster, boolean)
- RETURNS boolean AS
- $$ SELECT st_intersects($1, $2, 1, $3);
- $$ LANGUAGE 'SQL' IMMUTABLE;
+ RETURN FALSE;
+ END;
+ $$ LANGUAGE 'plpgsql' IMMUTABLE;
-- This function can not be STRICT
-CREATE OR REPLACE FUNCTION st_intersects(raster, boolean, geometry)
- RETURNS boolean AS
- $$ SELECT st_intersects($3, $1, 1, $2);
- $$ LANGUAGE 'SQL' IMMUTABLE;
+CREATE OR REPLACE FUNCTION st_intersects(geom geometry, rast raster, nband integer DEFAULT NULL)
+ RETURNS boolean AS
+ $$ SELECT $1 && $2 AND _st_intersects($1, $2, $3); $$
+ LANGUAGE 'SQL' IMMUTABLE;
-----------------------------------------------------------------------
-- _st_intersection(geom geometry, rast raster, band integer)
-- Raster nodata value areas are not vectorized and hence do not intersect
-- with any geometries.
-----------------------------------------------------------------------
-CREATE OR REPLACE FUNCTION ST_Intersection(geomin geometry, rast raster, band integer)
- RETURNS SETOF geomval AS
-$BODY$
- DECLARE
- intersects boolean := FALSE;
- BEGIN
- intersects := ST_Intersects(geomin, rast, band);
- IF intersects THEN
- -- Return the intersections of the geometry with the vectorized parts of
- -- the raster and the values associated with those parts, if really their
- -- intersection is not empty.
- RETURN QUERY
- SELECT intgeom,
- val
- FROM (SELECT ST_Intersection((gv).geom, geomin) AS intgeom,
- (gv).val
- FROM ST_DumpAsPolygons(rast, band) gv
- WHERE ST_Intersects((gv).geom, geomin)
- ) foo
- WHERE NOT ST_IsEmpty(intgeom);
- ELSE
- -- If the geometry does not intersect with the raster, return an empty
- -- geometry and a null value
- RETURN QUERY
- SELECT emptygeom,
- NULL::float8
- FROM ST_GeomCollFromText('GEOMETRYCOLLECTION EMPTY', ST_SRID($1)) emptygeom;
- END IF;
- END;
- $BODY$
- LANGUAGE 'plpgsql' IMMUTABLE STRICT;
+CREATE OR REPLACE FUNCTION st_intersection(geomin geometry, rast raster, band integer DEFAULT 1)
+ RETURNS SETOF geomval AS $$
+ DECLARE
+ intersects boolean := FALSE;
+ BEGIN
+ intersects := ST_Intersects(geomin, rast, band);
+ IF intersects THEN
+ -- Return the intersections of the geometry with the vectorized parts of
+ -- the raster and the values associated with those parts, if really their
+ -- intersection is not empty.
+ RETURN QUERY
+ SELECT
+ intgeom,
+ val
+ FROM (
+ SELECT
+ ST_Intersection((gv).geom, geomin) AS intgeom,
+ (gv).val
+ FROM ST_DumpAsPolygons(rast, band) gv
+ WHERE ST_Intersects((gv).geom, geomin)
+ ) foo
+ WHERE NOT ST_IsEmpty(intgeom);
+ ELSE
+ -- If the geometry does not intersect with the raster, return an empty
+ -- geometry and a null value
+ RETURN QUERY
+ SELECT
+ emptygeom,
+ NULL::float8
+ FROM ST_GeomCollFromText('GEOMETRYCOLLECTION EMPTY', ST_SRID($1)) emptygeom;
+ END IF;
+ END;
+ $$
+ LANGUAGE 'plpgsql' IMMUTABLE STRICT;
CREATE OR REPLACE FUNCTION st_intersection(rast raster, geom geometry)
- RETURNS SETOF geomval AS
- $$
- SELECT (gv).geom, (gv).val FROM st_intersection($2, $1, 1) gv;
- $$
- LANGUAGE SQL IMMUTABLE STRICT;
-
-CREATE OR REPLACE FUNCTION st_intersection(geom geometry, rast raster)
- RETURNS SETOF geomval AS
- $$
- SELECT (gv).geom, (gv).val FROM st_intersection($1, $2, 1) gv;
- $$
- LANGUAGE SQL IMMUTABLE STRICT;
+ RETURNS SETOF geomval
+ AS $$ SELECT (gv).geom, (gv).val FROM st_intersection($2, $1, 1) gv; $$
+ LANGUAGE SQL IMMUTABLE STRICT;
CREATE OR REPLACE FUNCTION st_intersection(rast raster, band integer, geom geometry)
- RETURNS SETOF geomval AS
- $$
- SELECT (gv).geom, (gv).val FROM st_intersection($3, $1, $2) gv;
- $$
- LANGUAGE SQL IMMUTABLE STRICT;
+ RETURNS SETOF geomval
+ AS $$ SELECT (gv).geom, (gv).val FROM st_intersection($3, $1, $2) gv; $$
+ LANGUAGE SQL IMMUTABLE STRICT;
------------------------------------------------------------------------------
-- RASTER_COLUMNS
qlls = NULL;
qlls_count = 0;
+ deepRelease(raster);
+
xmax = 100;
ymax = 100;
raster = rt_raster_new(xmax, ymax);
&scale_x, &scale_y,
NULL, NULL,
NULL, NULL,
- NULL, NULL
+ NULL, NULL,
+ NULL
);
free(wkb);
deepRelease(raster);
}
+static void testIntersects() {
+ rt_raster rast1;
+ rt_raster rast2;
+ rt_band band1;
+ rt_band band2;
+ double nodata;
+ int rtn;
+ int intersects;
+
+ /*
+ rast1
+
+ (-1, -1)
+ +-+-+
+ |1|1|
+ +-+-+
+ |1|1|
+ +-+-+
+ (1, 1)
+ */
+ rast1 = rt_raster_new(2, 2);
+ assert(rast1);
+ rt_raster_set_offsets(rast1, -1, -1);
+
+ band1 = addBand(rast1, PT_8BUI, 1, 0);
+ CHECK(band1);
+ rt_band_set_nodata(band1, 0);
+ rtn = rt_band_set_pixel(band1, 0, 0, 1);
+ rtn = rt_band_set_pixel(band1, 0, 1, 1);
+ rtn = rt_band_set_pixel(band1, 1, 0, 1);
+ rtn = rt_band_set_pixel(band1, 1, 1, 1);
+
+ nodata = rt_band_get_nodata(band1);
+ CHECK_EQUALS(nodata, 0);
+
+ /*
+ rast2
+
+ (0, 0)
+ +-+-+
+ |1|1|
+ +-+-+
+ |1|1|
+ +-+-+
+ (2, 2)
+ */
+ rast2 = rt_raster_new(2, 2);
+ assert(rast2);
+
+ band2 = addBand(rast2, PT_8BUI, 1, 0);
+ CHECK(band2);
+ rt_band_set_nodata(band2, 0);
+ rtn = rt_band_set_pixel(band2, 0, 0, 1);
+ rtn = rt_band_set_pixel(band2, 0, 1, 1);
+ rtn = rt_band_set_pixel(band2, 1, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 1, 1);
+
+ nodata = rt_band_get_nodata(band2);
+ CHECK_EQUALS(nodata, 0);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ rtn = rt_raster_intersects(
+ rast1, -1,
+ rast2, -1,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /*
+ rast2
+
+ (0, 0)
+ +-+-+
+ |0|1|
+ +-+-+
+ |1|1|
+ +-+-+
+ (2, 2)
+ */
+ rtn = rt_band_set_pixel(band2, 0, 0, 0);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /*
+ rast2
+
+ (0, 0)
+ +-+-+
+ |1|0|
+ +-+-+
+ |1|1|
+ +-+-+
+ (2, 2)
+ */
+ rtn = rt_band_set_pixel(band2, 0, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 0, 0);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /*
+ rast2
+
+ (0, 0)
+ +-+-+
+ |0|0|
+ +-+-+
+ |0|1|
+ +-+-+
+ (2, 2)
+ */
+ rtn = rt_band_set_pixel(band2, 0, 0, 0);
+ rtn = rt_band_set_pixel(band2, 1, 0, 0);
+ rtn = rt_band_set_pixel(band2, 0, 1, 0);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /*
+ rast2
+
+ (0, 0)
+ +-+-+
+ |0|0|
+ +-+-+
+ |0|0|
+ +-+-+
+ (2, 2)
+ */
+ rtn = rt_band_set_pixel(band2, 0, 0, 0);
+ rtn = rt_band_set_pixel(band2, 1, 0, 0);
+ rtn = rt_band_set_pixel(band2, 0, 1, 0);
+ rtn = rt_band_set_pixel(band2, 1, 1, 0);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects != 1));
+
+ /*
+ rast2
+
+ (2, 0)
+ +-+-+
+ |1|1|
+ +-+-+
+ |1|1|
+ +-+-+
+ (4, 2)
+ */
+ rt_raster_set_offsets(rast2, 2, 0);
+
+ rtn = rt_band_set_pixel(band2, 0, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 0, 1);
+ rtn = rt_band_set_pixel(band2, 0, 1, 1);
+ rtn = rt_band_set_pixel(band2, 1, 1, 1);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects != 1));
+
+ /*
+ rast2
+
+ (0.1, 0.1)
+ +-+-+
+ |1|1|
+ +-+-+
+ |1|1|
+ +-+-+
+ (0.9, 0.9)
+ */
+ rt_raster_set_offsets(rast2, 0.1, 0.1);
+ rt_raster_set_scale(rast2, 0.4, 0.4);
+
+ rtn = rt_band_set_pixel(band2, 0, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 0, 1);
+ rtn = rt_band_set_pixel(band2, 0, 1, 1);
+ rtn = rt_band_set_pixel(band2, 1, 1, 1);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /*
+ rast2
+
+ (-0.1, 0.1)
+ +-+-+
+ |1|1|
+ +-+-+
+ |1|1|
+ +-+-+
+ (0.9, 0.9)
+ */
+ rt_raster_set_offsets(rast2, -0.1, 0.1);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ deepRelease(rast2);
+
+ /*
+ rast2
+
+ (0, 0)
+ +-+-+-+
+ |1|1|1|
+ +-+-+-+
+ |1|1|1|
+ +-+-+-+
+ |1|1|1|
+ +-+-+-+
+ (3, 3)
+ */
+ rast2 = rt_raster_new(3, 3);
+ assert(rast2);
+
+ band2 = addBand(rast2, PT_8BUI, 1, 0);
+ CHECK(band2);
+ rt_band_set_nodata(band2, 0);
+ rtn = rt_band_set_pixel(band2, 0, 0, 1);
+ rtn = rt_band_set_pixel(band2, 0, 1, 1);
+ rtn = rt_band_set_pixel(band2, 0, 2, 1);
+ rtn = rt_band_set_pixel(band2, 1, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 1, 1);
+ rtn = rt_band_set_pixel(band2, 1, 2, 1);
+ rtn = rt_band_set_pixel(band2, 2, 0, 1);
+ rtn = rt_band_set_pixel(band2, 2, 1, 1);
+ rtn = rt_band_set_pixel(band2, 2, 2, 1);
+
+ nodata = rt_band_get_nodata(band2);
+ CHECK_EQUALS(nodata, 0);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /*
+ rast2
+
+ (-2, -2)
+ +-+-+-+
+ |1|1|1|
+ +-+-+-+
+ |1|1|1|
+ +-+-+-+
+ |1|1|1|
+ +-+-+-+
+ (1, 1)
+ */
+ rt_raster_set_offsets(rast2, -2, -2);
+
+ rtn = rt_band_set_pixel(band2, 0, 0, 1);
+ rtn = rt_band_set_pixel(band2, 0, 1, 1);
+ rtn = rt_band_set_pixel(band2, 0, 2, 1);
+ rtn = rt_band_set_pixel(band2, 1, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 1, 1);
+ rtn = rt_band_set_pixel(band2, 1, 2, 1);
+ rtn = rt_band_set_pixel(band2, 2, 0, 1);
+ rtn = rt_band_set_pixel(band2, 2, 1, 1);
+ rtn = rt_band_set_pixel(band2, 2, 2, 1);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /*
+ rast2
+
+ (-2, -2)
+ +-+-+-+
+ |0|1|1|
+ +-+-+-+
+ |1|0|1|
+ +-+-+-+
+ |1|1|0|
+ +-+-+-+
+ (1, 1)
+ */
+ rtn = rt_band_set_pixel(band2, 0, 0, 0);
+ rtn = rt_band_set_pixel(band2, 0, 1, 1);
+ rtn = rt_band_set_pixel(band2, 0, 2, 1);
+ rtn = rt_band_set_pixel(band2, 1, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 1, 0);
+ rtn = rt_band_set_pixel(band2, 1, 2, 1);
+ rtn = rt_band_set_pixel(band2, 2, 0, 1);
+ rtn = rt_band_set_pixel(band2, 2, 1, 1);
+ rtn = rt_band_set_pixel(band2, 2, 2, 0);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /*
+ rast2
+
+ (-2, -2)
+ +-+-+-+
+ |0|1|1|
+ +-+-+-+
+ |1|0|0|
+ +-+-+-+
+ |1|0|0|
+ +-+-+-+
+ (1, 1)
+ */
+ rtn = rt_band_set_pixel(band2, 0, 0, 0);
+ rtn = rt_band_set_pixel(band2, 0, 1, 1);
+ rtn = rt_band_set_pixel(band2, 0, 2, 1);
+ rtn = rt_band_set_pixel(band2, 1, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 1, 0);
+ rtn = rt_band_set_pixel(band2, 1, 2, 0);
+ rtn = rt_band_set_pixel(band2, 2, 0, 1);
+ rtn = rt_band_set_pixel(band2, 2, 1, 0);
+ rtn = rt_band_set_pixel(band2, 2, 2, 0);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /*
+ rast2
+
+ (-2, -2)
+ +-+-+-+
+ |0|1|0|
+ +-+-+-+
+ |1|0|0|
+ +-+-+-+
+ |0|0|0|
+ +-+-+-+
+ (1, 1)
+ */
+ rtn = rt_band_set_pixel(band2, 0, 0, 0);
+ rtn = rt_band_set_pixel(band2, 0, 1, 1);
+ rtn = rt_band_set_pixel(band2, 0, 2, 0);
+ rtn = rt_band_set_pixel(band2, 1, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 1, 0);
+ rtn = rt_band_set_pixel(band2, 1, 2, 0);
+ rtn = rt_band_set_pixel(band2, 2, 0, 0);
+ rtn = rt_band_set_pixel(band2, 2, 1, 0);
+ rtn = rt_band_set_pixel(band2, 2, 2, 0);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ deepRelease(rast2);
+
+ /* skew tests */
+ /* rast2 (skewed by -0.5, 0.5) */
+ rast2 = rt_raster_new(3, 3);
+ assert(rast2);
+ rt_raster_set_skews(rast2, -0.5, 0.5);
+
+ band2 = addBand(rast2, PT_8BUI, 1, 0);
+ CHECK(band2);
+ rt_band_set_nodata(band2, 0);
+ rtn = rt_band_set_pixel(band2, 0, 0, 1);
+ rtn = rt_band_set_pixel(band2, 0, 1, 2);
+ rtn = rt_band_set_pixel(band2, 0, 2, 3);
+ rtn = rt_band_set_pixel(band2, 1, 0, 1);
+ rtn = rt_band_set_pixel(band2, 1, 1, 2);
+ rtn = rt_band_set_pixel(band2, 1, 2, 3);
+ rtn = rt_band_set_pixel(band2, 2, 0, 1);
+ rtn = rt_band_set_pixel(band2, 2, 1, 2);
+ rtn = rt_band_set_pixel(band2, 2, 2, 3);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /* rast2 (skewed by -1, 1) */
+ rt_raster_set_skews(rast2, -1, 1);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ /* rast2 (skewed by 1, -1) */
+ rt_raster_set_skews(rast2, 1, -1);
+
+ rtn = rt_raster_intersects(
+ rast1, 0,
+ rast2, 0,
+ &intersects
+ );
+ CHECK((rtn != 0));
+ CHECK((intersects == 1));
+
+ deepRelease(rast2);
+ deepRelease(rast1);
+}
+
int
main()
{
testGDALRasterize();
printf("Successfully tested rt_raster_gdal_rasterize\n");
+ printf("Testing rt_raster_intersects\n");
+ testIntersects();
+ printf("Successfully tested rt_raster_intersects\n");
+
deepRelease(raster);
return EXIT_SUCCESS;
TEST_SREL = \
rt_spatial_relationship.sql \
+ rt_intersects.sql \
$(NULL)
TEST_BUGS = \
--- /dev/null
+DROP TABLE IF EXISTS raster_intersects_rast;
+DROP TABLE IF EXISTS raster_intersects_geom;
+CREATE TABLE raster_intersects_rast (
+ rid integer,
+ rast raster
+);
+CREATE TABLE raster_intersects_geom (
+ gid integer,
+ geom geometry
+);
+CREATE OR REPLACE FUNCTION make_test_raster(rid integer, width integer DEFAULT 2, height integer DEFAULT 2, ul_x double precision DEFAULT 0, ul_y double precision DEFAULT 0, skew_x double precision DEFAULT 0, skew_y double precision DEFAULT 0)
+ RETURNS void
+ AS $$
+ DECLARE
+ x int;
+ y int;
+ rast raster;
+ BEGIN
+ rast := ST_MakeEmptyRaster(width, height, ul_x, ul_y, 1, 1, skew_x, skew_y, -1);
+ rast := ST_AddBand(rast, 1, '8BUI', 1, 0);
+
+
+ INSERT INTO raster_intersects_rast VALUES (rid, rast);
+
+ RETURN;
+ END;
+ $$ LANGUAGE 'plpgsql';
+SELECT make_test_raster(0, 2, 2, -1, -1);
+SELECT make_test_raster(1, 2, 2);
+SELECT make_test_raster(2, 3, 3);
+DROP FUNCTION make_test_raster(integer, integer, integer, double precision, double precision, double precision, double precision);
+
+INSERT INTO raster_intersects_rast VALUES (10, (
+ SELECT
+ ST_SetValue(rast, 1, 1, 1, 0)
+ FROM raster_intersects_rast
+ WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (11, (
+ SELECT
+ ST_SetValue(rast, 1, 2, 1, 0)
+ FROM raster_intersects_rast
+ WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (12, (
+ SELECT
+ ST_SetValue(
+ ST_SetValue(
+ ST_SetValue(rast, 1, 1, 1, 0),
+ 1, 2, 1, 0
+ ),
+ 1, 1, 2, 0
+ )
+ FROM raster_intersects_rast
+ WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (13, (
+ SELECT
+ ST_SetValue(
+ ST_SetValue(
+ ST_SetValue(
+ ST_SetValue(rast, 1, 1, 1, 0),
+ 1, 2, 1, 0
+ ),
+ 1, 1, 2, 0
+ ),
+ 1, 2, 2, 0
+ )
+ FROM raster_intersects_rast
+ WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (14, (
+ SELECT
+ ST_SetUpperLeft(rast, 2, 0)
+ FROM raster_intersects_rast
+ WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (15, (
+ SELECT
+ ST_SetScale(
+ ST_SetUpperLeft(rast, 0.1, 0.1),
+ 0.4, 0.4
+ )
+ FROM raster_intersects_rast
+ WHERE rid = 1
+));
+INSERT INTO raster_intersects_rast VALUES (16, (
+ SELECT
+ ST_SetScale(
+ ST_SetUpperLeft(rast, -0.1, 0.1),
+ 0.4, 0.4
+ )
+ FROM raster_intersects_rast
+ WHERE rid = 1
+));
+
+INSERT INTO raster_intersects_rast VALUES (20, (
+ SELECT
+ ST_SetUpperLeft(rast, -2, -2)
+ FROM raster_intersects_rast
+ WHERE rid = 2
+));
+INSERT INTO raster_intersects_rast VALUES (21, (
+ SELECT
+ ST_SetValue(
+ ST_SetValue(
+ ST_SetValue(rast, 1, 1, 1, 0),
+ 1, 2, 2, 0
+ ),
+ 1, 3, 3, 0
+ )
+ FROM raster_intersects_rast
+ WHERE rid = 20
+));
+INSERT INTO raster_intersects_rast VALUES (22, (
+ SELECT
+ ST_SetValue(
+ ST_SetValue(
+ rast, 1, 3, 2, 0
+ ),
+ 1, 2, 3, 0
+ )
+ FROM raster_intersects_rast
+ WHERE rid = 21
+));
+INSERT INTO raster_intersects_rast VALUES (23, (
+ SELECT
+ ST_SetValue(
+ ST_SetValue(
+ rast, 1, 3, 1, 0
+ ),
+ 1, 1, 3, 0
+ )
+ FROM raster_intersects_rast
+ WHERE rid = 22
+));
+
+INSERT INTO raster_intersects_rast VALUES (30, (
+ SELECT
+ ST_SetSkew(rast, -0.5, 0.5)
+ FROM raster_intersects_rast
+ WHERE rid = 2
+));
+INSERT INTO raster_intersects_rast VALUES (31, (
+ SELECT
+ ST_SetSkew(rast, -1, 1)
+ FROM raster_intersects_rast
+ WHERE rid = 2
+));
+INSERT INTO raster_intersects_rast VALUES (32, (
+ SELECT
+ ST_SetSkew(rast, 1, -1)
+ FROM raster_intersects_rast
+ WHERE rid = 2
+));
+
+SELECT
+ '1.1',
+ r1.rid,
+ r2.rid,
+ ST_Intersects(r1.rast, r2.rast)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_rast r2
+ ON r1.rid != r2.rid
+WHERE r1.rid = 0;
+
+SELECT
+ '1.2',
+ r1.rid,
+ r2.rid,
+ ST_Intersects(r1.rast, 1, r2.rast, 1)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_rast r2
+ ON r1.rid != r2.rid
+WHERE r1.rid = 0;
+
+-- point
+INSERT INTO raster_intersects_geom VALUES (
+ 1, (
+ SELECT ST_MakePoint(0, 0)
+ )
+), (
+ 2, (
+ SELECT ST_MakePoint(0.1, 0.1)
+ )
+), (
+ 3, (
+ SELECT ST_MakePoint(-0.1, -0.1)
+ )
+), (
+ 4, (
+ SELECT ST_MakePoint(-1, -1)
+ )
+), (
+ 5, (
+ SELECT ST_MakePoint(-1.1, -1)
+ )
+), (
+ 6, (
+ SELECT ST_MakePoint(-1, -1.1)
+ )
+), (
+ 7, (
+ SELECT ST_MakePoint(-1.5, -1.5)
+ )
+), (
+ 8, (
+ SELECT ST_MakePoint(3, 3)
+ )
+);
+
+-- multipoint
+INSERT INTO raster_intersects_geom VALUES (
+ 11, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 1 AND 10
+ )
+), (
+ 12, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 3 AND 10
+ )
+), (
+ 13, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 4 AND 10
+ )
+), (
+ 14, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 5 AND 10
+ )
+), (
+ 15, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 6 AND 10
+ )
+);
+
+-- linestring
+INSERT INTO raster_intersects_geom VALUES (
+ 21, (
+ SELECT ST_MakeLine(ARRAY[
+ ST_MakePoint(1, 1),
+ ST_MakePoint(1, 0)
+ ])
+ )
+), (
+ 22, (
+ SELECT ST_MakeLine(ARRAY[
+ ST_MakePoint(-1, -1),
+ ST_MakePoint(1, 1),
+ ST_MakePoint(1, 0)
+ ])
+ )
+), (
+ 23, (
+ SELECT ST_MakeLine(ARRAY[
+ ST_MakePoint(-1, -1),
+ ST_MakePoint(-1, 1),
+ ST_MakePoint(1, 1),
+ ST_MakePoint(1, -1)
+ ])
+ )
+), (
+ 24, (
+ SELECT ST_MakeLine(ARRAY[
+ ST_MakePoint(-1.1, 1.1),
+ ST_MakePoint(1.1, 1.1),
+ ST_MakePoint(1.1, -1.1),
+ ST_MakePoint(-1.1, -1.1),
+ ST_MakePoint(-1.1, 1.1)
+ ])
+ )
+), (
+ 25, (
+ SELECT ST_MakeLine(ARRAY[
+ ST_MakePoint(-2, 1),
+ ST_MakePoint(1, 2),
+ ST_MakePoint(2, -1),
+ ST_MakePoint(-1, -2),
+ ST_MakePoint(-2, 1)
+ ])
+ )
+), (
+ 26, (
+ SELECT ST_MakeLine(ARRAY[
+ ST_MakePoint(-0.5, 0.5),
+ ST_MakePoint(0, 0.5),
+ ST_MakePoint(0, 0),
+ ST_MakePoint(0, -0.5),
+ ST_MakePoint(-0.5, 0.5)
+ ])
+ )
+), (
+ 27, (
+ SELECT ST_MakeLine(ARRAY[
+ ST_MakePoint(0.5, 0.5),
+ ST_MakePoint(1, 1),
+ ST_MakePoint(1, 0),
+ ST_MakePoint(0.5, 0.5)
+ ])
+ )
+), (
+ 28, (
+ SELECT ST_MakeLine(ARRAY[
+ ST_MakePoint(1, 1),
+ ST_MakePoint(0, 2),
+ ST_MakePoint(1, 2),
+ ST_MakePoint(1, 1)
+ ])
+ )
+), (
+ 29, (
+ SELECT ST_MakeLine(ARRAY[
+ ST_MakePoint(0, 2),
+ ST_MakePoint(1, 2),
+ ST_MakePoint(1, 4),
+ ST_MakePoint(0, 2)
+ ])
+ )
+);
+
+-- polygon
+INSERT INTO raster_intersects_geom VALUES (
+ 31, (
+ SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 24
+ )
+), (
+ 32, (
+ SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 25
+ )
+), (
+ 33, (
+ SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 26
+ )
+), (
+ 34, (
+ SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 27
+ )
+), (
+ 35, (
+ SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 28
+ )
+), (
+ 36, (
+ SELECT ST_MakePolygon(geom) FROM raster_intersects_geom WHERE gid = 29
+ )
+);
+
+-- multipolygon
+INSERT INTO raster_intersects_geom VALUES (
+ 41, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 31 and 40
+ )
+), (
+ 42, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 32 and 40
+ )
+), (
+ 43, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 33 and 40
+ )
+), (
+ 44, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 34 and 40
+ )
+), (
+ 45, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 35 and 40
+ )
+), (
+ 46, (
+ SELECT ST_Collect(geom) FROM raster_intersects_geom WHERE gid BETWEEN 36 and 40
+ )
+);
+
+SELECT
+ '2.1',
+ r1.rid,
+ g1.gid,
+ ST_GeometryType(g1.geom),
+ ST_Intersects(r1.rast, g1.geom)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+ ON 1 = 1
+WHERE r1.rid = 0;
+
+SELECT
+ '2.2',
+ r1.rid,
+ g1.gid,
+ ST_GeometryType(g1.geom),
+ ST_Intersects(g1.geom, r1.rast)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+ ON 1 = 1
+WHERE r1.rid = 0;
+
+SELECT
+ '2.3',
+ r1.rid,
+ g1.gid,
+ ST_GeometryType(g1.geom),
+ ST_Intersects(r1.rast, g1.geom)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+ ON 1 = 1
+WHERE r1.rid = 2;
+
+SELECT
+ '2.4',
+ r1.rid,
+ g1.gid,
+ ST_GeometryType(g1.geom),
+ ST_Intersects(g1.geom, r1.rast)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+ ON 1 = 1
+WHERE r1.rid = 2;
+
+SELECT
+ '2.5',
+ r1.rid,
+ g1.gid,
+ ST_GeometryType(g1.geom),
+ ST_Intersects(r1.rast, g1.geom, 1)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+ ON 1 = 1
+WHERE r1.rid = 0;
+
+SELECT
+ '2.6',
+ r1.rid,
+ g1.gid,
+ ST_GeometryType(g1.geom),
+ ST_Intersects(r1.rast, g1.geom, 1)
+FROM raster_intersects_rast r1
+JOIN raster_intersects_geom g1
+ ON 1 = 1
+WHERE r1.rid = 2;
--- /dev/null
+NOTICE: table "raster_intersects_rast" does not exist, skipping
+NOTICE: table "raster_intersects_geom" does not exist, skipping
+1.1|0|1|t
+1.1|0|2|t
+1.1|0|10|t
+1.1|0|11|t
+1.1|0|12|t
+1.1|0|13|t
+1.1|0|14|f
+1.1|0|15|t
+1.1|0|16|t
+1.1|0|20|t
+1.1|0|21|t
+1.1|0|22|t
+1.1|0|23|t
+1.1|0|30|t
+1.1|0|31|t
+1.1|0|32|t
+1.2|0|1|t
+1.2|0|2|t
+1.2|0|10|t
+1.2|0|11|t
+1.2|0|12|t
+1.2|0|13|f
+1.2|0|14|f
+1.2|0|15|t
+1.2|0|16|t
+1.2|0|20|t
+1.2|0|21|t
+1.2|0|22|t
+1.2|0|23|t
+1.2|0|30|t
+1.2|0|31|t
+1.2|0|32|t
+2.1|0|1|ST_Point|t
+2.1|0|2|ST_Point|t
+2.1|0|3|ST_Point|t
+2.1|0|4|ST_Point|t
+2.1|0|5|ST_Point|f
+2.1|0|6|ST_Point|f
+2.1|0|7|ST_Point|f
+2.1|0|8|ST_Point|f
+2.1|0|11|ST_MultiPoint|t
+2.1|0|12|ST_MultiPoint|t
+2.1|0|13|ST_MultiPoint|t
+2.1|0|14|ST_MultiPoint|f
+2.1|0|15|ST_MultiPoint|f
+2.1|0|21|ST_LineString|t
+2.1|0|22|ST_LineString|t
+2.1|0|23|ST_LineString|t
+2.1|0|24|ST_LineString|f
+2.1|0|25|ST_LineString|f
+2.1|0|26|ST_LineString|t
+2.1|0|27|ST_LineString|t
+2.1|0|28|ST_LineString|t
+2.1|0|29|ST_LineString|f
+2.1|0|31|ST_Polygon|t
+2.1|0|32|ST_Polygon|t
+2.1|0|33|ST_Polygon|t
+2.1|0|34|ST_Polygon|t
+2.1|0|35|ST_Polygon|t
+2.1|0|36|ST_Polygon|f
+2.1|0|41|ST_MultiPolygon|t
+2.1|0|42|ST_MultiPolygon|t
+2.1|0|43|ST_MultiPolygon|t
+2.1|0|44|ST_MultiPolygon|t
+2.1|0|45|ST_MultiPolygon|t
+2.1|0|46|ST_MultiPolygon|f
+2.2|0|1|ST_Point|t
+2.2|0|2|ST_Point|t
+2.2|0|3|ST_Point|t
+2.2|0|4|ST_Point|t
+2.2|0|5|ST_Point|f
+2.2|0|6|ST_Point|f
+2.2|0|7|ST_Point|f
+2.2|0|8|ST_Point|f
+2.2|0|11|ST_MultiPoint|t
+2.2|0|12|ST_MultiPoint|t
+2.2|0|13|ST_MultiPoint|t
+2.2|0|14|ST_MultiPoint|f
+2.2|0|15|ST_MultiPoint|f
+2.2|0|21|ST_LineString|t
+2.2|0|22|ST_LineString|t
+2.2|0|23|ST_LineString|t
+2.2|0|24|ST_LineString|f
+2.2|0|25|ST_LineString|f
+2.2|0|26|ST_LineString|t
+2.2|0|27|ST_LineString|t
+2.2|0|28|ST_LineString|t
+2.2|0|29|ST_LineString|f
+2.2|0|31|ST_Polygon|t
+2.2|0|32|ST_Polygon|t
+2.2|0|33|ST_Polygon|t
+2.2|0|34|ST_Polygon|t
+2.2|0|35|ST_Polygon|t
+2.2|0|36|ST_Polygon|f
+2.2|0|41|ST_MultiPolygon|t
+2.2|0|42|ST_MultiPolygon|t
+2.2|0|43|ST_MultiPolygon|t
+2.2|0|44|ST_MultiPolygon|t
+2.2|0|45|ST_MultiPolygon|t
+2.2|0|46|ST_MultiPolygon|f
+2.3|2|1|ST_Point|t
+2.3|2|2|ST_Point|t
+2.3|2|3|ST_Point|f
+2.3|2|4|ST_Point|f
+2.3|2|5|ST_Point|f
+2.3|2|6|ST_Point|f
+2.3|2|7|ST_Point|f
+2.3|2|8|ST_Point|t
+2.3|2|11|ST_MultiPoint|t
+2.3|2|12|ST_MultiPoint|t
+2.3|2|13|ST_MultiPoint|t
+2.3|2|14|ST_MultiPoint|t
+2.3|2|15|ST_MultiPoint|t
+2.3|2|21|ST_LineString|t
+2.3|2|22|ST_LineString|t
+2.3|2|23|ST_LineString|t
+2.3|2|24|ST_LineString|t
+2.3|2|25|ST_LineString|t
+2.3|2|26|ST_LineString|t
+2.3|2|27|ST_LineString|t
+2.3|2|28|ST_LineString|t
+2.3|2|29|ST_LineString|t
+2.3|2|31|ST_Polygon|t
+2.3|2|32|ST_Polygon|t
+2.3|2|33|ST_Polygon|t
+2.3|2|34|ST_Polygon|t
+2.3|2|35|ST_Polygon|t
+2.3|2|36|ST_Polygon|t
+2.3|2|41|ST_MultiPolygon|t
+2.3|2|42|ST_MultiPolygon|t
+2.3|2|43|ST_MultiPolygon|t
+2.3|2|44|ST_MultiPolygon|t
+2.3|2|45|ST_MultiPolygon|t
+2.3|2|46|ST_MultiPolygon|t
+2.4|2|1|ST_Point|t
+2.4|2|2|ST_Point|t
+2.4|2|3|ST_Point|f
+2.4|2|4|ST_Point|f
+2.4|2|5|ST_Point|f
+2.4|2|6|ST_Point|f
+2.4|2|7|ST_Point|f
+2.4|2|8|ST_Point|t
+2.4|2|11|ST_MultiPoint|t
+2.4|2|12|ST_MultiPoint|t
+2.4|2|13|ST_MultiPoint|t
+2.4|2|14|ST_MultiPoint|t
+2.4|2|15|ST_MultiPoint|t
+2.4|2|21|ST_LineString|t
+2.4|2|22|ST_LineString|t
+2.4|2|23|ST_LineString|t
+2.4|2|24|ST_LineString|t
+2.4|2|25|ST_LineString|t
+2.4|2|26|ST_LineString|t
+2.4|2|27|ST_LineString|t
+2.4|2|28|ST_LineString|t
+2.4|2|29|ST_LineString|t
+2.4|2|31|ST_Polygon|t
+2.4|2|32|ST_Polygon|t
+2.4|2|33|ST_Polygon|t
+2.4|2|34|ST_Polygon|t
+2.4|2|35|ST_Polygon|t
+2.4|2|36|ST_Polygon|t
+2.4|2|41|ST_MultiPolygon|t
+2.4|2|42|ST_MultiPolygon|t
+2.4|2|43|ST_MultiPolygon|t
+2.4|2|44|ST_MultiPolygon|t
+2.4|2|45|ST_MultiPolygon|t
+2.4|2|46|ST_MultiPolygon|t
+2.5|0|1|ST_Point|t
+2.5|0|2|ST_Point|t
+2.5|0|3|ST_Point|t
+2.5|0|4|ST_Point|t
+2.5|0|5|ST_Point|f
+2.5|0|6|ST_Point|f
+2.5|0|7|ST_Point|f
+2.5|0|8|ST_Point|f
+2.5|0|11|ST_MultiPoint|t
+2.5|0|12|ST_MultiPoint|t
+2.5|0|13|ST_MultiPoint|t
+2.5|0|14|ST_MultiPoint|f
+2.5|0|15|ST_MultiPoint|f
+2.5|0|21|ST_LineString|t
+2.5|0|22|ST_LineString|t
+2.5|0|23|ST_LineString|t
+2.5|0|24|ST_LineString|f
+2.5|0|25|ST_LineString|f
+2.5|0|26|ST_LineString|t
+2.5|0|27|ST_LineString|t
+2.5|0|28|ST_LineString|t
+2.5|0|29|ST_LineString|f
+2.5|0|31|ST_Polygon|t
+2.5|0|32|ST_Polygon|t
+2.5|0|33|ST_Polygon|t
+2.5|0|34|ST_Polygon|t
+2.5|0|35|ST_Polygon|f
+2.5|0|36|ST_Polygon|f
+2.5|0|41|ST_MultiPolygon|t
+2.5|0|42|ST_MultiPolygon|t
+2.5|0|43|ST_MultiPolygon|t
+2.5|0|44|ST_MultiPolygon|t
+2.5|0|45|ST_MultiPolygon|f
+2.5|0|46|ST_MultiPolygon|f
+2.6|2|1|ST_Point|t
+2.6|2|2|ST_Point|t
+2.6|2|3|ST_Point|f
+2.6|2|4|ST_Point|f
+2.6|2|5|ST_Point|f
+2.6|2|6|ST_Point|f
+2.6|2|7|ST_Point|f
+2.6|2|8|ST_Point|t
+2.6|2|11|ST_MultiPoint|t
+2.6|2|12|ST_MultiPoint|t
+2.6|2|13|ST_MultiPoint|t
+2.6|2|14|ST_MultiPoint|t
+2.6|2|15|ST_MultiPoint|t
+2.6|2|21|ST_LineString|t
+2.6|2|22|ST_LineString|t
+2.6|2|23|ST_LineString|t
+2.6|2|24|ST_LineString|t
+2.6|2|25|ST_LineString|t
+2.6|2|26|ST_LineString|t
+2.6|2|27|ST_LineString|t
+2.6|2|28|ST_LineString|t
+2.6|2|29|ST_LineString|t
+2.6|2|31|ST_Polygon|t
+2.6|2|32|ST_Polygon|t
+2.6|2|33|ST_Polygon|t
+2.6|2|34|ST_Polygon|t
+2.6|2|35|ST_Polygon|t
+2.6|2|36|ST_Polygon|t
+2.6|2|41|ST_MultiPolygon|t
+2.6|2|42|ST_MultiPolygon|t
+2.6|2|43|ST_MultiPolygon|t
+2.6|2|44|ST_MultiPolygon|t
+2.6|2|45|ST_MultiPolygon|t
+2.6|2|46|ST_MultiPolygon|t
-- Test 1 - st_intersect(geometry, raster)
-----------------------------------------------------------------------
-SELECT 'test 1.1', rid, geomid, st_intersects(geom, rast) AS intersect
+SELECT 'test 1.1', rid, geomid, ST_Intersects(geom, rast) AS intersect
FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
WHERE forrast = rid ORDER BY rid, geomid;
-SELECT 'test 1.2', rid, geomid, st_intersects(geom, rast) AS intersect
+SELECT 'test 1.2', rid, geomid, ST_Intersects(geom, rast) AS intersect
FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
WHERE forrast = rid AND ST_Intersects(geom, rast, 1) ORDER BY rid, geomid;
-SELECT 'test 1.3', rid, geomid, st_intersects(geom, st_setbandnodatavalue(rast, NULL))
+SELECT 'test 1.3', rid, geomid, ST_Intersects(geom, st_setbandnodatavalue(rast, NULL), 1)
FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
WHERE forrast = rid ORDER BY rid, geomid;
-SELECT 'test 1.4', rid, geomid, st_intersects(geom, st_setbandnodatavalue(rast, NULL))
+SELECT 'test 1.4', rid, geomid, ST_Intersects(geom, st_setbandnodatavalue(rast, NULL))
FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
WHERE forrast = rid AND ST_Intersects(geom, rast, 1) ORDER BY rid, geomid;
SELECT 'test 2.2', rid, geomid, ST_AsText((gv).geom), (gv).val
FROM (SELECT *, ST_Intersection(geom, rast) gv
FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
- WHERE forrast = rid AND ST_Intersects(geom, rast)
+ WHERE forrast = rid AND ST_Intersects(geom, rast, 1)
) foo;
SELECT 'test 2.3', rid, geomid, ST_AsText((gv).geom), (gv).val
SELECT 'test 2.4', rid, geomid, ST_AsText((gv).geom), (gv).val
FROM (SELECT *, ST_Intersection(geom, st_setbandnodatavalue(rast, NULL)) gv
FROM rt_spatial_relationship_test, rt_spatial_relationship_test_geom
- WHERE forrast = rid AND ST_Intersects(geom, st_setbandnodatavalue(rast, NULL))
+ WHERE forrast = rid AND ST_Intersects(geom, st_setbandnodatavalue(rast, NULL), 1)
) foo;
test 1.1|1|1|t
test 1.1|1|2|t
test 1.1|1|3|f
-test 1.1|1|4|f
+test 1.1|1|4|t
test 1.1|1|5|t
test 1.1|1|6|t
test 1.1|1|7|f
test 1.1|1|11|f
test 1.1|1|12|t
test 1.1|1|13|t
-test 1.1|1|14|f
+test 1.1|1|14|t
test 1.1|1|15|t
-test 1.1|1|16|f
+test 1.1|1|16|t
test 1.1|1|17|t
test 1.1|1|18|t
test 1.1|2|31|t
test 1.1|2|32|t
test 1.1|2|33|f
-test 1.1|2|34|f
-test 1.1|2|35|f
-test 1.1|2|36|f
+test 1.1|2|34|t
+test 1.1|2|35|t
+test 1.1|2|36|t
test 1.1|2|41|f
test 1.1|2|42|t
test 1.1|2|43|t
-test 1.1|2|44|f
+test 1.1|2|44|t
test 1.1|2|45|f
test 1.1|2|46|f
test 1.2|1|1|t
test 1.3|2|33|f
test 1.3|2|34|t
test 1.3|2|35|t
-test 1.3|2|36|f
+test 1.3|2|36|t
test 1.3|2|41|f
test 1.3|2|42|t
test 1.3|2|43|t
test 2.3|2|33|GEOMETRYCOLLECTION EMPTY|
test 2.3|2|34|POINT(-75.5518328537098 49.2814585505576)|4
test 2.3|2|35|POINT(-75.5523150760919 49.2818643977074)|4
-test 2.3|2|36|GEOMETRYCOLLECTION EMPTY|
test 2.3|2|41|GEOMETRYCOLLECTION EMPTY|
test 2.3|2|42|LINESTRING(-75.5533120424461 49.2823793618213,-75.5531972042327 49.2824942000347)|1
test 2.3|2|43|LINESTRING(-75.5529884291587 49.2811479840607,-75.5522356128797 49.2815620330142)|3
projects / postgis / commitdiff