PostGIS  3.7.0dev-r@@SVN_REVISION@@
rt_spatial_relationship.c
Go to the documentation of this file.
1 /*
2  *
3  * WKTRaster - Raster Types for PostGIS
4  * http://trac.osgeo.org/postgis/wiki/WKTRaster
5  *
6  * Copyright (C) 2011-2013 Regents of the University of California
7  * <bkpark@ucdavis.edu>
8  * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
9  * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
10  * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
11  * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
12  * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
13  *
14  * This program is free software; you can redistribute it and/or
15  * modify it under the terms of the GNU General Public License
16  * as published by the Free Software Foundation; either version 2
17  * of the License, or (at your option) any later version.
18  *
19  * This program is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with this program; if not, write to the Free Software Foundation,
26  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
27  *
28  */
29 
30 #include "librtcore.h"
31 #include "librtcore_internal.h"
32 
33 /*
34  * Return ES_ERROR if error occurred in function.
35  * Parameter aligned returns non-zero if two rasters are aligned
36  *
37  * @param rast1 : the first raster for alignment test
38  * @param rast2 : the second raster for alignment test
39  * @param *aligned : non-zero value if the two rasters are aligned
40  * @param *reason : reason why rasters are not aligned
41  *
42  * @return ES_NONE if success, ES_ERROR if error
43  */
44 rt_errorstate
45 rt_raster_same_alignment(
46  rt_raster rast1,
47  rt_raster rast2,
48  int *aligned, char **reason
49 ) {
50  double xr;
51  double yr;
52  double xw;
53  double yw;
54  int err = 0;
55 
56  assert(NULL != rast1);
57  assert(NULL != rast2);
58  assert(NULL != aligned);
59 
60  err = 0;
61  /* same srid */
62  if (rt_raster_get_srid(rast1) != rt_raster_get_srid(rast2)) {
63  if (reason != NULL) *reason = "The rasters have different SRIDs";
64  err = 1;
65  }
66  /* scales must match */
67  else if (FLT_NEQ(fabs(rast1->scaleX), fabs(rast2->scaleX))) {
68  if (reason != NULL) *reason = "The rasters have different scales on the X axis";
69  err = 1;
70  }
71  else if (FLT_NEQ(fabs(rast1->scaleY), fabs(rast2->scaleY))) {
72  if (reason != NULL) *reason = "The rasters have different scales on the Y axis";
73  err = 1;
74  }
75  /* skews must match */
76  else if (FLT_NEQ(rast1->skewX, rast2->skewX)) {
77  if (reason != NULL) *reason = "The rasters have different skews on the X axis";
78  err = 1;
79  }
80  else if (FLT_NEQ(rast1->skewY, rast2->skewY)) {
81  if (reason != NULL) *reason = "The rasters have different skews on the Y axis";
82  err = 1;
83  }
84 
85  if (err) {
86  *aligned = 0;
87  return ES_NONE;
88  }
89 
90  /* raster coordinates in context of second raster of first raster's upper-left corner */
91  if (rt_raster_geopoint_to_cell(
92  rast2,
93  rast1->ipX, rast1->ipY,
94  &xr, &yr,
95  NULL
96  ) != ES_NONE) {
97  rterror("rt_raster_same_alignment: Could not get raster coordinates of second raster from first raster's spatial coordinates");
98  *aligned = 0;
99  return ES_ERROR;
100  }
101 
102  /* spatial coordinates of raster coordinates from above */
103  if (rt_raster_cell_to_geopoint(
104  rast2,
105  xr, yr,
106  &xw, &yw,
107  NULL
108  ) != ES_NONE) {
109  rterror("rt_raster_same_alignment: Could not get spatial coordinates of second raster from raster coordinates");
110  *aligned = 0;
111  return ES_ERROR;
112  }
113 
114  RASTER_DEBUGF(4, "rast1(ipX, ipxY) = (%f, %f)", rast1->ipX, rast1->ipY);
115  RASTER_DEBUGF(4, "rast2(xr, yr) = (%f, %f)", xr, yr);
116  RASTER_DEBUGF(4, "rast2(xw, yw) = (%f, %f)", xw, yw);
117 
118  /* spatial coordinates are identical to that of first raster's upper-left corner */
119  if (FLT_EQ(xw, rast1->ipX) && FLT_EQ(yw, rast1->ipY)) {
120  if (reason != NULL) *reason = "The rasters are aligned";
121  *aligned = 1;
122  return ES_NONE;
123  }
124 
125  /* no alignment */
126  if (reason != NULL) *reason = "The rasters (pixel corner coordinates) are not aligned";
127 
128  *aligned = 0;
129  return ES_NONE;
130 }
131 
132 /******************************************************************************
133 * GEOS-based spatial relationship tests
134 ******************************************************************************/
135 
136 static
137 rt_errorstate rt_raster_geos_spatial_relationship(
138  rt_raster rast1, int nband1,
139  rt_raster rast2, int nband2,
140  rt_geos_spatial_test testtype,
141  int *testresult
142 ) {
143  LWMPOLY *surface1 = NULL;
144  LWMPOLY *surface2 = NULL;
145  GEOSGeometry *geom1 = NULL;
146  GEOSGeometry *geom2 = NULL;
147  int rtn = 0;
148  int flag = 0;
149 
150  RASTER_DEBUG(3, "Starting");
151 
152  assert(NULL != rast1);
153  assert(NULL != rast2);
154  assert(NULL != testresult);
155 
156  if (nband1 < 0 && nband2 < 0) {
157  nband1 = -1;
158  nband2 = -1;
159  }
160  else {
161  assert(nband1 >= 0 && nband1 < rt_raster_get_num_bands(rast1));
162  assert(nband2 >= 0 && nband2 < rt_raster_get_num_bands(rast2));
163  }
164 
165  /* initialize to zero, false result of spatial relationship test */
166  *testresult = 0;
167 
168  /* same srid */
169  if (rt_raster_get_srid(rast1) != rt_raster_get_srid(rast2)) {
170  rterror("rt_raster_geos_spatial_relationship: The two rasters provided have different SRIDs");
171  return ES_ERROR;
172  }
173 
174  initGEOS(rtinfo, lwgeom_geos_error);
175 
176  /* get LWMPOLY of each band */
177  if (rt_raster_surface(rast1, nband1, &surface1) != ES_NONE) {
178  rterror("rt_raster_geos_spatial_relationship: Could not get surface of the specified band from the first raster");
179  return ES_ERROR;
180  }
181  if (rt_raster_surface(rast2, nband2, &surface2) != ES_NONE) {
182  rterror("rt_raster_geos_spatial_relationship: Could not get surface of the specified band from the second raster");
183  lwmpoly_free(surface1);
184  return ES_ERROR;
185  }
186 
187  /* either surface is NULL, spatial relationship test is false */
188  if (surface1 == NULL || surface2 == NULL) {
189  if (surface1 != NULL) lwmpoly_free(surface1);
190  if (surface2 != NULL) lwmpoly_free(surface2);
191  return ES_NONE;
192  }
193 
194  /* convert LWMPOLY to GEOSGeometry */
195  geom1 = LWGEOM2GEOS(lwmpoly_as_lwgeom(surface1), 0);
196  lwmpoly_free(surface1);
197  if (geom1 == NULL) {
198  rterror("rt_raster_geos_spatial_relationship: Could not convert surface of the specified band from the first raster to a GEOSGeometry");
199  lwmpoly_free(surface2);
200  return ES_ERROR;
201  }
202 
203  geom2 = LWGEOM2GEOS(lwmpoly_as_lwgeom(surface2), 0);
204  lwmpoly_free(surface2);
205  if (geom2 == NULL) {
206  rterror("rt_raster_geos_spatial_relationship: Could not convert surface of the specified band from the second raster to a GEOSGeometry");
207  return ES_ERROR;
208  }
209 
210  flag = 0;
211  switch (testtype) {
212  case GSR_OVERLAPS:
213  rtn = GEOSOverlaps(geom1, geom2);
214  break;
215  case GSR_TOUCHES:
216  rtn = GEOSTouches(geom1, geom2);
217  break;
218  case GSR_CONTAINS:
219  rtn = GEOSContains(geom1, geom2);
220  break;
221  case GSR_CONTAINSPROPERLY:
222  rtn = GEOSRelatePattern(geom1, geom2, "T**FF*FF*");
223  break;
224  case GSR_COVERS:
225  rtn = GEOSRelatePattern(geom1, geom2, "******FF*");
226  break;
227  case GSR_COVEREDBY:
228  rtn = GEOSRelatePattern(geom1, geom2, "**F**F***");
229  break;
230  default:
231  rterror("rt_raster_geos_spatial_relationship: Unknown or unsupported GEOS spatial relationship test");
232  flag = -1;
233  break;
234  }
235  GEOSGeom_destroy(geom1);
236  GEOSGeom_destroy(geom2);
237 
238  /* something happened in the spatial relationship test */
239  if (rtn == 2) {
240  rterror("rt_raster_geos_spatial_relationship: Could not run the appropriate GEOS spatial relationship test");
241  flag = ES_ERROR;
242  }
243  /* spatial relationship test ran fine */
244  else if (flag >= 0) {
245  if (rtn != 0)
246  *testresult = 1;
247  flag = ES_NONE;
248  }
249  /* flag < 0 for when testtype is unknown */
250  else
251  flag = ES_ERROR;
252 
253  return flag;
254 }
255 
272 rt_errorstate rt_raster_overlaps(
273  rt_raster rast1, int nband1,
274  rt_raster rast2, int nband2,
275  int *overlaps
276 ) {
277  RASTER_DEBUG(3, "Starting");
278 
279  return rt_raster_geos_spatial_relationship(
280  rast1, nband1,
281  rast2, nband2,
282  GSR_OVERLAPS,
283  overlaps
284  );
285 }
286 
303 rt_errorstate rt_raster_touches(
304  rt_raster rast1, int nband1,
305  rt_raster rast2, int nband2,
306  int *touches
307 ) {
308  RASTER_DEBUG(3, "Starting");
309 
310  return rt_raster_geos_spatial_relationship(
311  rast1, nband1,
312  rast2, nband2,
313  GSR_TOUCHES,
314  touches
315  );
316 }
317 
334 rt_errorstate rt_raster_contains(
335  rt_raster rast1, int nband1,
336  rt_raster rast2, int nband2,
337  int *contains
338 ) {
339  RASTER_DEBUG(3, "Starting");
340 
341  return rt_raster_geos_spatial_relationship(
342  rast1, nband1,
343  rast2, nband2,
344  GSR_CONTAINS,
345  contains
346  );
347 }
348 
365 rt_errorstate rt_raster_contains_properly(
366  rt_raster rast1, int nband1,
367  rt_raster rast2, int nband2,
368  int *contains
369 ) {
370  RASTER_DEBUG(3, "Starting");
371 
372  return rt_raster_geos_spatial_relationship(
373  rast1, nband1,
374  rast2, nband2,
375  GSR_CONTAINSPROPERLY,
376  contains
377  );
378 }
379 
396 rt_errorstate rt_raster_covers(
397  rt_raster rast1, int nband1,
398  rt_raster rast2, int nband2,
399  int *covers
400 ) {
401  RASTER_DEBUG(3, "Starting");
402 
403  return rt_raster_geos_spatial_relationship(
404  rast1, nband1,
405  rast2, nband2,
406  GSR_COVERS,
407  covers
408  );
409 }
410 
427 rt_errorstate rt_raster_coveredby(
428  rt_raster rast1, int nband1,
429  rt_raster rast2, int nband2,
430  int *coveredby
431 ) {
432  RASTER_DEBUG(3, "Starting");
433 
434  return rt_raster_geos_spatial_relationship(
435  rast1, nband1,
436  rast2, nband2,
437  GSR_COVEREDBY,
438  coveredby
439  );
440 }
441 
460 rt_errorstate rt_raster_within_distance(
461  rt_raster rast1, int nband1,
462  rt_raster rast2, int nband2,
463  double distance,
464  int *dwithin
465 ) {
466  LWMPOLY *surface = NULL;
467  LWGEOM *surface1 = NULL;
468  LWGEOM *surface2 = NULL;
469  double mindist = 0;
470 
471  RASTER_DEBUG(3, "Starting");
472 
473  assert(NULL != rast1);
474  assert(NULL != rast2);
475  assert(NULL != dwithin);
476 
477  if (nband1 < 0 && nband2 < 0) {
478  nband1 = -1;
479  nband2 = -1;
480  }
481  else {
482  assert(nband1 >= 0 && nband1 < rt_raster_get_num_bands(rast1));
483  assert(nband2 >= 0 && nband2 < rt_raster_get_num_bands(rast2));
484  }
485 
486  /* initialize to zero, false result */
487  *dwithin = 0;
488 
489  /* same srid */
490  if (rt_raster_get_srid(rast1) != rt_raster_get_srid(rast2)) {
491  rterror("rt_raster_distance_within: The two rasters provided have different SRIDs");
492  return ES_ERROR;
493  }
494 
495  /* distance cannot be less than zero */
496  if (distance < 0) {
497  rterror("rt_raster_distance_within: Distance cannot be less than zero");
498  return ES_ERROR;
499  }
500 
501  /* get LWMPOLY of each band */
502  if (rt_raster_surface(rast1, nband1, &surface) != ES_NONE) {
503  rterror("rt_raster_distance_within: Could not get surface of the specified band from the first raster");
504  return ES_ERROR;
505  }
506  surface1 = lwmpoly_as_lwgeom(surface);
507 
508  if (rt_raster_surface(rast2, nband2, &surface) != ES_NONE) {
509  rterror("rt_raster_distance_within: Could not get surface of the specified band from the second raster");
510  lwgeom_free(surface1);
511  return ES_ERROR;
512  }
513  surface2 = lwmpoly_as_lwgeom(surface);
514 
515  /* either surface is NULL, test is false */
516  if (surface1 == NULL || surface2 == NULL) {
517  if (surface1 != NULL) lwgeom_free(surface1);
518  if (surface2 != NULL) lwgeom_free(surface2);
519  return ES_NONE;
520  }
521 
522  /* get the min distance between the two surfaces */
523  mindist = lwgeom_mindistance2d_tolerance(surface1, surface2, distance);
524 
525  lwgeom_free(surface1);
526  lwgeom_free(surface2);
527 
528  /* if distance >= mindist, true */
529  if (FLT_EQ(mindist, distance) || distance > mindist)
530  *dwithin = 1;
531 
532  RASTER_DEBUGF(3, "(mindist, distance) = (%f, %f, %d)", mindist, distance, *dwithin);
533 
534  return ES_NONE;
535 }
536 
555 rt_errorstate rt_raster_fully_within_distance(
556  rt_raster rast1, int nband1,
557  rt_raster rast2, int nband2,
558  double distance,
559  int *dfwithin
560 ) {
561  LWMPOLY *surface = NULL;
562  LWGEOM *surface1 = NULL;
563  LWGEOM *surface2 = NULL;
564  double maxdist = 0;
565 
566  RASTER_DEBUG(3, "Starting");
567 
568  assert(NULL != rast1);
569  assert(NULL != rast2);
570  assert(NULL != dfwithin);
571 
572  if (nband1 < 0 && nband2 < 0) {
573  nband1 = -1;
574  nband2 = -1;
575  }
576  else {
577  assert(nband1 >= 0 && nband1 < rt_raster_get_num_bands(rast1));
578  assert(nband2 >= 0 && nband2 < rt_raster_get_num_bands(rast2));
579  }
580 
581  /* initialize to zero, false result */
582  *dfwithin = 0;
583 
584  /* same srid */
585  if (rt_raster_get_srid(rast1) != rt_raster_get_srid(rast2)) {
586  rterror("rt_raster_fully_within_distance: The two rasters provided have different SRIDs");
587  return ES_ERROR;
588  }
589 
590  /* distance cannot be less than zero */
591  if (distance < 0) {
592  rterror("rt_raster_fully_within_distance: Distance cannot be less than zero");
593  return ES_ERROR;
594  }
595 
596  /* get LWMPOLY of each band */
597  if (rt_raster_surface(rast1, nband1, &surface) != ES_NONE) {
598  rterror("rt_raster_fully_within_distance: Could not get surface of the specified band from the first raster");
599  return ES_ERROR;
600  }
601  surface1 = lwmpoly_as_lwgeom(surface);
602 
603  if (rt_raster_surface(rast2, nband2, &surface) != ES_NONE) {
604  rterror("rt_raster_fully_within_distance: Could not get surface of the specified band from the second raster");
605  lwgeom_free(surface1);
606  return ES_ERROR;
607  }
608  surface2 = lwmpoly_as_lwgeom(surface);
609 
610  /* either surface is NULL, test is false */
611  if (surface1 == NULL || surface2 == NULL) {
612  if (surface1 != NULL) lwgeom_free(surface1);
613  if (surface2 != NULL) lwgeom_free(surface2);
614  return ES_NONE;
615  }
616 
617  /* get the maximum distance between the two surfaces */
618  maxdist = lwgeom_maxdistance2d_tolerance(surface1, surface2, distance);
619 
620  lwgeom_free(surface1);
621  lwgeom_free(surface2);
622 
623  /* if distance >= maxdist, true */
624  if (FLT_EQ(maxdist, distance) || distance > maxdist)
625  *dfwithin = 1;
626 
627  RASTER_DEBUGF(3, "(maxdist, distance, dfwithin) = (%f, %f, %d)", maxdist, distance, *dfwithin);
628 
629  return ES_NONE;
630 }
631 
632 /******************************************************************************
633 * rt_raster_intersects()
634 ******************************************************************************/
635 
636 static
637 int rt_raster_intersects_algorithm(
638  rt_raster rast1, rt_raster rast2,
639  rt_band band1, rt_band band2,
640  int hasnodata1, int hasnodata2,
641  double nodata1, double nodata2
642 ) {
643  int i;
644  int byHeight = 1;
645  uint32_t dimValue;
646 
647  uint32_t row;
648  uint32_t rowoffset;
649  uint32_t col;
650  uint32_t coloffset;
651 
652  enum line_points {X1, Y1, X2, Y2};
653  enum point {pX, pY};
654  double line1[4] = {0.};
655  double line2[4] = {0.};
656  double P[2] = {0.};
657  double Qw[2] = {0.};
658  double Qr[2] = {0.};
659  double gt1[6] = {0.};
660  double gt2[6] = {0.};
661  double igt1[6] = {0};
662  double igt2[6] = {0};
663  double d;
664  double val1;
665  int noval1;
666  int isnodata1;
667  double val2;
668  int noval2;
669  int isnodata2;
670  uint32_t adjacent[8] = {0};
671 
672  double xscale;
673  double yscale;
674 
675  uint16_t width1;
676  uint16_t height1;
677  uint16_t width2;
678  uint16_t height2;
679 
680  width1 = rt_raster_get_width(rast1);
681  height1 = rt_raster_get_height(rast1);
682  width2 = rt_raster_get_width(rast2);
683  height2 = rt_raster_get_height(rast2);
684 
685  /* sampling scale */
686  xscale = fmin(rt_raster_get_x_scale(rast1), rt_raster_get_x_scale(rast2)) / 10.;
687  yscale = fmin(rt_raster_get_y_scale(rast1), rt_raster_get_y_scale(rast2)) / 10.;
688 
689  /* see if skew made rast2's rows are parallel to rast1's cols */
690  rt_raster_cell_to_geopoint(
691  rast1,
692  0, 0,
693  &(line1[X1]), &(line1[Y1]),
694  gt1
695  );
696 
697  rt_raster_cell_to_geopoint(
698  rast1,
699  0, height1,
700  &(line1[X2]), &(line1[Y2]),
701  gt1
702  );
703 
704  rt_raster_cell_to_geopoint(
705  rast2,
706  0, 0,
707  &(line2[X1]), &(line2[Y1]),
708  gt2
709  );
710 
711  rt_raster_cell_to_geopoint(
712  rast2,
713  width2, 0,
714  &(line2[X2]), &(line2[Y2]),
715  gt2
716  );
717 
718  /* Parallel lines */
719  if (FLT_EQ(((line1[X2] - line1[X1]) * (line2[Y2] - line2[Y1])),
720  ((line2[X2] - line2[X1]) * (line1[Y2] - line1[Y1]))))
721  byHeight = 0;
722 
723  if (byHeight)
724  dimValue = height2;
725  else
726  dimValue = width2;
727  RASTER_DEBUGF(4, "byHeight: %d, dimValue: %d", byHeight, dimValue);
728 
729  /* 3 x 3 search */
730  for (coloffset = 0; coloffset < 3; coloffset++) {
731  for (rowoffset = 0; rowoffset < 3; rowoffset++) {
732  /* smaller raster */
733  for (col = coloffset; col <= width1; col += 3) {
734 
735  rt_raster_cell_to_geopoint(
736  rast1,
737  col, 0,
738  &(line1[X1]), &(line1[Y1]),
739  gt1
740  );
741 
742  rt_raster_cell_to_geopoint(
743  rast1,
744  col, height1,
745  &(line1[X2]), &(line1[Y2]),
746  gt1
747  );
748 
749  /* larger raster */
750  for (row = rowoffset; row <= dimValue; row += 3) {
751 
752  if (byHeight) {
753  rt_raster_cell_to_geopoint(
754  rast2,
755  0, row,
756  &(line2[X1]), &(line2[Y1]),
757  gt2
758  );
759 
760  rt_raster_cell_to_geopoint(
761  rast2,
762  width2, row,
763  &(line2[X2]), &(line2[Y2]),
764  gt2
765  );
766  }
767  else {
768  rt_raster_cell_to_geopoint(
769  rast2,
770  row, 0,
771  &(line2[X1]), &(line2[Y1]),
772  gt2
773  );
774 
775  rt_raster_cell_to_geopoint(
776  rast2,
777  row, height2,
778  &(line2[X2]), &(line2[Y2]),
779  gt2
780  );
781  }
782 
783  RASTER_DEBUGF(4, "(col, row) = (%d, %d)", col, row);
784  RASTER_DEBUGF(4, "line1(x1, y1, x2, y2) = (%f, %f, %f, %f)",
785  line1[X1], line1[Y1], line1[X2], line1[Y2]);
786  RASTER_DEBUGF(4, "line2(x1, y1, x2, y2) = (%f, %f, %f, %f)",
787  line2[X1], line2[Y1], line2[X2], line2[Y2]);
788 
789  /* intersection */
790  /* http://en.wikipedia.org/wiki/Line-line_intersection */
791  d = ((line1[X1] - line1[X2]) * (line2[Y1] - line2[Y2])) - ((line1[Y1] - line1[Y2]) * (line2[X1] - line2[X2]));
792  /* no intersection */
793  if (FLT_EQ(d, 0.)) {
794  continue;
795  }
796 
797  P[pX] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[X1] - line2[X2])) -
798  ((line1[X1] - line1[X2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
799  P[pX] = P[pX] / d;
800 
801  P[pY] = (((line1[X1] * line1[Y2]) - (line1[Y1] * line1[X2])) * (line2[Y1] - line2[Y2])) -
802  ((line1[Y1] - line1[Y2]) * ((line2[X1] * line2[Y2]) - (line2[Y1] * line2[X2])));
803  P[pY] = P[pY] / d;
804 
805  RASTER_DEBUGF(4, "P(x, y) = (%f, %f)", P[pX], P[pY]);
806 
807  /* intersection within bounds */
808  if ((
809  (FLT_EQ(P[pX], line1[X1]) || FLT_EQ(P[pX], line1[X2])) ||
810  (P[pX] > fmin(line1[X1], line1[X2]) && (P[pX] < fmax(line1[X1], line1[X2])))
811  ) && (
812  (FLT_EQ(P[pY], line1[Y1]) || FLT_EQ(P[pY], line1[Y2])) ||
813  (P[pY] > fmin(line1[Y1], line1[Y2]) && (P[pY] < fmax(line1[Y1], line1[Y2])))
814  ) && (
815  (FLT_EQ(P[pX], line2[X1]) || FLT_EQ(P[pX], line2[X2])) ||
816  (P[pX] > fmin(line2[X1], line2[X2]) && (P[pX] < fmax(line2[X1], line2[X2])))
817  ) && (
818  (FLT_EQ(P[pY], line2[Y1]) || FLT_EQ(P[pY], line2[Y2])) ||
819  (P[pY] > fmin(line2[Y1], line2[Y2]) && (P[pY] < fmax(line2[Y1], line2[Y2])))
820  )) {
821  RASTER_DEBUG(4, "within bounds");
822 
823  for (i = 0; i < 8; i++) adjacent[i] = 0;
824 
825  /* test points around intersection */
826  for (i = 0; i < 8; i++) {
827  switch (i) {
828  case 7:
829  Qw[pX] = P[pX] - xscale;
830  Qw[pY] = P[pY] + yscale;
831  break;
832  /* 270 degrees = 09:00 */
833  case 6:
834  Qw[pX] = P[pX] - xscale;
835  Qw[pY] = P[pY];
836  break;
837  case 5:
838  Qw[pX] = P[pX] - xscale;
839  Qw[pY] = P[pY] - yscale;
840  break;
841  /* 180 degrees = 06:00 */
842  case 4:
843  Qw[pX] = P[pX];
844  Qw[pY] = P[pY] - yscale;
845  break;
846  case 3:
847  Qw[pX] = P[pX] + xscale;
848  Qw[pY] = P[pY] - yscale;
849  break;
850  /* 90 degrees = 03:00 */
851  case 2:
852  Qw[pX] = P[pX] + xscale;
853  Qw[pY] = P[pY];
854  break;
855  /* 45 degrees */
856  case 1:
857  Qw[pX] = P[pX] + xscale;
858  Qw[pY] = P[pY] + yscale;
859  break;
860  /* 0 degrees = 00:00 */
861  case 0:
862  Qw[pX] = P[pX];
863  Qw[pY] = P[pY] + yscale;
864  break;
865  }
866 
867  /* unable to convert point to cell */
868  noval1 = 0;
869  if (rt_raster_geopoint_to_cell(
870  rast1,
871  Qw[pX], Qw[pY],
872  &(Qr[pX]), &(Qr[pY]),
873  igt1
874  ) != ES_NONE) {
875  noval1 = 1;
876  }
877  /* cell is outside bounds of grid */
878  else if ((Qr[pX] < 0 || Qr[pX] >= width1) ||
879  (Qr[pY] < 0 || Qr[pY] >= height1))
880  {
881  noval1 = 1;
882  }
883  else if (hasnodata1 == FALSE)
884  val1 = 1;
885  /* unable to get value at cell */
886  else if (rt_band_get_pixel(band1, Qr[pX], Qr[pY], &val1, &isnodata1) != ES_NONE)
887  noval1 = 1;
888 
889  /* unable to convert point to cell */
890  noval2 = 0;
891  if (rt_raster_geopoint_to_cell(
892  rast2,
893  Qw[pX], Qw[pY],
894  &(Qr[pX]), &(Qr[pY]),
895  igt2
896  ) != ES_NONE) {
897  noval2 = 1;
898  }
899  /* cell is outside bounds of grid */
900  else if ((Qr[pX] < 0 || Qr[pX] >= width2) ||
901  (Qr[pY] < 0 || Qr[pY] >= height2))
902  {
903  noval2 = 1;
904  }
905  else if (hasnodata2 == FALSE)
906  val2 = 1;
907  /* unable to get value at cell */
908  else if (rt_band_get_pixel(band2, Qr[pX], Qr[pY], &val2, &isnodata2) != ES_NONE)
909  noval2 = 1;
910 
911  if (!noval1) {
912  RASTER_DEBUGF(4, "val1 = %f", val1);
913  }
914  if (!noval2) {
915  RASTER_DEBUGF(4, "val2 = %f", val2);
916  }
917 
918  /* pixels touch */
919  if (!noval1 && (
920  (hasnodata1 == FALSE) || !isnodata1
921  )) {
922  adjacent[i]++;
923  }
924  if (!noval2 && (
925  (hasnodata2 == FALSE) || !isnodata2
926  )) {
927  adjacent[i] += 3;
928  }
929 
930  /* two pixel values not present */
931  if (noval1 || noval2) {
932  RASTER_DEBUGF(4, "noval1 = %d, noval2 = %d", noval1, noval2);
933  continue;
934  }
935 
936  /* pixels valid, so intersect */
937  if (
938  ((hasnodata1 == FALSE) || !isnodata1) &&
939  ((hasnodata2 == FALSE) || !isnodata2)
940  ) {
941  RASTER_DEBUG(3, "The two rasters do intersect");
942 
943  return 1;
944  }
945  }
946 
947  /* pixels touch */
948  for (i = 0; i < 4; i++) {
949  RASTER_DEBUGF(4, "adjacent[%d] = %d, adjacent[%d] = %d"
950  , i, adjacent[i], i + 4, adjacent[i + 4]);
951  if (adjacent[i] == 0) continue;
952 
953  if (adjacent[i] + adjacent[i + 4] == 4) {
954  RASTER_DEBUG(3, "The two rasters touch");
955 
956  return 1;
957  }
958  }
959  }
960  else {
961  RASTER_DEBUG(4, "outside of bounds");
962  }
963  }
964  }
965  }
966  }
967 
968  return 0;
969 }
970 
987 rt_errorstate
988 rt_raster_intersects(
989  rt_raster rast1, int nband1,
990  rt_raster rast2, int nband2,
991  int *intersects
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 }
1295 
1296 #if POSTGIS_GEOS_VERSION >= 31400
1297 
1298 rt_raster
1299 rt_raster_intersection_fractions(
1300  const rt_raster rast_in,
1301  const LWGEOM *geom
1302 ) {
1303  rt_raster rast_out = NULL;
1304  rt_envelope rast_env;
1305  rt_band band_out = NULL;
1306  int band_num;
1307  unsigned int width, height;
1308  GEOSGeometry *geos_geom = NULL;
1309  int geos_result;
1310 
1311  /* Validate inputs */
1312  if (!rast_in || !geom)
1313  {
1314  rterror("%s: NULL input raster or geometry provided", __func__);
1315  return NULL;
1316  }
1317 
1318  /* Extract raster properties */
1319  width = rt_raster_get_width(rast_in);
1320  height = rt_raster_get_height(rast_in);
1321 
1322  if (width == 0 || height == 0)
1323  {
1324  rterror("%s: Input raster has zero width or height", __func__);
1325  return NULL;
1326  }
1327 
1328  /* Get the full bounding box of the raster in world coordinates */
1329  if (rt_raster_get_envelope(rast_in, &rast_env) != ES_NONE)
1330  {
1331  rterror("%s: Could not determine raster envelope", __func__);
1332  return NULL;
1333  }
1334 
1335  /* Shallow clone a new raster with no bands */
1336  rast_out = rt_raster_clone(rast_in, 0);
1337  /* Add a float4 band to the empty raster */
1338  band_num = rt_raster_generate_new_band(
1339  rast_out, /* rast */
1340  PT_32BF, /* pixel type */
1341  0.0, /* initial value */
1342  0, 0.0, /* hasnodata, nodataval */
1343  0); /* index */
1344 
1345  /* Get a reference to the new band */
1346  band_out = rt_raster_get_band(rast_out, band_num);
1347  if (!band_out)
1348  {
1349  rt_raster_destroy(rast_out);
1350  rterror("%s: Failed to create output band", __func__);
1351  return NULL;
1352  }
1353 
1354  /* Prepare for GEOS call */
1355  /* Initialize GEOS to handle notices and errors */
1356  initGEOS(rtinfo, lwgeom_geos_error);
1357 
1358  /* Convert the PostGIS LWGEOM into a GEOS-compatible geometry */
1359  /* Note: lwgeom_to_geos is a common function within PostGIS for this task */
1360  geos_geom = LWGEOM2GEOS(geom, 0);
1361  if (!geos_geom)
1362  {
1363  rt_raster_destroy(rast_out);
1364  rterror("%s: Could not convert LWGEOM to GEOSGeometry", __func__);
1365  return NULL;
1366  }
1367 
1368  /* Call the core GEOS function */
1369  geos_result = GEOSGridIntersectionFractions(
1370  geos_geom,
1371  rast_env.MinX, rast_env.MinY,
1372  rast_env.MaxX, rast_env.MaxY,
1373  width, height,
1374  rt_band_get_data(band_out)
1375  );
1376 
1377  /* The GEOS geometry is no longer needed after this call, so we can clean it up */
1378  GEOSGeom_destroy(geos_geom);
1379 
1380  /* Check if the GEOS function executed successfully (it returns 1 on success) */
1381  if (geos_result != 1)
1382  {
1383  rt_raster_destroy(rast_out);
1384  rterror("%s: GEOSGridIntersectionFractions call failed", __func__);
1385  return NULL;
1386  }
1387 
1388  return rast_out; /* Return the newly created raster and band*/
1389 }
1390 #endif /* POSTGIS_GEOS_VERSION >= 31400 */
1391 
FALSE
#define FALSE
Definition: dbfopen.c:72
LWGEOM2GEOS
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
Definition: liblwgeom/lwgeom_geos.c:400
lwgeom_geos_error
void lwgeom_geos_error(const char *fmt,...)
Definition: liblwgeom/lwgeom_geos.c:57
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1218
lwmpoly_free
void lwmpoly_free(LWMPOLY *mpoly)
Definition: lwmpoly.c:53
lwmpoly_as_lwgeom
LWGEOM * lwmpoly_as_lwgeom(const LWMPOLY *obj)
Definition: lwgeom.c:294
lwgeom_mindistance2d_tolerance
double lwgeom_mindistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling min distance calculations and dwithin calculations.
Definition: measures.c:210
lwgeom_maxdistance2d_tolerance
double lwgeom_maxdistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling max distance calculations and dfullywithin calculations.
Definition: measures.c:180
rterror
void rterror(const char *fmt,...) __attribute__((format(printf
Wrappers used for reporting errors and info.
FLT_NEQ
#define FLT_NEQ(x, y)
Definition: librtcore.h:2423
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:302
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:637
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:360
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:306
rt_raster_generate_new_band
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
Definition: rt_raster.c:489
rtinfo
void void rtinfo(const char *fmt,...) __attribute__((format(printf
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:825
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:686
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:1527
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:86
PT_32BF
@ PT_32BF
Definition: librtcore.h:197
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:865
rt_raster_surface
rt_errorstate rt_raster_surface(rt_raster raster, int nband, LWMPOLY **surface)
Get a raster as a surface (multipolygon).
Definition: rt_geometry.c:355
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2424
rt_raster_get_x_scale
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:154
rt_errorstate
rt_errorstate
Enum definitions.
Definition: librtcore.h:181
ES_NONE
@ ES_NONE
Definition: librtcore.h:182
ES_ERROR
@ ES_ERROR
Definition: librtcore.h:183
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:376
rt_raster_clone
rt_raster rt_raster_clone(rt_raster raster, uint8_t deep)
Clone an existing raster.
Definition: rt_raster.c:1446
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:838
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:133
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:2038
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:125
rt_band_get_data
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition: rt_band.c:551
rt_geos_spatial_test
rt_geos_spatial_test
GEOS spatial relationship tests available.
Definition: librtcore.h:219
GSR_TOUCHES
@ GSR_TOUCHES
Definition: librtcore.h:221
GSR_COVERS
@ GSR_COVERS
Definition: librtcore.h:224
GSR_COVEREDBY
@ GSR_COVEREDBY
Definition: librtcore.h:225
GSR_CONTAINSPROPERLY
@ GSR_CONTAINSPROPERLY
Definition: librtcore.h:223
GSR_OVERLAPS
@ GSR_OVERLAPS
Definition: librtcore.h:220
GSR_CONTAINS
@ GSR_CONTAINS
Definition: librtcore.h:222
rt_raster_get_y_scale
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:163
rt_raster_get_envelope
rt_errorstate rt_raster_get_envelope(rt_raster raster, rt_envelope *env)
Get raster's envelope.
Definition: rt_raster.c:782
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:385
librtcore.h
This library is the generic raster handling section of PostGIS.
librtcore_internal.h
contains
Datum contains(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:496
coveredby
Datum coveredby(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:806
covers
Datum covers(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:722
touches
Datum touches(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:229
overlaps
Datum overlaps(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:361
within
Datum within(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:577
distance
static double distance(double x1, double y1, double x2, double y2)
Definition: lwtree.c:1032
rt_raster_contains
rt_errorstate rt_raster_contains(rt_raster rast1, int nband1, rt_raster rast2, int nband2, int *contains)
Return ES_ERROR if error occurred in function.
Definition: rt_spatial_relationship.c:334
rt_raster_fully_within_distance
rt_errorstate rt_raster_fully_within_distance(rt_raster rast1, int nband1, rt_raster rast2, int nband2, double distance, int *dfwithin)
Return ES_ERROR if error occurred in function.
Definition: rt_spatial_relationship.c:555
rt_raster_contains_properly
rt_errorstate rt_raster_contains_properly(rt_raster rast1, int nband1, rt_raster rast2, int nband2, int *contains)
Return ES_ERROR if error occurred in function.
Definition: rt_spatial_relationship.c:365
rt_raster_intersects
rt_errorstate rt_raster_intersects(rt_raster rast1, int nband1, rt_raster rast2, int nband2, int *intersects)
Return zero if error occurred in function.
Definition: rt_spatial_relationship.c:988
rt_raster_coveredby
rt_errorstate rt_raster_coveredby(rt_raster rast1, int nband1, rt_raster rast2, int nband2, int *coveredby)
Return ES_ERROR if error occurred in function.
Definition: rt_spatial_relationship.c:427
rt_raster_touches
rt_errorstate rt_raster_touches(rt_raster rast1, int nband1, rt_raster rast2, int nband2, int *touches)
Return ES_ERROR if error occurred in function.
Definition: rt_spatial_relationship.c:303
rt_raster_geos_spatial_relationship
static rt_errorstate rt_raster_geos_spatial_relationship(rt_raster rast1, int nband1, rt_raster rast2, int nband2, rt_geos_spatial_test testtype, int *testresult)
Definition: rt_spatial_relationship.c:137
rt_raster_overlaps
rt_errorstate rt_raster_overlaps(rt_raster rast1, int nband1, rt_raster rast2, int nband2, int *overlaps)
Return ES_ERROR if error occurred in function.
Definition: rt_spatial_relationship.c:272
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_raster_within_distance
rt_errorstate rt_raster_within_distance(rt_raster rast1, int nband1, rt_raster rast2, int nband2, double distance, int *dwithin)
Return ES_ERROR if error occurred in function.
Definition: rt_spatial_relationship.c:460
rt_raster_same_alignment
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
Definition: rt_spatial_relationship.c:45
rt_raster_covers
rt_errorstate rt_raster_covers(rt_raster rast1, int nband1, rt_raster rast2, int nband2, int *covers)
Return ES_ERROR if error occurred in function.
Definition: rt_spatial_relationship.c:396
rt_band_t
Definition: librtcore.h:2503
rt_envelope::MinX
double MinX
Definition: librtcore.h:167
rt_envelope::MaxX
double MaxX
Definition: librtcore.h:168
rt_envelope::MinY
double MinY
Definition: librtcore.h:169
rt_envelope::MaxY
double MaxY
Definition: librtcore.h:170
rt_envelope
Definition: librtcore.h:166
rt_raster_t::scaleX
double scaleX
Definition: librtcore.h:2483
rt_raster_t::skewY
double skewY
Definition: librtcore.h:2488
rt_raster_t::ipX
double ipX
Definition: librtcore.h:2485
rt_raster_t::scaleY
double scaleY
Definition: librtcore.h:2484
rt_raster_t::ipY
double ipY
Definition: librtcore.h:2486
rt_raster_t::skewX
double skewX
Definition: librtcore.h:2487
rt_raster_t
Definition: librtcore.h:2474