PostGIS  2.2.7dev-r@@SVN_REVISION@@
rt_errorstate rt_raster_intersects ( rt_raster  rast1,
int  nband1,
rt_raster  rast2,
int  nband2,
int *  intersects 
)

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

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 992 of file rt_spatial_relationship.c.

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(), and rtinfo().

Referenced by RASTER_intersects(), and test_raster_intersects().

996  {
997  int i;
998  int j;
999  int within = 0;
1000 
1001  LWGEOM *hull[2] = {NULL};
1002  GEOSGeometry *ghull[2] = {NULL};
1003 
1004  uint16_t width1;
1005  uint16_t height1;
1006  uint16_t width2;
1007  uint16_t height2;
1008  double area1;
1009  double area2;
1010  double pixarea1;
1011  double pixarea2;
1012  rt_raster rastS = NULL;
1013  rt_raster rastL = NULL;
1014  uint16_t *widthS = NULL;
1015  uint16_t *heightS = NULL;
1016  uint16_t *widthL = NULL;
1017  uint16_t *heightL = NULL;
1018  int nbandS;
1019  int nbandL;
1020  rt_band bandS = NULL;
1021  rt_band bandL = NULL;
1022  int hasnodataS = FALSE;
1023  int hasnodataL = FALSE;
1024  double nodataS = 0;
1025  double nodataL = 0;
1026  int isnodataS = 0;
1027  int isnodataL = 0;
1028  double gtS[6] = {0};
1029  double igtL[6] = {0};
1030 
1031  uint32_t row;
1032  uint32_t rowoffset;
1033  uint32_t col;
1034  uint32_t coloffset;
1035 
1036  enum line_points {X1, Y1, X2, Y2};
1037  enum point {pX, pY};
1038  double lineS[4];
1039  double Qr[2];
1040  double valS;
1041  double valL;
1042 
1043  RASTER_DEBUG(3, "Starting");
1044 
1045  assert(NULL != rast1);
1046  assert(NULL != rast2);
1047  assert(NULL != intersects);
1048 
1049  if (nband1 < 0 && nband2 < 0) {
1050  nband1 = -1;
1051  nband2 = -1;
1052  }
1053  else {
1054  assert(nband1 >= 0 && nband1 < rt_raster_get_num_bands(rast1));
1055  assert(nband2 >= 0 && nband2 < rt_raster_get_num_bands(rast2));
1056  }
1057 
1058  /* same srid */
1059  if (rt_raster_get_srid(rast1) != rt_raster_get_srid(rast2)) {
1060  rterror("rt_raster_intersects: The two rasters provided have different SRIDs");
1061  *intersects = 0;
1062  return ES_ERROR;
1063  }
1064 
1065  /* raster extents need to intersect */
1066  do {
1067  int rtn;
1068 
1069  initGEOS(rtinfo, lwgeom_geos_error);
1070 
1071  rtn = 1;
1072  for (i = 0; i < 2; i++) {
1073  if ((rt_raster_get_convex_hull(i < 1 ? rast1 : rast2, &(hull[i])) != ES_NONE) || NULL == hull[i]) {
1074  for (j = 0; j < i; j++) {
1075  GEOSGeom_destroy(ghull[j]);
1076  lwgeom_free(hull[j]);
1077  }
1078  rtn = 0;
1079  break;
1080  }
1081  ghull[i] = (GEOSGeometry *) LWGEOM2GEOS(hull[i], 0);
1082  if (NULL == ghull[i]) {
1083  for (j = 0; j < i; j++) {
1084  GEOSGeom_destroy(ghull[j]);
1085  lwgeom_free(hull[j]);
1086  }
1087  lwgeom_free(hull[i]);
1088  rtn = 0;
1089  break;
1090  }
1091  }
1092  if (!rtn) break;
1093 
1094  /* test to see if raster within the other */
1095  within = 0;
1096  if (GEOSWithin(ghull[0], ghull[1]) == 1)
1097  within = -1;
1098  else if (GEOSWithin(ghull[1], ghull[0]) == 1)
1099  within = 1;
1100 
1101  if (within != 0)
1102  rtn = 1;
1103  else
1104  rtn = GEOSIntersects(ghull[0], ghull[1]);
1105 
1106  for (i = 0; i < 2; i++) {
1107  GEOSGeom_destroy(ghull[i]);
1108  lwgeom_free(hull[i]);
1109  }
1110 
1111  if (rtn != 2) {
1112  RASTER_DEBUGF(4, "convex hulls of rasters do %sintersect", rtn != 1 ? "NOT " : "");
1113  if (rtn != 1) {
1114  *intersects = 0;
1115  return ES_NONE;
1116  }
1117  /* band isn't specified */
1118  else if (nband1 < 0) {
1119  *intersects = 1;
1120  return ES_NONE;
1121  }
1122  }
1123  else {
1124  RASTER_DEBUG(4, "GEOSIntersects() returned a 2!!!!");
1125  }
1126  }
1127  while (0);
1128 
1129  /* smaller raster by area or width */
1130  width1 = rt_raster_get_width(rast1);
1131  height1 = rt_raster_get_height(rast1);
1132  width2 = rt_raster_get_width(rast2);
1133  height2 = rt_raster_get_height(rast2);
1134  pixarea1 = fabs(rt_raster_get_x_scale(rast1) * rt_raster_get_y_scale(rast1));
1135  pixarea2 = fabs(rt_raster_get_x_scale(rast2) * rt_raster_get_y_scale(rast2));
1136  area1 = fabs(width1 * height1 * pixarea1);
1137  area2 = fabs(width2 * height2 * pixarea2);
1138  RASTER_DEBUGF(4, "pixarea1, pixarea2, area1, area2 = %f, %f, %f, %f",
1139  pixarea1, pixarea2, area1, area2);
1140  if (
1141  (within <= 0) ||
1142  (area1 < area2) ||
1143  FLT_EQ(area1, area2) ||
1144  (area1 < pixarea2) || /* area of rast1 smaller than pixel area of rast2 */
1145  FLT_EQ(area1, pixarea2)
1146  ) {
1147  rastS = rast1;
1148  nbandS = nband1;
1149  widthS = &width1;
1150  heightS = &height1;
1151 
1152  rastL = rast2;
1153  nbandL = nband2;
1154  widthL = &width2;
1155  heightL = &height2;
1156  }
1157  else {
1158  rastS = rast2;
1159  nbandS = nband2;
1160  widthS = &width2;
1161  heightS = &height2;
1162 
1163  rastL = rast1;
1164  nbandL = nband1;
1165  widthL = &width1;
1166  heightL = &height1;
1167  }
1168 
1169  /* no band to use, set band to zero */
1170  if (nband1 < 0) {
1171  nbandS = 0;
1172  nbandL = 0;
1173  }
1174 
1175  RASTER_DEBUGF(4, "rast1 @ %p", rast1);
1176  RASTER_DEBUGF(4, "rast2 @ %p", rast2);
1177  RASTER_DEBUGF(4, "rastS @ %p", rastS);
1178  RASTER_DEBUGF(4, "rastL @ %p", rastL);
1179 
1180  /* load band of smaller raster */
1181  bandS = rt_raster_get_band(rastS, nbandS);
1182  if (NULL == bandS) {
1183  rterror("rt_raster_intersects: Could not get band %d of the first raster", nbandS);
1184  *intersects = 0;
1185  return ES_ERROR;
1186  }
1187 
1188  hasnodataS = rt_band_get_hasnodata_flag(bandS);
1189  if (hasnodataS != FALSE)
1190  rt_band_get_nodata(bandS, &nodataS);
1191 
1192  /* load band of larger raster */
1193  bandL = rt_raster_get_band(rastL, nbandL);
1194  if (NULL == bandL) {
1195  rterror("rt_raster_intersects: Could not get band %d of the first raster", nbandL);
1196  *intersects = 0;
1197  return ES_ERROR;
1198  }
1199 
1200  hasnodataL = rt_band_get_hasnodata_flag(bandL);
1201  if (hasnodataL != FALSE)
1202  rt_band_get_nodata(bandL, &nodataL);
1203 
1204  /* no band to use, ignore nodata */
1205  if (nband1 < 0) {
1206  hasnodataS = FALSE;
1207  hasnodataL = FALSE;
1208  }
1209 
1210  /* hasnodata(S|L) = TRUE and one of the two bands is isnodata */
1211  if (
1212  (hasnodataS && rt_band_get_isnodata_flag(bandS)) ||
1213  (hasnodataL && rt_band_get_isnodata_flag(bandL))
1214  ) {
1215  RASTER_DEBUG(3, "One of the two raster bands is NODATA. The two rasters do not intersect");
1216  *intersects = 0;
1217  return ES_NONE;
1218  }
1219 
1220  /* special case where a raster can fit inside another raster's pixel */
1221  if (within != 0 && ((pixarea1 > area2) || (pixarea2 > area1))) {
1222  RASTER_DEBUG(4, "Using special case of raster fitting into another raster's pixel");
1223  /* 3 x 3 search */
1224  for (coloffset = 0; coloffset < 3; coloffset++) {
1225  for (rowoffset = 0; rowoffset < 3; rowoffset++) {
1226  for (col = coloffset; col < *widthS; col += 3) {
1227  for (row = rowoffset; row < *heightS; row += 3) {
1228  if (hasnodataS == FALSE)
1229  valS = 1;
1230  else if (rt_band_get_pixel(bandS, col, row, &valS, &isnodataS) != ES_NONE)
1231  continue;
1232 
1233  if ((hasnodataS == FALSE) || !isnodataS) {
1235  rastS,
1236  col, row,
1237  &(lineS[X1]), &(lineS[Y1]),
1238  gtS
1239  );
1240 
1242  rastL,
1243  lineS[X1], lineS[Y1],
1244  &(Qr[pX]), &(Qr[pY]),
1245  igtL
1246  ) != ES_NONE) {
1247  continue;
1248  }
1249 
1250  if (
1251  (Qr[pX] < 0 || Qr[pX] > *widthL || FLT_EQ(Qr[pX], *widthL)) ||
1252  (Qr[pY] < 0 || Qr[pY] > *heightL || FLT_EQ(Qr[pY], *heightL))
1253  ) {
1254  continue;
1255  }
1256 
1257  if (hasnodataS == FALSE)
1258  valL = 1;
1259  else if (rt_band_get_pixel(bandL, Qr[pX], Qr[pY], &valL, &isnodataL) != ES_NONE)
1260  continue;
1261 
1262  if ((hasnodataL == FALSE) || !isnodataL) {
1263  RASTER_DEBUG(3, "The two rasters do intersect");
1264  *intersects = 1;
1265  return ES_NONE;
1266  }
1267  }
1268  }
1269  }
1270  }
1271  }
1272  RASTER_DEBUG(4, "Smaller raster not in the other raster's pixel. Continuing");
1273  }
1274 
1275  RASTER_DEBUG(4, "Testing smaller raster vs larger raster");
1277  rastS, rastL,
1278  bandS, bandL,
1279  hasnodataS, hasnodataL,
1280  nodataS, nodataL
1281  );
1282 
1283  if (*intersects) return ES_NONE;
1284 
1285  RASTER_DEBUG(4, "Testing larger raster vs smaller raster");
1287  rastL, rastS,
1288  bandL, bandS,
1289  hasnodataL, hasnodataS,
1290  nodataL, nodataS
1291  );
1292 
1293  if (*intersects) return ES_NONE;
1294 
1295  RASTER_DEBUG(3, "The two rasters do not intersect");
1296 
1297  *intersects = 0;
1298  return ES_NONE;
1299 }
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1050
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)
#define FLT_EQ(x, y)
Definition: librtcore.h:2197
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition: rt_geometry.c:822
rt_errorstate rt_raster_cell_to_geopoint(rt_raster raster, double xr, double yr, double *xw, double *yw, double *gt)
Convert an xr, yr raster point to an xw, yw point on map.
Definition: rt_raster.c:755
void lwgeom_geos_error(const char *fmt,...)
Datum intersects(PG_FUNCTION_ARGS)
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1088
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
void rtinfo(const char *fmt,...)
Definition: rt_context.c:211
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:311
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:356
rt_errorstate rt_raster_geopoint_to_cell(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw, yw map point to a xr, yr raster point.
Definition: rt_raster.c:806
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
#define FALSE
Definition: dbfopen.c:168
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:307
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:581
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129

Here is the call graph for this function:

Here is the caller graph for this function: