PostGIS  3.7.0dev-r@@SVN_REVISION@@
rtpg_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 <postgres.h> /* for palloc */
31 #include <fmgr.h>
32 #include <utils/builtins.h>
33 
34 #include "../../postgis_config.h"
35 
36 #include "lwgeom_pg.h"
37 #include "rtpostgis.h"
38 
39 /* determine if two rasters intersect */
40 Datum RASTER_intersects(PG_FUNCTION_ARGS);
41 
42 /* determine if two rasters overlap */
43 Datum RASTER_overlaps(PG_FUNCTION_ARGS);
44 
45 /* determine if two rasters touch */
46 Datum RASTER_touches(PG_FUNCTION_ARGS);
47 
48 /* determine if the first raster contains the second raster */
49 Datum RASTER_contains(PG_FUNCTION_ARGS);
50 
51 /* determine if the first raster contains properly the second raster */
52 Datum RASTER_containsProperly(PG_FUNCTION_ARGS);
53 
54 /* determine if the first raster covers the second raster */
55 Datum RASTER_covers(PG_FUNCTION_ARGS);
56 
57 /* determine if the first raster is covered by the second raster */
58 Datum RASTER_coveredby(PG_FUNCTION_ARGS);
59 
60 /* determine if the two rasters are within the specified distance of each other */
61 Datum RASTER_dwithin(PG_FUNCTION_ARGS);
62 
63 /* determine if the two rasters are fully within the specified distance of each other */
64 Datum RASTER_dfullywithin(PG_FUNCTION_ARGS);
65 
66 /* determine if two rasters are aligned */
67 Datum RASTER_sameAlignment(PG_FUNCTION_ARGS);
68 Datum RASTER_notSameAlignmentReason(PG_FUNCTION_ARGS);
69 
70 /* calculate the proportion of each cell in a grid covered */
71 /* by an input geometry */
72 Datum RASTER_intersectionFractions(PG_FUNCTION_ARGS);
73 
77 PG_FUNCTION_INFO_V1(RASTER_intersects);
78 Datum RASTER_intersects(PG_FUNCTION_ARGS)
79 {
80  const uint32_t set_count = 2;
81  rt_pgraster *pgrast[2];
82  int pgrastpos[2] = {-1, -1};
83  rt_raster rast[2] = {NULL};
84  uint32_t bandindex[2] = {0};
85  uint32_t hasbandindex[2] = {0};
86 
87  uint32_t i;
88  uint32_t j;
89  uint32_t k;
90  uint32_t numBands;
91  int rtn;
92  int result;
93 
94  for (i = 0, j = 0; i < set_count; i++) {
95  /* pgrast is null, return null */
96  if (PG_ARGISNULL(j)) {
97  for (k = 0; k < i; k++) {
98  rt_raster_destroy(rast[k]);
99  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
100  }
101  PG_RETURN_NULL();
102  }
103  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
104  pgrastpos[i] = j;
105  j++;
106 
107  /* raster */
108  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
109  if (!rast[i]) {
110  for (k = 0; k <= i; k++) {
111  if (k < i)
112  rt_raster_destroy(rast[k]);
113  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
114  }
115  elog(ERROR, "RASTER_intersects: Could not deserialize the %s raster", i < 1 ? "first" : "second");
116  PG_RETURN_NULL();
117  }
118 
119  /* numbands */
120  numBands = rt_raster_get_num_bands(rast[i]);
121  if (numBands < 1) {
122  elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
123  if (i > 0) i++;
124  for (k = 0; k < i; k++) {
125  rt_raster_destroy(rast[k]);
126  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
127  }
128  PG_RETURN_NULL();
129  }
130 
131  /* band index */
132  if (!PG_ARGISNULL(j)) {
133  bandindex[i] = PG_GETARG_INT32(j);
134  if (bandindex[i] < 1 || bandindex[i] > numBands) {
135  elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
136  if (i > 0) i++;
137  for (k = 0; k < i; k++) {
138  rt_raster_destroy(rast[k]);
139  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
140  }
141  PG_RETURN_NULL();
142  }
143  hasbandindex[i] = 1;
144  }
145  else
146  hasbandindex[i] = 0;
147  POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
148  POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
149  j++;
150  }
151 
152  /* hasbandindex must be balanced */
153  if (
154  (hasbandindex[0] && !hasbandindex[1]) ||
155  (!hasbandindex[0] && hasbandindex[1])
156  ) {
157  elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
158  for (k = 0; k < set_count; k++) {
159  rt_raster_destroy(rast[k]);
160  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
161  }
162  PG_RETURN_NULL();
163  }
164 
165  /* SRID must match */
166  if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
167  for (k = 0; k < set_count; k++) {
168  rt_raster_destroy(rast[k]);
169  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
170  }
171  elog(ERROR, "The two rasters provided have different SRIDs");
172  PG_RETURN_NULL();
173  }
174 
175  rtn = rt_raster_intersects(
176  rast[0], (hasbandindex[0] ? (int)bandindex[0] - 1 : -1),
177  rast[1], (hasbandindex[1] ? (int)bandindex[1] - 1 : -1),
178  &result
179  );
180  for (k = 0; k < set_count; k++) {
181  rt_raster_destroy(rast[k]);
182  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
183  }
184 
185  if (rtn != ES_NONE) {
186  elog(ERROR, "RASTER_intersects: Could not test for intersection on the two rasters");
187  PG_RETURN_NULL();
188  }
189 
190  PG_RETURN_BOOL(result);
191 }
192 
196 PG_FUNCTION_INFO_V1(RASTER_overlaps);
197 Datum RASTER_overlaps(PG_FUNCTION_ARGS)
198 {
199  const uint32_t set_count = 2;
200  rt_pgraster *pgrast[2];
201  int pgrastpos[2] = {-1, -1};
202  rt_raster rast[2] = {NULL};
203  uint32_t bandindex[2] = {0};
204  uint32_t hasbandindex[2] = {0};
205 
206  uint32_t i;
207  uint32_t j;
208  uint32_t k;
209  uint32_t numBands;
210  int rtn;
211  int result;
212 
213  for (i = 0, j = 0; i < set_count; i++) {
214  /* pgrast is null, return null */
215  if (PG_ARGISNULL(j)) {
216  for (k = 0; k < i; k++) {
217  rt_raster_destroy(rast[k]);
218  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
219  }
220  PG_RETURN_NULL();
221  }
222  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
223  pgrastpos[i] = j;
224  j++;
225 
226  /* raster */
227  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
228  if (!rast[i]) {
229  for (k = 0; k <= i; k++) {
230  if (k < i)
231  rt_raster_destroy(rast[k]);
232  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
233  }
234  elog(ERROR, "RASTER_overlaps: Could not deserialize the %s raster", i < 1 ? "first" : "second");
235  PG_RETURN_NULL();
236  }
237 
238  /* numbands */
239  numBands = rt_raster_get_num_bands(rast[i]);
240  if (numBands < 1) {
241  elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
242  if (i > 0) i++;
243  for (k = 0; k < i; k++) {
244  rt_raster_destroy(rast[k]);
245  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
246  }
247  PG_RETURN_NULL();
248  }
249 
250  /* band index */
251  if (!PG_ARGISNULL(j)) {
252  bandindex[i] = PG_GETARG_INT32(j);
253  if (bandindex[i] < 1 || bandindex[i] > numBands) {
254  elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
255  if (i > 0) i++;
256  for (k = 0; k < i; k++) {
257  rt_raster_destroy(rast[k]);
258  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
259  }
260  PG_RETURN_NULL();
261  }
262  hasbandindex[i] = 1;
263  }
264  else
265  hasbandindex[i] = 0;
266  POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
267  POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
268  j++;
269  }
270 
271  /* hasbandindex must be balanced */
272  if (
273  (hasbandindex[0] && !hasbandindex[1]) ||
274  (!hasbandindex[0] && hasbandindex[1])
275  ) {
276  elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
277  for (k = 0; k < set_count; k++) {
278  rt_raster_destroy(rast[k]);
279  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
280  }
281  PG_RETURN_NULL();
282  }
283 
284  /* SRID must match */
285  if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
286  for (k = 0; k < set_count; k++) {
287  rt_raster_destroy(rast[k]);
288  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
289  }
290  elog(ERROR, "The two rasters provided have different SRIDs");
291  PG_RETURN_NULL();
292  }
293 
294  rtn = rt_raster_overlaps(
295  rast[0], (hasbandindex[0] ? (int)bandindex[0] - 1 : -1),
296  rast[1], (hasbandindex[1] ? (int)bandindex[1] - 1 : -1),
297  &result
298  );
299  for (k = 0; k < set_count; k++) {
300  rt_raster_destroy(rast[k]);
301  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
302  }
303 
304  if (rtn != ES_NONE) {
305  elog(ERROR, "RASTER_overlaps: Could not test for overlap on the two rasters");
306  PG_RETURN_NULL();
307  }
308 
309  PG_RETURN_BOOL(result);
310 }
311 
315 PG_FUNCTION_INFO_V1(RASTER_touches);
316 Datum RASTER_touches(PG_FUNCTION_ARGS)
317 {
318  const uint32_t set_count = 2;
319  rt_pgraster *pgrast[2];
320  int pgrastpos[2] = {-1, -1};
321  rt_raster rast[2] = {NULL};
322  uint32_t bandindex[2] = {0};
323  uint32_t hasbandindex[2] = {0};
324 
325  uint32_t i;
326  uint32_t j;
327  uint32_t k;
328  uint32_t numBands;
329  int rtn;
330  int result;
331 
332  for (i = 0, j = 0; i < set_count; i++) {
333  /* pgrast is null, return null */
334  if (PG_ARGISNULL(j)) {
335  for (k = 0; k < i; k++) {
336  rt_raster_destroy(rast[k]);
337  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
338  }
339  PG_RETURN_NULL();
340  }
341  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
342  pgrastpos[i] = j;
343  j++;
344 
345  /* raster */
346  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
347  if (!rast[i]) {
348  for (k = 0; k <= i; k++) {
349  if (k < i)
350  rt_raster_destroy(rast[k]);
351  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
352  }
353  elog(ERROR, "RASTER_touches: Could not deserialize the %s raster", i < 1 ? "first" : "second");
354  PG_RETURN_NULL();
355  }
356 
357  /* numbands */
358  numBands = rt_raster_get_num_bands(rast[i]);
359  if (numBands < 1) {
360  elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
361  if (i > 0) i++;
362  for (k = 0; k < i; k++) {
363  rt_raster_destroy(rast[k]);
364  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
365  }
366  PG_RETURN_NULL();
367  }
368 
369  /* band index */
370  if (!PG_ARGISNULL(j)) {
371  bandindex[i] = PG_GETARG_INT32(j);
372  if (bandindex[i] < 1 || bandindex[i] > numBands) {
373  elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
374  if (i > 0) i++;
375  for (k = 0; k < i; k++) {
376  rt_raster_destroy(rast[k]);
377  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
378  }
379  PG_RETURN_NULL();
380  }
381  hasbandindex[i] = 1;
382  }
383  else
384  hasbandindex[i] = 0;
385  POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
386  POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
387  j++;
388  }
389 
390  /* hasbandindex must be balanced */
391  if (
392  (hasbandindex[0] && !hasbandindex[1]) ||
393  (!hasbandindex[0] && hasbandindex[1])
394  ) {
395  elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
396  for (k = 0; k < set_count; k++) {
397  rt_raster_destroy(rast[k]);
398  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
399  }
400  PG_RETURN_NULL();
401  }
402 
403  /* SRID must match */
404  if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
405  for (k = 0; k < set_count; k++) {
406  rt_raster_destroy(rast[k]);
407  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
408  }
409  elog(ERROR, "The two rasters provided have different SRIDs");
410  PG_RETURN_NULL();
411  }
412 
413  rtn = rt_raster_touches(
414  rast[0], (hasbandindex[0] ? (int)bandindex[0] - 1 : -1),
415  rast[1], (hasbandindex[1] ? (int)bandindex[1] - 1 : -1),
416  &result
417  );
418  for (k = 0; k < set_count; k++) {
419  rt_raster_destroy(rast[k]);
420  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
421  }
422 
423  if (rtn != ES_NONE) {
424  elog(ERROR, "RASTER_touches: Could not test for touch on the two rasters");
425  PG_RETURN_NULL();
426  }
427 
428  PG_RETURN_BOOL(result);
429 }
430 
434 PG_FUNCTION_INFO_V1(RASTER_contains);
435 Datum RASTER_contains(PG_FUNCTION_ARGS)
436 {
437  const uint32_t set_count = 2;
438  rt_pgraster *pgrast[2];
439  int pgrastpos[2] = {-1, -1};
440  rt_raster rast[2] = {NULL};
441  uint32_t bandindex[2] = {0};
442  uint32_t hasbandindex[2] = {0};
443 
444  uint32_t i;
445  uint32_t j;
446  uint32_t k;
447  uint32_t numBands;
448  int rtn;
449  int result;
450 
451  for (i = 0, j = 0; i < set_count; i++) {
452  /* pgrast is null, return null */
453  if (PG_ARGISNULL(j)) {
454  for (k = 0; k < i; k++) {
455  rt_raster_destroy(rast[k]);
456  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
457  }
458  PG_RETURN_NULL();
459  }
460  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
461  pgrastpos[i] = j;
462  j++;
463 
464  /* raster */
465  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
466  if (!rast[i]) {
467  for (k = 0; k <= i; k++) {
468  if (k < i)
469  rt_raster_destroy(rast[k]);
470  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
471  }
472  elog(ERROR, "RASTER_contains: Could not deserialize the %s raster", i < 1 ? "first" : "second");
473  PG_RETURN_NULL();
474  }
475 
476  /* numbands */
477  numBands = rt_raster_get_num_bands(rast[i]);
478  if (numBands < 1) {
479  elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
480  if (i > 0) i++;
481  for (k = 0; k < i; k++) {
482  rt_raster_destroy(rast[k]);
483  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
484  }
485  PG_RETURN_NULL();
486  }
487 
488  /* band index */
489  if (!PG_ARGISNULL(j)) {
490  bandindex[i] = PG_GETARG_INT32(j);
491  if (bandindex[i] < 1 || bandindex[i] > numBands) {
492  elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
493  if (i > 0) i++;
494  for (k = 0; k < i; k++) {
495  rt_raster_destroy(rast[k]);
496  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
497  }
498  PG_RETURN_NULL();
499  }
500  hasbandindex[i] = 1;
501  }
502  else
503  hasbandindex[i] = 0;
504  POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
505  POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
506  j++;
507  }
508 
509  /* hasbandindex must be balanced */
510  if (
511  (hasbandindex[0] && !hasbandindex[1]) ||
512  (!hasbandindex[0] && hasbandindex[1])
513  ) {
514  elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
515  for (k = 0; k < set_count; k++) {
516  rt_raster_destroy(rast[k]);
517  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
518  }
519  PG_RETURN_NULL();
520  }
521 
522  /* SRID must match */
523  if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
524  for (k = 0; k < set_count; k++) {
525  rt_raster_destroy(rast[k]);
526  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
527  }
528  elog(ERROR, "The two rasters provided have different SRIDs");
529  PG_RETURN_NULL();
530  }
531 
532  rtn = rt_raster_contains(
533  rast[0], (hasbandindex[0] ? (int)bandindex[0] - 1 : -1),
534  rast[1], (hasbandindex[1] ? (int)bandindex[1] - 1 : -1),
535  &result
536  );
537  for (k = 0; k < set_count; k++) {
538  rt_raster_destroy(rast[k]);
539  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
540  }
541 
542  if (rtn != ES_NONE) {
543  elog(ERROR, "RASTER_contains: Could not test that the first raster contains the second raster");
544  PG_RETURN_NULL();
545  }
546 
547  PG_RETURN_BOOL(result);
548 }
549 
553 PG_FUNCTION_INFO_V1(RASTER_containsProperly);
554 Datum RASTER_containsProperly(PG_FUNCTION_ARGS)
555 {
556  const uint32_t set_count = 2;
557  rt_pgraster *pgrast[2];
558  int pgrastpos[2] = {-1, -1};
559  rt_raster rast[2] = {NULL};
560  uint32_t bandindex[2] = {0};
561  uint32_t hasbandindex[2] = {0};
562 
563  uint32_t i;
564  uint32_t j;
565  uint32_t k;
566  uint32_t numBands;
567  int rtn;
568  int result;
569 
570  for (i = 0, j = 0; i < set_count; i++) {
571  /* pgrast is null, return null */
572  if (PG_ARGISNULL(j)) {
573  for (k = 0; k < i; k++) {
574  rt_raster_destroy(rast[k]);
575  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
576  }
577  PG_RETURN_NULL();
578  }
579  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
580  pgrastpos[i] = j;
581  j++;
582 
583  /* raster */
584  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
585  if (!rast[i]) {
586  for (k = 0; k <= i; k++) {
587  if (k < i)
588  rt_raster_destroy(rast[k]);
589  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
590  }
591  elog(ERROR, "RASTER_containsProperly: Could not deserialize the %s raster", i < 1 ? "first" : "second");
592  PG_RETURN_NULL();
593  }
594 
595  /* numbands */
596  numBands = rt_raster_get_num_bands(rast[i]);
597  if (numBands < 1) {
598  elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
599  if (i > 0) i++;
600  for (k = 0; k < i; k++) {
601  rt_raster_destroy(rast[k]);
602  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
603  }
604  PG_RETURN_NULL();
605  }
606 
607  /* band index */
608  if (!PG_ARGISNULL(j)) {
609  bandindex[i] = PG_GETARG_INT32(j);
610  if (bandindex[i] < 1 || bandindex[i] > numBands) {
611  elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
612  if (i > 0) i++;
613  for (k = 0; k < i; k++) {
614  rt_raster_destroy(rast[k]);
615  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
616  }
617  PG_RETURN_NULL();
618  }
619  hasbandindex[i] = 1;
620  }
621  else
622  hasbandindex[i] = 0;
623  POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
624  POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
625  j++;
626  }
627 
628  /* hasbandindex must be balanced */
629  if (
630  (hasbandindex[0] && !hasbandindex[1]) ||
631  (!hasbandindex[0] && hasbandindex[1])
632  ) {
633  elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
634  for (k = 0; k < set_count; k++) {
635  rt_raster_destroy(rast[k]);
636  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
637  }
638  PG_RETURN_NULL();
639  }
640 
641  /* SRID must match */
642  if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
643  for (k = 0; k < set_count; k++) {
644  rt_raster_destroy(rast[k]);
645  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
646  }
647  elog(ERROR, "The two rasters provided have different SRIDs");
648  PG_RETURN_NULL();
649  }
650 
651  rtn = rt_raster_contains_properly(
652  rast[0], (hasbandindex[0] ? (int)bandindex[0] - 1 : -1),
653  rast[1], (hasbandindex[1] ? (int)bandindex[1] - 1 : -1),
654  &result
655  );
656  for (k = 0; k < set_count; k++) {
657  rt_raster_destroy(rast[k]);
658  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
659  }
660 
661  if (rtn != ES_NONE) {
662  elog(ERROR, "RASTER_containsProperly: Could not test that the first raster contains properly the second raster");
663  PG_RETURN_NULL();
664  }
665 
666  PG_RETURN_BOOL(result);
667 }
668 
672 PG_FUNCTION_INFO_V1(RASTER_covers);
673 Datum RASTER_covers(PG_FUNCTION_ARGS)
674 {
675  const uint32_t set_count = 2;
676  rt_pgraster *pgrast[2];
677  int pgrastpos[2] = {-1, -1};
678  rt_raster rast[2] = {NULL};
679  uint32_t bandindex[2] = {0};
680  uint32_t hasbandindex[2] = {0};
681 
682  uint32_t i;
683  uint32_t j;
684  uint32_t k;
685  uint32_t numBands;
686  int rtn;
687  int result;
688 
689  for (i = 0, j = 0; i < set_count; i++) {
690  /* pgrast is null, return null */
691  if (PG_ARGISNULL(j)) {
692  for (k = 0; k < i; k++) {
693  rt_raster_destroy(rast[k]);
694  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
695  }
696  PG_RETURN_NULL();
697  }
698  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
699  pgrastpos[i] = j;
700  j++;
701 
702  /* raster */
703  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
704  if (!rast[i]) {
705  for (k = 0; k <= i; k++) {
706  if (k < i)
707  rt_raster_destroy(rast[k]);
708  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
709  }
710  elog(ERROR, "RASTER_covers: Could not deserialize the %s raster", i < 1 ? "first" : "second");
711  PG_RETURN_NULL();
712  }
713 
714  /* numbands */
715  numBands = rt_raster_get_num_bands(rast[i]);
716  if (numBands < 1) {
717  elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
718  if (i > 0) i++;
719  for (k = 0; k < i; k++) {
720  rt_raster_destroy(rast[k]);
721  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
722  }
723  PG_RETURN_NULL();
724  }
725 
726  /* band index */
727  if (!PG_ARGISNULL(j)) {
728  bandindex[i] = PG_GETARG_INT32(j);
729  if (bandindex[i] < 1 || bandindex[i] > numBands) {
730  elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
731  if (i > 0) i++;
732  for (k = 0; k < i; k++) {
733  rt_raster_destroy(rast[k]);
734  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
735  }
736  PG_RETURN_NULL();
737  }
738  hasbandindex[i] = 1;
739  }
740  else
741  hasbandindex[i] = 0;
742  POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
743  POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
744  j++;
745  }
746 
747  /* hasbandindex must be balanced */
748  if (
749  (hasbandindex[0] && !hasbandindex[1]) ||
750  (!hasbandindex[0] && hasbandindex[1])
751  ) {
752  elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
753  for (k = 0; k < set_count; k++) {
754  rt_raster_destroy(rast[k]);
755  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
756  }
757  PG_RETURN_NULL();
758  }
759 
760  /* SRID must match */
761  if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
762  for (k = 0; k < set_count; k++) {
763  rt_raster_destroy(rast[k]);
764  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
765  }
766  elog(ERROR, "The two rasters provided have different SRIDs");
767  PG_RETURN_NULL();
768  }
769 
770  rtn = rt_raster_covers(
771  rast[0], (hasbandindex[0] ? (int)bandindex[0] - 1 : -1),
772  rast[1], (hasbandindex[1] ? (int)bandindex[1] - 1 : -1),
773  &result
774  );
775  for (k = 0; k < set_count; k++) {
776  rt_raster_destroy(rast[k]);
777  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
778  }
779 
780  if (rtn != ES_NONE) {
781  elog(ERROR, "RASTER_covers: Could not test that the first raster covers the second raster");
782  PG_RETURN_NULL();
783  }
784 
785  PG_RETURN_BOOL(result);
786 }
787 
791 PG_FUNCTION_INFO_V1(RASTER_coveredby);
792 Datum RASTER_coveredby(PG_FUNCTION_ARGS)
793 {
794  const uint32_t set_count = 2;
795  rt_pgraster *pgrast[2];
796  int pgrastpos[2] = {-1, -1};
797  rt_raster rast[2] = {NULL};
798  uint32_t bandindex[2] = {0};
799  uint32_t hasbandindex[2] = {0};
800 
801  uint32_t i;
802  uint32_t j;
803  uint32_t k;
804  uint32_t numBands;
805  int rtn;
806  int result;
807 
808  for (i = 0, j = 0; i < set_count; i++) {
809  /* pgrast is null, return null */
810  if (PG_ARGISNULL(j)) {
811  for (k = 0; k < i; k++) {
812  rt_raster_destroy(rast[k]);
813  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
814  }
815  PG_RETURN_NULL();
816  }
817  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
818  pgrastpos[i] = j;
819  j++;
820 
821  /* raster */
822  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
823  if (!rast[i]) {
824  for (k = 0; k <= i; k++) {
825  if (k < i)
826  rt_raster_destroy(rast[k]);
827  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
828  }
829  elog(ERROR, "RASTER_coveredby: Could not deserialize the %s raster", i < 1 ? "first" : "second");
830  PG_RETURN_NULL();
831  }
832 
833  /* numbands */
834  numBands = rt_raster_get_num_bands(rast[i]);
835  if (numBands < 1) {
836  elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
837  if (i > 0) i++;
838  for (k = 0; k < i; k++) {
839  rt_raster_destroy(rast[k]);
840  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
841  }
842  PG_RETURN_NULL();
843  }
844 
845  /* band index */
846  if (!PG_ARGISNULL(j)) {
847  bandindex[i] = PG_GETARG_INT32(j);
848  if (bandindex[i] < 1 || bandindex[i] > numBands) {
849  elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
850  if (i > 0) i++;
851  for (k = 0; k < i; k++) {
852  rt_raster_destroy(rast[k]);
853  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
854  }
855  PG_RETURN_NULL();
856  }
857  hasbandindex[i] = 1;
858  }
859  else
860  hasbandindex[i] = 0;
861  POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
862  POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
863  j++;
864  }
865 
866  /* hasbandindex must be balanced */
867  if (
868  (hasbandindex[0] && !hasbandindex[1]) ||
869  (!hasbandindex[0] && hasbandindex[1])
870  ) {
871  elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
872  for (k = 0; k < set_count; k++) {
873  rt_raster_destroy(rast[k]);
874  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
875  }
876  PG_RETURN_NULL();
877  }
878 
879  /* SRID must match */
880  if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
881  for (k = 0; k < set_count; k++) {
882  rt_raster_destroy(rast[k]);
883  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
884  }
885  elog(ERROR, "The two rasters provided have different SRIDs");
886  PG_RETURN_NULL();
887  }
888 
889  rtn = rt_raster_coveredby(
890  rast[0], (hasbandindex[0] ? (int)bandindex[0] - 1 : -1),
891  rast[1], (hasbandindex[1] ? (int)bandindex[1] - 1 : -1),
892  &result
893  );
894  for (k = 0; k < set_count; k++) {
895  rt_raster_destroy(rast[k]);
896  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
897  }
898 
899  if (rtn != ES_NONE) {
900  elog(ERROR, "RASTER_coveredby: Could not test that the first raster is covered by the second raster");
901  PG_RETURN_NULL();
902  }
903 
904  PG_RETURN_BOOL(result);
905 }
906 
910 PG_FUNCTION_INFO_V1(RASTER_dwithin);
911 Datum RASTER_dwithin(PG_FUNCTION_ARGS)
912 {
913  const uint32_t set_count = 2;
914  rt_pgraster *pgrast[2];
915  int pgrastpos[2] = {-1, -1};
916  rt_raster rast[2] = {NULL};
917  uint32_t bandindex[2] = {0};
918  uint32_t hasbandindex[2] = {0};
919  double distance = 0;
920 
921  uint32_t i;
922  uint32_t j;
923  uint32_t k;
924  uint32_t numBands;
925  int rtn;
926  int result;
927 
928  for (i = 0, j = 0; i < set_count; i++) {
929  /* pgrast is null, return null */
930  if (PG_ARGISNULL(j)) {
931  for (k = 0; k < i; k++) {
932  rt_raster_destroy(rast[k]);
933  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
934  }
935  PG_RETURN_NULL();
936  }
937  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
938  pgrastpos[i] = j;
939  j++;
940 
941  /* raster */
942  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
943  if (!rast[i]) {
944  for (k = 0; k <= i; k++) {
945  if (k < i)
946  rt_raster_destroy(rast[k]);
947  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
948  }
949  elog(ERROR, "RASTER_dwithin: Could not deserialize the %s raster", i < 1 ? "first" : "second");
950  PG_RETURN_NULL();
951  }
952 
953  /* numbands */
954  numBands = rt_raster_get_num_bands(rast[i]);
955  if (numBands < 1) {
956  elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
957  if (i > 0) i++;
958  for (k = 0; k < i; k++) {
959  rt_raster_destroy(rast[k]);
960  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
961  }
962  PG_RETURN_NULL();
963  }
964 
965  /* band index */
966  if (!PG_ARGISNULL(j)) {
967  bandindex[i] = PG_GETARG_INT32(j);
968  if (bandindex[i] < 1 || bandindex[i] > numBands) {
969  elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
970  if (i > 0) i++;
971  for (k = 0; k < i; k++) {
972  rt_raster_destroy(rast[k]);
973  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
974  }
975  PG_RETURN_NULL();
976  }
977  hasbandindex[i] = 1;
978  }
979  else
980  hasbandindex[i] = 0;
981  POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
982  POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
983  j++;
984  }
985 
986  /* distance */
987  if (PG_ARGISNULL(4)) {
988  elog(NOTICE, "Distance cannot be NULL. Returning NULL");
989  for (k = 0; k < set_count; k++) {
990  rt_raster_destroy(rast[k]);
991  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
992  }
993  PG_RETURN_NULL();
994  }
995 
996  distance = PG_GETARG_FLOAT8(4);
997  if (distance < 0) {
998  elog(NOTICE, "Distance cannot be less than zero. Returning NULL");
999  for (k = 0; k < set_count; k++) {
1000  rt_raster_destroy(rast[k]);
1001  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1002  }
1003  PG_RETURN_NULL();
1004  }
1005 
1006  /* hasbandindex must be balanced */
1007  if (
1008  (hasbandindex[0] && !hasbandindex[1]) ||
1009  (!hasbandindex[0] && hasbandindex[1])
1010  ) {
1011  elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
1012  for (k = 0; k < set_count; k++) {
1013  rt_raster_destroy(rast[k]);
1014  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1015  }
1016  PG_RETURN_NULL();
1017  }
1018 
1019  /* SRID must match */
1020  if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
1021  for (k = 0; k < set_count; k++) {
1022  rt_raster_destroy(rast[k]);
1023  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1024  }
1025  elog(ERROR, "The two rasters provided have different SRIDs");
1026  PG_RETURN_NULL();
1027  }
1028 
1029  rtn = rt_raster_within_distance(
1030  rast[0], (hasbandindex[0] ? (int)bandindex[0] - 1 : -1),
1031  rast[1], (hasbandindex[1] ? (int)bandindex[1] - 1 : -1),
1032  distance,
1033  &result
1034  );
1035  for (k = 0; k < set_count; k++) {
1036  rt_raster_destroy(rast[k]);
1037  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1038  }
1039 
1040  if (rtn != ES_NONE) {
1041  elog(ERROR, "RASTER_dwithin: Could not test that the two rasters are within the specified distance of each other");
1042  PG_RETURN_NULL();
1043  }
1044 
1045  PG_RETURN_BOOL(result);
1046 }
1047 
1051 PG_FUNCTION_INFO_V1(RASTER_dfullywithin);
1052 Datum RASTER_dfullywithin(PG_FUNCTION_ARGS)
1053 {
1054  const uint32_t set_count = 2;
1055  rt_pgraster *pgrast[2];
1056  int pgrastpos[2] = {-1, -1};
1057  rt_raster rast[2] = {NULL};
1058  uint32_t bandindex[2] = {0};
1059  uint32_t hasbandindex[2] = {0};
1060  double distance = 0;
1061 
1062  uint32_t i;
1063  uint32_t j;
1064  uint32_t k;
1065  uint32_t numBands;
1066  int rtn;
1067  int result;
1068 
1069  for (i = 0, j = 0; i < set_count; i++) {
1070  /* pgrast is null, return null */
1071  if (PG_ARGISNULL(j)) {
1072  for (k = 0; k < i; k++) {
1073  rt_raster_destroy(rast[k]);
1074  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1075  }
1076  PG_RETURN_NULL();
1077  }
1078  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
1079  pgrastpos[i] = j;
1080  j++;
1081 
1082  /* raster */
1083  rast[i] = rt_raster_deserialize(pgrast[i], FALSE);
1084  if (!rast[i]) {
1085  for (k = 0; k <= i; k++) {
1086  if (k < i)
1087  rt_raster_destroy(rast[k]);
1088  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1089  }
1090  elog(ERROR, "RASTER_dfullywithin: Could not deserialize the %s raster", i < 1 ? "first" : "second");
1091  PG_RETURN_NULL();
1092  }
1093 
1094  /* numbands */
1095  numBands = rt_raster_get_num_bands(rast[i]);
1096  if (numBands < 1) {
1097  elog(NOTICE, "The %s raster provided has no bands", i < 1 ? "first" : "second");
1098  if (i > 0) i++;
1099  for (k = 0; k < i; k++) {
1100  rt_raster_destroy(rast[k]);
1101  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1102  }
1103  PG_RETURN_NULL();
1104  }
1105 
1106  /* band index */
1107  if (!PG_ARGISNULL(j)) {
1108  bandindex[i] = PG_GETARG_INT32(j);
1109  if (bandindex[i] < 1 || bandindex[i] > numBands) {
1110  elog(NOTICE, "Invalid band index (must use 1-based) for the %s raster. Returning NULL", i < 1 ? "first" : "second");
1111  if (i > 0) i++;
1112  for (k = 0; k < i; k++) {
1113  rt_raster_destroy(rast[k]);
1114  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1115  }
1116  PG_RETURN_NULL();
1117  }
1118  hasbandindex[i] = 1;
1119  }
1120  else
1121  hasbandindex[i] = 0;
1122  POSTGIS_RT_DEBUGF(4, "hasbandindex[%d] = %d", i, hasbandindex[i]);
1123  POSTGIS_RT_DEBUGF(4, "bandindex[%d] = %d", i, bandindex[i]);
1124  j++;
1125  }
1126 
1127  /* distance */
1128  if (PG_ARGISNULL(4)) {
1129  elog(NOTICE, "Distance cannot be NULL. Returning NULL");
1130  for (k = 0; k < set_count; k++) {
1131  rt_raster_destroy(rast[k]);
1132  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1133  }
1134  PG_RETURN_NULL();
1135  }
1136 
1137  distance = PG_GETARG_FLOAT8(4);
1138  if (distance < 0) {
1139  elog(NOTICE, "Distance cannot be less than zero. Returning NULL");
1140  for (k = 0; k < set_count; k++) {
1141  rt_raster_destroy(rast[k]);
1142  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1143  }
1144  PG_RETURN_NULL();
1145  }
1146 
1147  /* hasbandindex must be balanced */
1148  if (
1149  (hasbandindex[0] && !hasbandindex[1]) ||
1150  (!hasbandindex[0] && hasbandindex[1])
1151  ) {
1152  elog(NOTICE, "Missing band index. Band indices must be provided for both rasters if any one is provided");
1153  for (k = 0; k < set_count; k++) {
1154  rt_raster_destroy(rast[k]);
1155  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1156  }
1157  PG_RETURN_NULL();
1158  }
1159 
1160  /* SRID must match */
1161  if (rt_raster_get_srid(rast[0]) != rt_raster_get_srid(rast[1])) {
1162  for (k = 0; k < set_count; k++) {
1163  rt_raster_destroy(rast[k]);
1164  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1165  }
1166  elog(ERROR, "The two rasters provided have different SRIDs");
1167  PG_RETURN_NULL();
1168  }
1169 
1170  rtn = rt_raster_fully_within_distance(
1171  rast[0], (hasbandindex[0] ? (int)bandindex[0] - 1 : -1),
1172  rast[1], (hasbandindex[1] ? (int)bandindex[1] - 1 : -1),
1173  distance,
1174  &result
1175  );
1176  for (k = 0; k < set_count; k++) {
1177  rt_raster_destroy(rast[k]);
1178  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1179  }
1180 
1181  if (rtn != ES_NONE) {
1182  elog(ERROR, "RASTER_dfullywithin: Could not test that the two rasters are fully within the specified distance of each other");
1183  PG_RETURN_NULL();
1184  }
1185 
1186  PG_RETURN_BOOL(result);
1187 }
1188 
1192 PG_FUNCTION_INFO_V1(RASTER_sameAlignment);
1193 Datum RASTER_sameAlignment(PG_FUNCTION_ARGS)
1194 {
1195  const uint32_t set_count = 2;
1196  rt_pgraster *pgrast[2];
1197  int pgrastpos[2] = {-1, -1};
1198  rt_raster rast[2] = {NULL};
1199 
1200  uint32_t i;
1201  uint32_t j;
1202  uint32_t k;
1203  int rtn;
1204  int aligned = 0;
1205  char *reason = NULL;
1206 
1207  for (i = 0, j = 0; i < set_count; i++) {
1208  /* pgrast is null, return null */
1209  if (PG_ARGISNULL(j)) {
1210  for (k = 0; k < i; k++) {
1211  rt_raster_destroy(rast[k]);
1212  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1213  }
1214  PG_RETURN_NULL();
1215  }
1216  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(j), 0, sizeof(struct rt_raster_serialized_t));
1217  pgrastpos[i] = j;
1218  j++;
1219 
1220  /* raster */
1221  rast[i] = rt_raster_deserialize(pgrast[i], TRUE);
1222  if (!rast[i]) {
1223  for (k = 0; k <= i; k++) {
1224  if (k < i)
1225  rt_raster_destroy(rast[k]);
1226  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1227  }
1228  elog(ERROR, "RASTER_sameAlignment: Could not deserialize the %s raster", i < 1 ? "first" : "second");
1229  PG_RETURN_NULL();
1230  }
1231  }
1232 
1233  rtn = rt_raster_same_alignment(
1234  rast[0],
1235  rast[1],
1236  &aligned,
1237  &reason
1238  );
1239  for (k = 0; k < set_count; k++) {
1240  rt_raster_destroy(rast[k]);
1241  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1242  }
1243 
1244  if (rtn != ES_NONE) {
1245  elog(ERROR, "RASTER_sameAlignment: Could not test for alignment on the two rasters");
1246  PG_RETURN_NULL();
1247  }
1248 
1249  /* only output reason if not aligned */
1250  if (reason != NULL && !aligned)
1251  elog(NOTICE, "%s", reason);
1252 
1253  PG_RETURN_BOOL(aligned);
1254 }
1255 
1259 PG_FUNCTION_INFO_V1(RASTER_notSameAlignmentReason);
1260 Datum RASTER_notSameAlignmentReason(PG_FUNCTION_ARGS)
1261 {
1262  const uint32_t set_count = 2;
1263  rt_pgraster *pgrast[2];
1264  int pgrastpos[2] = {-1, -1};
1265  rt_raster rast[2] = {NULL};
1266 
1267  uint32_t i;
1268  uint32_t j;
1269  uint32_t k;
1270  int rtn;
1271  int aligned = 0;
1272  char *reason = NULL;
1273  text *result = NULL;
1274 
1275  for (i = 0, j = 0; i < set_count; i++) {
1276  /* pgrast is null, return null */
1277  if (PG_ARGISNULL(j)) {
1278  for (k = 0; k < i; k++) {
1279  rt_raster_destroy(rast[k]);
1280  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1281  }
1282  PG_RETURN_NULL();
1283  }
1284  pgrast[i] = (rt_pgraster *) PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(j), 0, sizeof(struct rt_raster_serialized_t));
1285  pgrastpos[i] = j;
1286  j++;
1287 
1288  /* raster */
1289  rast[i] = rt_raster_deserialize(pgrast[i], TRUE);
1290  if (!rast[i]) {
1291  for (k = 0; k <= i; k++) {
1292  if (k < i)
1293  rt_raster_destroy(rast[k]);
1294  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1295  }
1296  elog(ERROR, "RASTER_notSameAlignmentReason: Could not deserialize the %s raster", i < 1 ? "first" : "second");
1297  PG_RETURN_NULL();
1298  }
1299  }
1300 
1301  rtn = rt_raster_same_alignment(
1302  rast[0],
1303  rast[1],
1304  &aligned,
1305  &reason
1306  );
1307  for (k = 0; k < set_count; k++) {
1308  rt_raster_destroy(rast[k]);
1309  PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
1310  }
1311 
1312  if (rtn != ES_NONE) {
1313  elog(ERROR, "RASTER_notSameAlignmentReason: Could not test for alignment on the two rasters");
1314  PG_RETURN_NULL();
1315  }
1316 
1317  result = cstring_to_text(reason);
1318  PG_RETURN_TEXT_P(result);
1319 }
1320 
1321 
1322 /* calculate the proportion of each cell in a grid covered */
1323 /* by an input geometry */
1324 PG_FUNCTION_INFO_V1(RASTER_intersectionFractions);
1325 Datum RASTER_intersectionFractions(PG_FUNCTION_ARGS)
1326 {
1327 #if POSTGIS_GEOS_VERSION < 31400
1328  elog(ERROR, "The GEOS version this PostGIS binary "
1329  "was compiled against (%d) does not include the "
1330  "'GEOSGridIntersectionFractions' function (3.14.0+ required)",
1331  POSTGIS_GEOS_VERSION);
1332  PG_RETURN_NULL();
1333 #else
1334  rt_pgraster *pgrast_in = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1335  GSERIALIZED *gser = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
1336  LWGEOM *lwg = lwgeom_from_gserialized(gser);
1337  rt_raster rast_out = NULL;
1338  rt_pgraster *pgrast_out = NULL;
1339  rt_raster rast_in = rt_raster_deserialize(pgrast_in, FALSE);
1340 
1341  int gtype = lwgeom_get_type(lwg);
1342  if (gtype != POLYGONTYPE && gtype != MULTIPOLYGONTYPE &&
1343  gtype != LINETYPE && gtype != MULTILINETYPE)
1344  {
1345  elog(ERROR, "ST_IntersectionFractions: Unsupported geometry type '%s'", lwtype_name(gtype));
1346  }
1347 
1348  rast_out = rt_raster_intersection_fractions(rast_in, lwg);
1349  if (!rast_out)
1350  elog(ERROR, "ST_IntersectionFractions: calculation returned NULL");
1351 
1352  pgrast_out = rt_raster_serialize(rast_out);
1353  rt_raster_destroy(rast_out);
1354  if (!pgrast_out)
1355  PG_RETURN_NULL();
1356 
1357  SET_VARSIZE(pgrast_out, pgrast_out->size);
1358  PG_RETURN_POINTER(pgrast_out);
1359 
1360 #endif /* POSTGIS_GEOS_VERSION < 31400 */
1361 }
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition: cu_print.c:267
TRUE
#define TRUE
Definition: dbfopen.c:73
FALSE
#define FALSE
Definition: dbfopen.c:72
lwgeom_from_gserialized
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Definition: gserialized.c:268
MULTILINETYPE
#define MULTILINETYPE
Definition: liblwgeom.h:106
LINETYPE
#define LINETYPE
Definition: liblwgeom.h:103
MULTIPOLYGONTYPE
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:107
POLYGONTYPE
#define POLYGONTYPE
Definition: liblwgeom.h:104
lwtype_name
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:216
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_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:360
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_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:86
rt_raster_intersects
rt_errorstate rt_raster_intersects(rt_raster rast1, int nband1, rt_raster rast2, int nband2, int *intersects)
Return ES_ERROR if error occurred in function.
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_serialize
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
Definition: rt_serialize.c:521
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
ES_NONE
@ ES_NONE
Definition: librtcore.h:182
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_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:376
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_raster_deserialize
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:725
lwgeom_get_type
static uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwinline.h:141
distance
static double distance(double x1, double y1, double x2, double y2)
Definition: lwtree.c:1032
PG_FUNCTION_INFO_V1
PG_FUNCTION_INFO_V1(RASTER_intersects)
See if two rasters intersect.
RASTER_coveredby
Datum RASTER_coveredby(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:792
RASTER_sameAlignment
Datum RASTER_sameAlignment(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:1193
RASTER_containsProperly
Datum RASTER_containsProperly(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:554
RASTER_dfullywithin
Datum RASTER_dfullywithin(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:1052
RASTER_notSameAlignmentReason
Datum RASTER_notSameAlignmentReason(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:1260
RASTER_dwithin
Datum RASTER_dwithin(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:911
RASTER_contains
Datum RASTER_contains(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:435
RASTER_intersectionFractions
Datum RASTER_intersectionFractions(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:1325
RASTER_intersects
Datum RASTER_intersects(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:78
RASTER_covers
Datum RASTER_covers(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:673
RASTER_touches
Datum RASTER_touches(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:316
RASTER_overlaps
Datum RASTER_overlaps(PG_FUNCTION_ARGS)
Definition: rtpg_spatial_relationship.c:197
POSTGIS_RT_DEBUGF
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:69
POSTGIS_GEOS_VERSION
#define POSTGIS_GEOS_VERSION
Definition: sqldefines.h:11
rt_raster_serialized_t::size
uint32_t size
Definition: librtcore.h:2442
rt_raster_serialized_t
Struct definitions.
Definition: librtcore.h:2440
rt_raster_t
Definition: librtcore.h:2474