PostGIS  3.0.6dev-r@@SVN_REVISION@@

◆ rt_raster_intersects()

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

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

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

Referenced by RASTER_intersects(), and test_raster_intersects().

Here is the call graph for this function:
Here is the caller graph for this function: