PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches
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 */
44rt_errorstate
45rt_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
136static
137rt_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
272rt_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
303rt_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
334rt_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
365rt_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
396rt_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
427rt_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
460rt_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
555rt_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
636static
637int 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
987rt_errorstate
988rt_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
1298rt_raster
1299rt_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:1246
lwmpoly_free
void lwmpoly_free(LWMPOLY *mpoly)
Definition lwmpoly.c:53
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:222
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:192
lwmpoly_as_lwgeom
LWGEOM * lwmpoly_as_lwgeom(const LWMPOLY *obj)
Definition lwgeom.c:322
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:2435
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition librtcore.h:304
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:308
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:833
rt_band_get_data
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition rt_band.c:559
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:1551
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:199
rt_band_get_isnodata_flag
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition rt_band.c:873
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:356
FLT_EQ
#define FLT_EQ(x, y)
Definition librtcore.h:2436
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:839
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:2067
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition rt_raster.c:125
rt_geos_spatial_test
rt_geos_spatial_test
GEOS spatial relationship tests available.
Definition librtcore.h:221
GSR_TOUCHES
@ GSR_TOUCHES
Definition librtcore.h:223
GSR_COVERS
@ GSR_COVERS
Definition librtcore.h:226
GSR_COVEREDBY
@ GSR_COVEREDBY
Definition librtcore.h:227
GSR_CONTAINSPROPERLY
@ GSR_CONTAINSPROPERLY
Definition librtcore.h:225
GSR_OVERLAPS
@ GSR_OVERLAPS
Definition librtcore.h:222
GSR_CONTAINS
@ GSR_CONTAINS
Definition librtcore.h:224
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:2515
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:2495
rt_raster_t::skewY
double skewY
Definition librtcore.h:2500
rt_raster_t::ipX
double ipX
Definition librtcore.h:2497
rt_raster_t::scaleY
double scaleY
Definition librtcore.h:2496
rt_raster_t::ipY
double ipY
Definition librtcore.h:2498
rt_raster_t::skewX
double skewX
Definition librtcore.h:2499
rt_raster_t
Definition librtcore.h:2486