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

Return zero if error occurred in function.

Return ES_ERROR if error occurred in function.

Parameter intersects returns non-zero if two rasters intersect

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

Definition at line 11953 of file rt_api.c.

References ES_ERROR, ES_NONE, FALSE, FLT_EQ, LWGEOM2GEOS(), lwgeom_free(), lwgeom_geos_error(), lwnotice(), 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(), and rterror().

Referenced by RASTER_intersects(), and test_raster_intersects().

11957  {
11958  int i;
11959  int j;
11960  int within = 0;
11961 
11962  LWGEOM *hull[2] = {NULL};
11963  GEOSGeometry *ghull[2] = {NULL};
11964 
11965  uint16_t width1;
11966  uint16_t height1;
11967  uint16_t width2;
11968  uint16_t height2;
11969  double area1;
11970  double area2;
11971  double pixarea1;
11972  double pixarea2;
11973  rt_raster rastS = NULL;
11974  rt_raster rastL = NULL;
11975  uint16_t *widthS = NULL;
11976  uint16_t *heightS = NULL;
11977  uint16_t *widthL = NULL;
11978  uint16_t *heightL = NULL;
11979  int nbandS;
11980  int nbandL;
11981  rt_band bandS = NULL;
11982  rt_band bandL = NULL;
11983  int hasnodataS = FALSE;
11984  int hasnodataL = FALSE;
11985  double nodataS = 0;
11986  double nodataL = 0;
11987  int isnodataS = 0;
11988  int isnodataL = 0;
11989  double gtS[6] = {0};
11990  double igtL[6] = {0};
11991 
11992  uint32_t row;
11993  uint32_t rowoffset;
11994  uint32_t col;
11995  uint32_t coloffset;
11996 
11997  enum line_points {X1, Y1, X2, Y2};
11998  enum point {pX, pY};
11999  double lineS[4];
12000  double Qr[2];
12001  double valS;
12002  double valL;
12003 
12004  RASTER_DEBUG(3, "Starting");
12005 
12006  assert(NULL != rast1);
12007  assert(NULL != rast2);
12008  assert(NULL != intersects);
12009 
12010  if (nband1 < 0 && nband2 < 0) {
12011  nband1 = -1;
12012  nband2 = -1;
12013  }
12014  else {
12015  assert(nband1 >= 0 && nband1 < rt_raster_get_num_bands(rast1));
12016  assert(nband2 >= 0 && nband2 < rt_raster_get_num_bands(rast2));
12017  }
12018 
12019  /* same srid */
12020  if (rt_raster_get_srid(rast1) != rt_raster_get_srid(rast2)) {
12021  rterror("rt_raster_intersects: The two rasters provided have different SRIDs");
12022  *intersects = 0;
12023  return ES_ERROR;
12024  }
12025 
12026  /* raster extents need to intersect */
12027  do {
12028  int rtn;
12029 
12030  initGEOS(lwnotice, lwgeom_geos_error);
12031 
12032  rtn = 1;
12033  for (i = 0; i < 2; i++) {
12034  if ((rt_raster_get_convex_hull(i < 1 ? rast1 : rast2, &(hull[i])) != ES_NONE) || NULL == hull[i]) {
12035  for (j = 0; j < i; j++) {
12036  GEOSGeom_destroy(ghull[j]);
12037  lwgeom_free(hull[j]);
12038  }
12039  rtn = 0;
12040  break;
12041  }
12042  ghull[i] = (GEOSGeometry *) LWGEOM2GEOS(hull[i]);
12043  if (NULL == ghull[i]) {
12044  for (j = 0; j < i; j++) {
12045  GEOSGeom_destroy(ghull[j]);
12046  lwgeom_free(hull[j]);
12047  }
12048  lwgeom_free(hull[i]);
12049  rtn = 0;
12050  break;
12051  }
12052  }
12053  if (!rtn) break;
12054 
12055  /* test to see if raster within the other */
12056  within = 0;
12057  if (GEOSWithin(ghull[0], ghull[1]) == 1)
12058  within = -1;
12059  else if (GEOSWithin(ghull[1], ghull[0]) == 1)
12060  within = 1;
12061 
12062  if (within != 0)
12063  rtn = 1;
12064  else
12065  rtn = GEOSIntersects(ghull[0], ghull[1]);
12066 
12067  for (i = 0; i < 2; i++) {
12068  GEOSGeom_destroy(ghull[i]);
12069  lwgeom_free(hull[i]);
12070  }
12071 
12072  if (rtn != 2) {
12073  RASTER_DEBUGF(4, "convex hulls of rasters do %sintersect", rtn != 1 ? "NOT " : "");
12074  if (rtn != 1) {
12075  *intersects = 0;
12076  return ES_NONE;
12077  }
12078  /* band isn't specified */
12079  else if (nband1 < 0) {
12080  *intersects = 1;
12081  return ES_NONE;
12082  }
12083  }
12084  else {
12085  RASTER_DEBUG(4, "GEOSIntersects() returned a 2!!!!");
12086  }
12087  }
12088  while (0);
12089 
12090  /* smaller raster by area or width */
12091  width1 = rt_raster_get_width(rast1);
12092  height1 = rt_raster_get_height(rast1);
12093  width2 = rt_raster_get_width(rast2);
12094  height2 = rt_raster_get_height(rast2);
12095  pixarea1 = fabs(rt_raster_get_x_scale(rast1) * rt_raster_get_y_scale(rast1));
12096  pixarea2 = fabs(rt_raster_get_x_scale(rast2) * rt_raster_get_y_scale(rast2));
12097  area1 = fabs(width1 * height1 * pixarea1);
12098  area2 = fabs(width2 * height2 * pixarea2);
12099  RASTER_DEBUGF(4, "pixarea1, pixarea2, area1, area2 = %f, %f, %f, %f",
12100  pixarea1, pixarea2, area1, area2);
12101  if (
12102  (within <= 0) ||
12103  (area1 < area2) ||
12104  FLT_EQ(area1, area2) ||
12105  (area1 < pixarea2) || /* area of rast1 smaller than pixel area of rast2 */
12106  FLT_EQ(area1, pixarea2)
12107  ) {
12108  rastS = rast1;
12109  nbandS = nband1;
12110  widthS = &width1;
12111  heightS = &height1;
12112 
12113  rastL = rast2;
12114  nbandL = nband2;
12115  widthL = &width2;
12116  heightL = &height2;
12117  }
12118  else {
12119  rastS = rast2;
12120  nbandS = nband2;
12121  widthS = &width2;
12122  heightS = &height2;
12123 
12124  rastL = rast1;
12125  nbandL = nband1;
12126  widthL = &width1;
12127  heightL = &height1;
12128  }
12129 
12130  /* no band to use, set band to zero */
12131  if (nband1 < 0) {
12132  nbandS = 0;
12133  nbandL = 0;
12134  }
12135 
12136  RASTER_DEBUGF(4, "rast1 @ %p", rast1);
12137  RASTER_DEBUGF(4, "rast2 @ %p", rast2);
12138  RASTER_DEBUGF(4, "rastS @ %p", rastS);
12139  RASTER_DEBUGF(4, "rastL @ %p", rastL);
12140 
12141  /* load band of smaller raster */
12142  bandS = rt_raster_get_band(rastS, nbandS);
12143  if (NULL == bandS) {
12144  rterror("rt_raster_intersects: Could not get band %d of the first raster", nbandS);
12145  *intersects = 0;
12146  return ES_ERROR;
12147  }
12148 
12149  hasnodataS = rt_band_get_hasnodata_flag(bandS);
12150  if (hasnodataS != FALSE)
12151  rt_band_get_nodata(bandS, &nodataS);
12152 
12153  /* load band of larger raster */
12154  bandL = rt_raster_get_band(rastL, nbandL);
12155  if (NULL == bandL) {
12156  rterror("rt_raster_intersects: Could not get band %d of the first raster", nbandL);
12157  *intersects = 0;
12158  return ES_ERROR;
12159  }
12160 
12161  hasnodataL = rt_band_get_hasnodata_flag(bandL);
12162  if (hasnodataL != FALSE)
12163  rt_band_get_nodata(bandL, &nodataL);
12164 
12165  /* no band to use, ignore nodata */
12166  if (nband1 < 0) {
12167  hasnodataS = FALSE;
12168  hasnodataL = FALSE;
12169  }
12170 
12171  /* hasnodata(S|L) = TRUE and one of the two bands is isnodata */
12172  if (
12173  (hasnodataS && rt_band_get_isnodata_flag(bandS)) ||
12174  (hasnodataL && rt_band_get_isnodata_flag(bandL))
12175  ) {
12176  RASTER_DEBUG(3, "One of the two raster bands is NODATA. The two rasters do not intersect");
12177  *intersects = 0;
12178  return ES_NONE;
12179  }
12180 
12181  /* special case where a raster can fit inside another raster's pixel */
12182  if (within != 0 && ((pixarea1 > area2) || (pixarea2 > area1))) {
12183  RASTER_DEBUG(4, "Using special case of raster fitting into another raster's pixel");
12184  /* 3 x 3 search */
12185  for (coloffset = 0; coloffset < 3; coloffset++) {
12186  for (rowoffset = 0; rowoffset < 3; rowoffset++) {
12187  for (col = coloffset; col < *widthS; col += 3) {
12188  for (row = rowoffset; row < *heightS; row += 3) {
12189  if (hasnodataS == FALSE)
12190  valS = 1;
12191  else if (rt_band_get_pixel(bandS, col, row, &valS, &isnodataS) != ES_NONE)
12192  continue;
12193 
12194  if ((hasnodataS == FALSE) || !isnodataS) {
12196  rastS,
12197  col, row,
12198  &(lineS[X1]), &(lineS[Y1]),
12199  gtS
12200  );
12201 
12203  rastL,
12204  lineS[X1], lineS[Y1],
12205  &(Qr[pX]), &(Qr[pY]),
12206  igtL
12207  ) != ES_NONE) {
12208  continue;
12209  }
12210 
12211  if (
12212  (Qr[pX] < 0 || Qr[pX] > *widthL || FLT_EQ(Qr[pX], *widthL)) ||
12213  (Qr[pY] < 0 || Qr[pY] > *heightL || FLT_EQ(Qr[pY], *heightL))
12214  ) {
12215  continue;
12216  }
12217 
12218  if (hasnodataS == FALSE)
12219  valL = 1;
12220  else if (rt_band_get_pixel(bandL, Qr[pX], Qr[pY], &valL, &isnodataL) != ES_NONE)
12221  continue;
12222 
12223  if ((hasnodataL == FALSE) || !isnodataL) {
12224  RASTER_DEBUG(3, "The two rasters do intersect");
12225  *intersects = 1;
12226  return ES_NONE;
12227  }
12228  }
12229  }
12230  }
12231  }
12232  }
12233  RASTER_DEBUG(4, "Smaller raster not in the other raster's pixel. Continuing");
12234  }
12235 
12236  RASTER_DEBUG(4, "Testing smaller raster vs larger raster");
12238  rastS, rastL,
12239  bandS, bandL,
12240  hasnodataS, hasnodataL,
12241  nodataS, nodataL
12242  );
12243 
12244  if (*intersects) return ES_NONE;
12245 
12246  RASTER_DEBUG(4, "Testing larger raster vs smaller raster");
12248  rastL, rastS,
12249  bandL, bandS,
12250  hasnodataL, hasnodataS,
12251  nodataL, nodataS
12252  );
12253 
12254  if (*intersects) return ES_NONE;
12255 
12256  RASTER_DEBUG(3, "The two rasters do not intersect");
12257 
12258  *intersects = 0;
12259  return ES_NONE;
12260 }
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_api.c:5677
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1006
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom)
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_api.c:5661
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_api.c:2042
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_api.c:6054
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:54
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
void lwgeom_geos_error(const char *fmt,...)
Datum intersects(PG_FUNCTION_ARGS)
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_api.c:3058
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_api.c:2549
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_api.c:5455
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_api.c:5434
void rterror(const char *fmt,...)
Raster core error and info handlers.
Definition: rt_api.c:895
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
#define FALSE
Definition: dbfopen.c:169
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_api.c:5464
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_api.c:6105
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition: rt_api.c:6556
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_api.c:11598
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_api.c:2002
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_api.c:5426

Here is the call graph for this function:

Here is the caller graph for this function: