PostGIS  2.5.7dev-r@@SVN_REVISION@@
rt_geometry.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 #include "rt_serialize.h"
33 
34 /******************************************************************************
35 * rt_raster_perimeter()
36 ******************************************************************************/
37 
38 static rt_errorstate
40  uint16_t width = 0;
41  uint16_t height = 0;
42  int x = 0;
43  int y = 0;
44  int offset = 0;
45  int done[4] = {0};
46  double value = 0;
47  int nodata = 0;
48 
49  assert(band != NULL);
50  assert(band->raster != NULL);
51  assert(trim != NULL);
52 
53  memset(trim, 0, sizeof(uint16_t) * 4);
54 
55  width = rt_band_get_width(band);
56  height = rt_band_get_height(band);
57 
58  /* top */
59  for (y = 0; y < height; y++) {
60  for (offset = 0; offset < 3; offset++) {
61  /* every third pixel */
62  for (x = offset; x < width; x += 3) {
63  if (rt_band_get_pixel(band, x, y, &value, &nodata) != ES_NONE) {
64  rterror("_rti_raster_get_band_perimeter: Could not get band pixel");
65  return ES_ERROR;
66  }
67 
68  RASTER_DEBUGF(4, "top (x, y, value, nodata) = (%d, %d, %f, %d)", x, y, value, nodata);
69  if (!nodata) {
70  trim[0] = y;
71  done[0] = 1;
72  break;
73  }
74  }
75 
76  if (done[0])
77  break;
78  }
79 
80  if (done[0])
81  break;
82  }
83 
84  /* right */
85  for (x = width - 1; x >= 0; x--) {
86  for (offset = 0; offset < 3; offset++) {
87  /* every third pixel */
88  for (y = offset; y < height; y += 3) {
89  if (rt_band_get_pixel(band, x, y, &value, &nodata) != ES_NONE) {
90  rterror("_rti_raster_get_band_perimeter: Could not get band pixel");
91  return ES_ERROR;
92  }
93 
94  RASTER_DEBUGF(4, "right (x, y, value, nodata) = (%d, %d, %f, %d)", x, y, value, nodata);
95  if (!nodata) {
96  trim[1] = width - (x + 1);
97  done[1] = 1;
98  break;
99  }
100  }
101 
102  if (done[1])
103  break;
104  }
105 
106  if (done[1])
107  break;
108  }
109 
110  /* bottom */
111  for (y = height - 1; y >= 0; y--) {
112  for (offset = 0; offset < 3; offset++) {
113  /* every third pixel */
114  for (x = offset; x < width; x += 3) {
115  if (rt_band_get_pixel(band, x, y, &value, &nodata) != ES_NONE) {
116  rterror("_rti_raster_get_band_perimeter: Could not get band pixel");
117  return ES_ERROR;
118  }
119 
120  RASTER_DEBUGF(4, "bottom (x, y, value, nodata) = (%d, %d, %f, %d)", x, y, value, nodata);
121  if (!nodata) {
122  trim[2] = height - (y + 1);
123  done[2] = 1;
124  break;
125  }
126  }
127 
128  if (done[2])
129  break;
130  }
131 
132  if (done[2])
133  break;
134  }
135 
136  /* left */
137  for (x = 0; x < width; x++) {
138  for (offset = 0; offset < 3; offset++) {
139  /* every third pixel */
140  for (y = offset; y < height; y += 3) {
141  if (rt_band_get_pixel(band, x, y, &value, &nodata) != ES_NONE) {
142  rterror("_rti_raster_get_band_perimeter: Could not get band pixel");
143  return ES_ERROR;
144  }
145 
146  RASTER_DEBUGF(4, "left (x, , value, nodata) = (%d, %d, %f, %d)", x, y, value, nodata);
147  if (!nodata) {
148  trim[3] = x;
149  done[3] = 1;
150  break;
151  }
152  }
153 
154  if (done[3])
155  break;
156  }
157 
158  if (done[3])
159  break;
160  }
161 
162  RASTER_DEBUGF(4, "trim = (%d, %d, %d, %d)",
163  trim[0], trim[1], trim[2], trim[3]);
164 
165  return ES_NONE;
166 }
167 
183  rt_raster raster, int nband,
184  LWGEOM **perimeter
185 ) {
186  rt_band band = NULL;
187  int numband = 0;
188  uint16_t *_nband = NULL;
189  int i = 0;
190  int j = 0;
191  uint16_t _trim[4] = {0};
192  uint16_t trim[4] = {0}; /* top, right, bottom, left */
193  int isset[4] = {0};
194  double gt[6] = {0.0};
195  int srid = SRID_UNKNOWN;
196 
197  POINTARRAY *pts = NULL;
198  POINT4D p4d;
199  POINTARRAY **rings = NULL;
200  LWPOLY* poly = NULL;
201 
202  assert(perimeter != NULL);
203 
204  *perimeter = NULL;
205 
206  /* empty raster, no perimeter */
208  return ES_NONE;
209 
210  /* raster metadata */
211  srid = rt_raster_get_srid(raster);
213  numband = rt_raster_get_num_bands(raster);
214 
215  RASTER_DEBUGF(3, "rt_raster_get_perimeter: raster is %dx%d", raster->width, raster->height);
216 
217  /* nband < 0 means all bands */
218  if (nband >= 0) {
219  if (nband >= numband) {
220  rterror("rt_raster_get_boundary: Band %d not found for raster", nband);
221  return ES_ERROR;
222  }
223 
224  numband = 1;
225  }
226  else
227  nband = -1;
228 
229  RASTER_DEBUGF(3, "rt_raster_get_perimeter: nband, numband = %d, %d", nband, numband);
230 
231  _nband = rtalloc(sizeof(uint16_t) * numband);
232  if (_nband == NULL) {
233  rterror("rt_raster_get_boundary: Could not allocate memory for band indices");
234  return ES_ERROR;
235  }
236 
237  if (nband < 0) {
238  for (i = 0; i < numband; i++)
239  _nband[i] = i;
240  }
241  else
242  _nband[0] = nband;
243 
244  for (i = 0; i < numband; i++) {
245  band = rt_raster_get_band(raster, _nband[i]);
246  if (band == NULL) {
247  rterror("rt_raster_get_boundary: Could not get band at index %d", _nband[i]);
248  rtdealloc(_nband);
249  return ES_ERROR;
250  }
251 
252  /* band is nodata */
254  continue;
255 
257  rterror("rt_raster_get_boundary: Could not get band perimeter");
258  rtdealloc(_nband);
259  return ES_ERROR;
260  }
261 
262  for (j = 0; j < 4; j++) {
263  if (!isset[j] || trim[j] < _trim[j]) {
264  _trim[j] = trim[j];
265  isset[j] = 1;
266  }
267  }
268  }
269 
270  /* no longer needed */
271  rtdealloc(_nband);
272 
273  /* check isset, just need to check one element */
274  if (!isset[0]) {
275  /* return NULL as bands are empty */
276  return ES_NONE;
277  }
278 
279  RASTER_DEBUGF(4, "trim = (%d, %d, %d, %d)",
280  trim[0], trim[1], trim[2], trim[3]);
281 
282  /* only one ring */
283  rings = (POINTARRAY **) rtalloc(sizeof (POINTARRAY*));
284  if (!rings) {
285  rterror("rt_raster_get_perimeter: Could not allocate memory for polygon ring");
286  return ES_ERROR;
287  }
288  rings[0] = ptarray_construct(0, 0, 5);
289  if (!rings[0]) {
290  rterror("rt_raster_get_perimeter: Could not construct point array");
291  return ES_ERROR;
292  }
293  pts = rings[0];
294 
295  /* Upper-left corner (first and last points) */
297  raster,
298  _trim[3], _trim[0],
299  &p4d.x, &p4d.y,
300  gt
301  );
302  ptarray_set_point4d(pts, 0, &p4d);
303  ptarray_set_point4d(pts, 4, &p4d);
304 
305  /* Upper-right corner (we go clockwise) */
307  raster,
308  raster->width - _trim[1], _trim[0],
309  &p4d.x, &p4d.y,
310  gt
311  );
312  ptarray_set_point4d(pts, 1, &p4d);
313 
314  /* Lower-right corner */
316  raster,
317  raster->width - _trim[1], raster->height - _trim[2],
318  &p4d.x, &p4d.y,
319  gt
320  );
321  ptarray_set_point4d(pts, 2, &p4d);
322 
323  /* Lower-left corner */
325  raster,
326  _trim[3], raster->height - _trim[2],
327  &p4d.x, &p4d.y,
328  gt
329  );
330  ptarray_set_point4d(pts, 3, &p4d);
331 
332  poly = lwpoly_construct(srid, 0, 1, rings);
333  *perimeter = lwpoly_as_lwgeom(poly);
334 
335  return ES_NONE;
336 }
337 
338 /******************************************************************************
339 * rt_raster_surface()
340 ******************************************************************************/
341 
356  rt_band band = NULL;
357  LWGEOM *mpoly = NULL;
358  LWGEOM *tmp = NULL;
359  LWGEOM *clone = NULL;
360  rt_geomval gv = NULL;
361  int gvcount = 0;
362  GEOSGeometry *gc = NULL;
363  GEOSGeometry *gunion = NULL;
364  GEOSGeometry **geoms = NULL;
365  int geomscount = 0;
366  int i = 0;
367 
368  assert(surface != NULL);
369 
370  /* init *surface to NULL */
371  *surface = NULL;
372 
373  /* raster is empty, surface = NULL */
375  return ES_NONE;
376 
377  /* if nband < 0, return the convex hull as a multipolygon */
378  if (nband < 0) {
379  /*
380  lwgeom_as_multi() only does a shallow clone internally
381  so input and output geometries may share memory
382  hence the deep clone of the output geometry for returning
383  is the only way to guarentee the memory isn't shared
384  */
385  if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
386  rterror("rt_raster_surface: Could not get convex hull of raster");
387  return ES_ERROR;
388  }
389  mpoly = lwgeom_as_multi(tmp);
390  clone = lwgeom_clone_deep(mpoly);
391  lwgeom_free(tmp);
392  lwgeom_free(mpoly);
393 
394  *surface = lwgeom_as_lwmpoly(clone);
395  return ES_NONE;
396  }
397  /* check that nband is valid */
398  else if (nband >= rt_raster_get_num_bands(raster)) {
399  rterror("rt_raster_surface: The band index %d is invalid", nband);
400  return ES_ERROR;
401  }
402 
403  /* get band */
405  if (band == NULL) {
406  rterror("rt_raster_surface: Error getting band %d from raster", nband);
407  return ES_ERROR;
408  }
409 
410  /* band does not have a NODATA flag, return convex hull */
412  /*
413  lwgeom_as_multi() only does a shallow clone internally
414  so input and output geometries may share memory
415  hence the deep clone of the output geometry for returning
416  is the only way to guarentee the memory isn't shared
417  */
418  if (rt_raster_get_convex_hull(raster, &tmp) != ES_NONE) {
419  rterror("rt_raster_surface: Could not get convex hull of raster");
420  return ES_ERROR;
421  }
422  mpoly = lwgeom_as_multi(tmp);
423  clone = lwgeom_clone_deep(mpoly);
424  lwgeom_free(tmp);
425  lwgeom_free(mpoly);
426 
427  *surface = lwgeom_as_lwmpoly(clone);
428  return ES_NONE;
429  }
430  /* band is NODATA, return NULL */
431  else if (rt_band_get_isnodata_flag(band)) {
432  RASTER_DEBUG(3, "Band is NODATA. Returning NULL");
433  return ES_NONE;
434  }
435 
436  /* initialize GEOS */
437  initGEOS(rtinfo, lwgeom_geos_error);
438 
439  /* use gdal polygonize */
440  gv = rt_raster_gdal_polygonize(raster, nband, 1, &gvcount);
441  /* no polygons returned */
442  if (gvcount < 1) {
443  RASTER_DEBUG(3, "All pixels of band are NODATA. Returning NULL");
444  if (gv != NULL) rtdealloc(gv);
445  return ES_NONE;
446  }
447  /* more than 1 polygon */
448  else if (gvcount > 1) {
449  /* convert LWPOLY to GEOSGeometry */
450  geomscount = gvcount;
451  geoms = rtalloc(sizeof(GEOSGeometry *) * geomscount);
452  if (geoms == NULL) {
453  rterror("rt_raster_surface: Could not allocate memory for pixel polygons as GEOSGeometry");
454  for (i = 0; i < gvcount; i++) lwpoly_free(gv[i].geom);
455  rtdealloc(gv);
456  return ES_ERROR;
457  }
458  for (i = 0; i < gvcount; i++) {
459 #if POSTGIS_DEBUG_LEVEL > 3
460  {
461  char *wkt = lwgeom_to_wkt(lwpoly_as_lwgeom(gv[i].geom), WKT_ISO, DBL_DIG, NULL);
462  RASTER_DEBUGF(4, "geom %d = %s", i, wkt);
463  rtdealloc(wkt);
464  }
465 #endif
466 
467  geoms[i] = LWGEOM2GEOS(lwpoly_as_lwgeom(gv[i].geom), 0);
468  lwpoly_free(gv[i].geom);
469  }
470  rtdealloc(gv);
471 
472  /* create geometry collection */
473  gc = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, geomscount);
474 
475  if (gc == NULL) {
476  rterror("rt_raster_surface: Could not create GEOS GEOMETRYCOLLECTION from set of pixel polygons");
477 
478  for (i = 0; i < geomscount; i++)
479  GEOSGeom_destroy(geoms[i]);
480  rtdealloc(geoms);
481  return ES_ERROR;
482  }
483 
484  /* run the union */
485  gunion = GEOSUnaryUnion(gc);
486 
487  GEOSGeom_destroy(gc);
488  rtdealloc(geoms);
489 
490  if (gunion == NULL) {
491  rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnaryUnion()");
492  return ES_ERROR;
493  }
494 
495  /* convert union result from GEOSGeometry to LWGEOM */
496  mpoly = GEOS2LWGEOM(gunion, 0);
497 
498  /*
499  is geometry valid?
500  if not, try to make valid
501  */
502  do {
503  LWGEOM *mpolyValid = NULL;
504 
505  if (GEOSisValid(gunion))
506  break;
507 
508  /* make geometry valid */
509  mpolyValid = lwgeom_make_valid(mpoly);
510  if (mpolyValid == NULL) {
511  rtwarn("Cannot fix invalid geometry");
512  break;
513  }
514 
515  lwgeom_free(mpoly);
516  mpoly = mpolyValid;
517  }
518  while (0);
519 
520  GEOSGeom_destroy(gunion);
521  }
522  else {
523  mpoly = lwpoly_as_lwgeom(gv[0].geom);
524  rtdealloc(gv);
525 
526 #if POSTGIS_DEBUG_LEVEL > 3
527  {
528  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
529  RASTER_DEBUGF(4, "geom 0 = %s", wkt);
530  rtdealloc(wkt);
531  }
532 #endif
533  }
534 
535  /* specify SRID */
537 
538  if (mpoly != NULL) {
539  /* convert to multi */
540  if (!lwgeom_is_collection(mpoly)) {
541  tmp = mpoly;
542 
543 #if POSTGIS_DEBUG_LEVEL > 3
544  {
545  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
546  RASTER_DEBUGF(4, "before multi = %s", wkt);
547  rtdealloc(wkt);
548  }
549 #endif
550 
551  RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
552 
553  /*
554  lwgeom_as_multi() only does a shallow clone internally
555  so input and output geometries may share memory
556  hence the deep clone of the output geometry for returning
557  is the only way to guarentee the memory isn't shared
558  */
559  mpoly = lwgeom_as_multi(tmp);
560  clone = lwgeom_clone_deep(mpoly);
561  lwgeom_free(tmp);
562  lwgeom_free(mpoly);
563  mpoly = clone;
564 
565  RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
566 
567 #if POSTGIS_DEBUG_LEVEL > 3
568  {
569  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
570  RASTER_DEBUGF(4, "after multi = %s", wkt);
571  rtdealloc(wkt);
572  }
573 #endif
574  }
575 
576 #if POSTGIS_DEBUG_LEVEL > 3
577  {
578  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
579  RASTER_DEBUGF(4, "returning geometry = %s", wkt);
580  rtdealloc(wkt);
581  }
582 #endif
583 
584  *surface = lwgeom_as_lwmpoly(mpoly);
585  return ES_NONE;
586  }
587 
588  return ES_NONE;
589 }
590 
591 /******************************************************************************
592 * rt_raster_pixel_as_polygon()
593 ******************************************************************************/
594 
609 LWPOLY*
611 {
612  double scale_x, scale_y;
613  double skew_x, skew_y;
614  double ul_x, ul_y;
615  int srid;
616  POINTARRAY **points;
617  POINT4D p, p0;
618  LWPOLY *poly;
619 
620  assert(rast != NULL);
621 
622  scale_x = rt_raster_get_x_scale(rast);
623  scale_y = rt_raster_get_y_scale(rast);
624  skew_x = rt_raster_get_x_skew(rast);
625  skew_y = rt_raster_get_y_skew(rast);
628  srid = rt_raster_get_srid(rast);
629 
630  points = rtalloc(sizeof(POINTARRAY *)*1);
631  points[0] = ptarray_construct(0, 0, 5);
632 
633  p0.x = scale_x * x + skew_x * y + ul_x;
634  p0.y = scale_y * y + skew_y * x + ul_y;
635  ptarray_set_point4d(points[0], 0, &p0);
636 
637  p.x = p0.x + scale_x;
638  p.y = p0.y + skew_y;
639  ptarray_set_point4d(points[0], 1, &p);
640 
641  p.x = p0.x + scale_x + skew_x;
642  p.y = p0.y + scale_y + skew_y;
643  ptarray_set_point4d(points[0], 2, &p);
644 
645  p.x = p0.x + skew_x;
646  p.y = p0.y + scale_y;
647  ptarray_set_point4d(points[0], 3, &p);
648 
649  /* close it */
650  ptarray_set_point4d(points[0], 4, &p0);
651 
652  poly = lwpoly_construct(srid, NULL, 1, points);
653 
654  return poly;
655 }
656 
657 /******************************************************************************
658 * rt_raster_get_envelope_geom()
659 ******************************************************************************/
660 
671  double gt[6] = {0.0};
672  int srid = SRID_UNKNOWN;
673 
674  POINTARRAY *pts = NULL;
675  POINT4D p4d;
676 
677  assert(env != NULL);
678  *env = NULL;
679 
680  /* raster is NULL, envelope is NULL */
681  if (raster == NULL)
682  return ES_NONE;
683 
684  /* raster metadata */
685  srid = rt_raster_get_srid(raster);
687 
689  3,
690  "rt_raster_get_envelope: raster is %dx%d",
691  raster->width,
692  raster->height
693  );
694 
695  /* return point or line since at least one of the two dimensions is 0 */
696  if ((!raster->width) || (!raster->height)) {
697  p4d.x = gt[0];
698  p4d.y = gt[3];
699 
700  /* return point */
701  if (!raster->width && !raster->height) {
702  LWPOINT *point = lwpoint_make2d(srid, p4d.x, p4d.y);
703  *env = lwpoint_as_lwgeom(point);
704  }
705  /* return linestring */
706  else {
707  LWLINE *line = NULL;
708  pts = ptarray_construct_empty(0, 0, 2);
709 
710  /* first point of line */
711  ptarray_append_point(pts, &p4d, LW_TRUE);
712 
713  /* second point of line */
715  raster,
717  &p4d.x, &p4d.y,
718  gt
719  ) != ES_NONE) {
720  rterror("rt_raster_get_envelope: Could not get second point for linestring");
721  return ES_ERROR;
722  }
723  ptarray_append_point(pts, &p4d, LW_TRUE);
724  line = lwline_construct(srid, NULL, pts);
725 
726  *env = lwline_as_lwgeom(line);
727  }
728 
729  return ES_NONE;
730  }
731  else {
732  rt_envelope rtenv;
733  int err = ES_NONE;
734  POINTARRAY **rings = NULL;
735  LWPOLY* poly = NULL;
736 
737  /* only one ring */
738  rings = (POINTARRAY **) rtalloc(sizeof (POINTARRAY*));
739  if (!rings) {
740  rterror("rt_raster_get_envelope_geom: Could not allocate memory for polygon ring");
741  return ES_ERROR;
742  }
743  rings[0] = ptarray_construct(0, 0, 5);
744  if (!rings[0]) {
745  rterror("rt_raster_get_envelope_geom: Could not construct point array");
746  return ES_ERROR;
747  }
748  pts = rings[0];
749 
750  err = rt_raster_get_envelope(raster, &rtenv);
751  if (err != ES_NONE) {
752  rterror("rt_raster_get_envelope_geom: Could not get raster envelope");
753  return err;
754  }
755 
756  /* build ring */
757 
758  /* minx, maxy */
759  p4d.x = rtenv.MinX;
760  p4d.y = rtenv.MaxY;
761  ptarray_set_point4d(pts, 0, &p4d);
762  ptarray_set_point4d(pts, 4, &p4d);
763 
764  /* maxx, maxy */
765  p4d.x = rtenv.MaxX;
766  p4d.y = rtenv.MaxY;
767  ptarray_set_point4d(pts, 1, &p4d);
768 
769  /* maxx, miny */
770  p4d.x = rtenv.MaxX;
771  p4d.y = rtenv.MinY;
772  ptarray_set_point4d(pts, 2, &p4d);
773 
774  /* minx, miny */
775  p4d.x = rtenv.MinX;
776  p4d.y = rtenv.MinY;
777  ptarray_set_point4d(pts, 3, &p4d);
778 
779  poly = lwpoly_construct(srid, 0, 1, rings);
780  *env = lwpoly_as_lwgeom(poly);
781  }
782 
783  return ES_NONE;
784 }
785 
786 /******************************************************************************
787 * rt_raster_get_convex_hull()
788 ******************************************************************************/
789 
804  double gt[6] = {0.0};
805  int srid = SRID_UNKNOWN;
806 
807  POINTARRAY *pts = NULL;
808  POINT4D p4d;
809 
810  assert(hull != NULL);
811  *hull = NULL;
812 
813  /* raster is NULL, convex hull is NULL */
814  if (raster == NULL)
815  return ES_NONE;
816 
817  /* raster metadata */
818  srid = rt_raster_get_srid(raster);
820 
821  RASTER_DEBUGF(3, "rt_raster_get_convex_hull: raster is %dx%d", raster->width, raster->height);
822 
823  /* return point or line since at least one of the two dimensions is 0 */
824  if ((!raster->width) || (!raster->height)) {
825  p4d.x = gt[0];
826  p4d.y = gt[3];
827 
828  /* return point */
829  if (!raster->width && !raster->height) {
830  LWPOINT *point = lwpoint_make2d(srid, p4d.x, p4d.y);
831  *hull = lwpoint_as_lwgeom(point);
832  }
833  /* return linestring */
834  else {
835  LWLINE *line = NULL;
836  pts = ptarray_construct_empty(0, 0, 2);
837 
838  /* first point of line */
839  ptarray_append_point(pts, &p4d, LW_TRUE);
840 
841  /* second point of line */
843  raster,
845  &p4d.x, &p4d.y,
846  gt
847  ) != ES_NONE) {
848  rterror("rt_raster_get_convex_hull: Could not get second point for linestring");
849  return ES_ERROR;
850  }
851  ptarray_append_point(pts, &p4d, LW_TRUE);
852  line = lwline_construct(srid, NULL, pts);
853 
854  *hull = lwline_as_lwgeom(line);
855  }
856 
857  return ES_NONE;
858  }
859  else {
860  POINTARRAY **rings = NULL;
861  LWPOLY* poly = NULL;
862 
863  /* only one ring */
864  rings = (POINTARRAY **) rtalloc(sizeof (POINTARRAY*));
865  if (!rings) {
866  rterror("rt_raster_get_convex_hull: Could not allocate memory for polygon ring");
867  return ES_ERROR;
868  }
869  rings[0] = ptarray_construct(0, 0, 5);
870  /* TODO: handle error on ptarray construction */
871  /* XXX jorgearevalo: the error conditions aren't managed in ptarray_construct */
872  if (!rings[0]) {
873  rterror("rt_raster_get_convex_hull: Could not construct point array");
874  return ES_ERROR;
875  }
876  pts = rings[0];
877 
878  /* Upper-left corner (first and last points) */
879  p4d.x = gt[0];
880  p4d.y = gt[3];
881  ptarray_set_point4d(pts, 0, &p4d);
882  ptarray_set_point4d(pts, 4, &p4d);
883 
884  /* Upper-right corner (we go clockwise) */
886  raster,
887  raster->width, 0,
888  &p4d.x, &p4d.y,
889  gt
890  );
891  ptarray_set_point4d(pts, 1, &p4d);
892 
893  /* Lower-right corner */
895  raster,
896  raster->width, raster->height,
897  &p4d.x, &p4d.y,
898  gt
899  );
900  ptarray_set_point4d(pts, 2, &p4d);
901 
902  /* Lower-left corner */
904  raster,
905  0, raster->height,
906  &p4d.x, &p4d.y,
907  gt
908  );
909  ptarray_set_point4d(pts, 3, &p4d);
910 
911  poly = lwpoly_construct(srid, 0, 1, rings);
912  *hull = lwpoly_as_lwgeom(poly);
913  }
914 
915  return ES_NONE;
916 }
917 
918 /******************************************************************************
919 * rt_raster_gdal_polygonize()
920 ******************************************************************************/
921 
941  rt_raster raster, int nband,
942  int exclude_nodata_value,
943  int *pnElements
944 ) {
945  CPLErr cplerr = CE_None;
946  char *pszQuery;
947  long j;
948  OGRSFDriverH ogr_drv = NULL;
949  GDALDriverH gdal_drv = NULL;
950  int destroy_gdal_drv = 0;
951  GDALDatasetH memdataset = NULL;
952  GDALRasterBandH gdal_band = NULL;
953  OGRDataSourceH memdatasource = NULL;
954  rt_geomval pols = NULL;
955  OGRLayerH hLayer = NULL;
956  OGRFeatureH hFeature = NULL;
957  OGRGeometryH hGeom = NULL;
958  OGRFieldDefnH hFldDfn = NULL;
959  unsigned char *wkb = NULL;
960  int wkbsize = 0;
961  LWGEOM *lwgeom = NULL;
962  int nFeatureCount = 0;
963  rt_band band = NULL;
964  int iPixVal = -1;
965  double dValue = 0.0;
966  int iBandHasNodataValue = FALSE;
967  double dBandNoData = 0.0;
968 
969  /* for checking that a geometry is valid */
970  GEOSGeometry *ggeom = NULL;
971  int isValid;
972  LWGEOM *lwgeomValid = NULL;
973 
974  uint32_t bandNums[1] = {nband};
975  int excludeNodataValues[1] = {exclude_nodata_value};
976 
977  /* checks */
978  assert(NULL != raster);
979  assert(NULL != pnElements);
980 
981  RASTER_DEBUG(2, "In rt_raster_gdal_polygonize");
982 
983  *pnElements = 0;
984 
985  /*******************************
986  * Get band
987  *******************************/
989  if (NULL == band) {
990  rterror("rt_raster_gdal_polygonize: Error getting band %d from raster", nband);
991  return NULL;
992  }
993 
994  if (exclude_nodata_value) {
995 
996  /* band is NODATA */
998  RASTER_DEBUG(3, "Band is NODATA. Returning null");
999  *pnElements = 0;
1000  return NULL;
1001  }
1002 
1003  iBandHasNodataValue = rt_band_get_hasnodata_flag(band);
1004  if (iBandHasNodataValue)
1005  rt_band_get_nodata(band, &dBandNoData);
1006  else
1007  exclude_nodata_value = FALSE;
1008  }
1009 
1010  /*****************************************************
1011  * Convert raster to GDAL MEM dataset
1012  *****************************************************/
1013  memdataset = rt_raster_to_gdal_mem(raster, NULL, bandNums, excludeNodataValues, 1, &gdal_drv, &destroy_gdal_drv);
1014  if (NULL == memdataset) {
1015  rterror("rt_raster_gdal_polygonize: Couldn't convert raster to GDAL MEM dataset");
1016  return NULL;
1017  }
1018 
1019  /*****************************
1020  * Register ogr mem driver
1021  *****************************/
1022 #ifdef GDAL_DCAP_RASTER
1023  /* in GDAL 2.0, OGRRegisterAll() is an alias to GDALAllRegister() */
1025 #else
1026  OGRRegisterAll();
1027 #endif
1028 
1029  RASTER_DEBUG(3, "creating OGR MEM vector");
1030 
1031  /*****************************************************
1032  * Create an OGR in-memory vector for layers
1033  *****************************************************/
1034  ogr_drv = OGRGetDriverByName("Memory");
1035  memdatasource = OGR_Dr_CreateDataSource(ogr_drv, "", NULL);
1036  if (NULL == memdatasource) {
1037  rterror("rt_raster_gdal_polygonize: Couldn't create a OGR Datasource to store pols");
1038  GDALClose(memdataset);
1039  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1040  return NULL;
1041  }
1042 
1043  /* Can MEM driver create new layers? */
1044  if (!OGR_DS_TestCapability(memdatasource, ODsCCreateLayer)) {
1045  rterror("rt_raster_gdal_polygonize: MEM driver can't create new layers, aborting");
1046 
1047  /* xxx jorgearevalo: what should we do now? */
1048  GDALClose(memdataset);
1049  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1050  OGRReleaseDataSource(memdatasource);
1051 
1052  return NULL;
1053  }
1054 
1055  RASTER_DEBUG(3, "polygonizying GDAL MEM raster band");
1056 
1057  /*****************************
1058  * Polygonize the raster band
1059  *****************************/
1060 
1066  hLayer = OGR_DS_CreateLayer(memdatasource, "PolygonizedLayer", NULL, wkbPolygon, NULL);
1067 
1068  if (NULL == hLayer) {
1069  rterror("rt_raster_gdal_polygonize: Couldn't create layer to store polygons");
1070 
1071  GDALClose(memdataset);
1072  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1073  OGRReleaseDataSource(memdatasource);
1074 
1075  return NULL;
1076  }
1077 
1082  /* First, create a field definition to create the field */
1083  hFldDfn = OGR_Fld_Create("PixelValue", OFTReal);
1084 
1085  /* Second, create the field */
1086  if (OGR_L_CreateField(hLayer, hFldDfn, TRUE) != OGRERR_NONE) {
1087  rtwarn("Couldn't create a field in OGR Layer. The polygons generated won't be able to store the pixel value");
1088  iPixVal = -1;
1089  }
1090  else {
1091  /* Index to the new field created in the layer */
1092  iPixVal = 0;
1093  }
1094 
1095  /* Get GDAL raster band */
1096  gdal_band = GDALGetRasterBand(memdataset, 1);
1097  if (NULL == gdal_band) {
1098  rterror("rt_raster_gdal_polygonize: Couldn't get GDAL band to polygonize");
1099 
1100  GDALClose(memdataset);
1101  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1102  OGR_Fld_Destroy(hFldDfn);
1103  OGR_DS_DeleteLayer(memdatasource, 0);
1104  OGRReleaseDataSource(memdatasource);
1105 
1106  return NULL;
1107  }
1108 
1109  /* We don't need a raster mask band. Each band has a nodata value. */
1110  cplerr = GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1111 
1112  if (cplerr != CE_None) {
1113  rterror("rt_raster_gdal_polygonize: Could not polygonize GDAL band");
1114 
1115  GDALClose(memdataset);
1116  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1117  OGR_Fld_Destroy(hFldDfn);
1118  OGR_DS_DeleteLayer(memdatasource, 0);
1119  OGRReleaseDataSource(memdatasource);
1120 
1121  return NULL;
1122  }
1123 
1130  if (iBandHasNodataValue) {
1131  pszQuery = (char *) rtalloc(50 * sizeof (char));
1132  sprintf(pszQuery, "PixelValue != %f", dBandNoData );
1133  OGRErr e = OGR_L_SetAttributeFilter(hLayer, pszQuery);
1134  if (e != OGRERR_NONE) {
1135  rtwarn("Error filtering NODATA values for band. All values will be treated as data values");
1136  }
1137  }
1138  else {
1139  pszQuery = NULL;
1140  }
1141 
1142  /*********************************************************************
1143  * Transform OGR layers to WKB polygons
1144  * XXX jorgearevalo: GDALPolygonize does not set the coordinate system
1145  * on the output layer. Application code should do this when the layer
1146  * is created, presumably matching the raster coordinate system.
1147  *********************************************************************/
1148  nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);
1149 
1150  /* Allocate memory for pols */
1151  pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t));
1152 
1153  if (NULL == pols) {
1154  rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set");
1155 
1156  GDALClose(memdataset);
1157  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1158  OGR_Fld_Destroy(hFldDfn);
1159  OGR_DS_DeleteLayer(memdatasource, 0);
1160  if (NULL != pszQuery)
1161  rtdealloc(pszQuery);
1162  OGRReleaseDataSource(memdatasource);
1163 
1164  return NULL;
1165  }
1166 
1167  /* initialize GEOS */
1168  initGEOS(rtinfo, lwgeom_geos_error);
1169 
1170  RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount);
1171 
1172  /* Reset feature reading to start in the first feature */
1173  OGR_L_ResetReading(hLayer);
1174 
1175  for (j = 0; j < nFeatureCount; j++) {
1176  hFeature = OGR_L_GetNextFeature(hLayer);
1177  dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
1178 
1179  hGeom = OGR_F_GetGeometryRef(hFeature);
1180  wkbsize = OGR_G_WkbSize(hGeom);
1181 
1182  /* allocate wkb buffer */
1183  wkb = rtalloc(sizeof(unsigned char) * wkbsize);
1184  if (wkb == NULL) {
1185  rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer");
1186 
1187  OGR_F_Destroy(hFeature);
1188  GDALClose(memdataset);
1189  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1190  OGR_Fld_Destroy(hFldDfn);
1191  OGR_DS_DeleteLayer(memdatasource, 0);
1192  if (NULL != pszQuery)
1193  rtdealloc(pszQuery);
1194  OGRReleaseDataSource(memdatasource);
1195 
1196  return NULL;
1197  }
1198 
1199  /* export WKB with LSB byte order */
1200  OGR_G_ExportToWkb(hGeom, wkbNDR, wkb);
1201 
1202  /* convert WKB to LWGEOM */
1203  lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE);
1204 
1205 #if POSTGIS_DEBUG_LEVEL > 3
1206  {
1207  char *wkt = NULL;
1208  OGR_G_ExportToWkt(hGeom, &wkt);
1209  RASTER_DEBUGF(4, "GDAL wkt = %s", wkt);
1210  CPLFree(wkt);
1211 
1212  d_print_binary_hex("GDAL wkb", wkb, wkbsize);
1213  }
1214 #endif
1215 
1216  /* cleanup unnecessary stuff */
1217  rtdealloc(wkb);
1218  wkb = NULL;
1219  wkbsize = 0;
1220 
1221  OGR_F_Destroy(hFeature);
1222 
1223  /* specify SRID */
1225 
1226  /*
1227  is geometry valid?
1228  if not, try to make valid
1229  */
1230  do {
1231  ggeom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
1232  if (ggeom == NULL) {
1233  rtwarn("Cannot test geometry for validity");
1234  break;
1235  }
1236 
1237  isValid = GEOSisValid(ggeom);
1238 
1239  GEOSGeom_destroy(ggeom);
1240  ggeom = NULL;
1241 
1242  /* geometry is valid */
1243  if (isValid)
1244  break;
1245 
1246  RASTER_DEBUG(3, "fixing invalid geometry");
1247 
1248  /* make geometry valid */
1249  lwgeomValid = lwgeom_make_valid(lwgeom);
1250  if (lwgeomValid == NULL) {
1251  rtwarn("Cannot fix invalid geometry");
1252  break;
1253  }
1254 
1255  lwgeom_free(lwgeom);
1256  lwgeom = lwgeomValid;
1257  }
1258  while (0);
1259 
1260  /* save lwgeom */
1261  pols[j].geom = lwgeom_as_lwpoly(lwgeom);
1262 
1263 #if POSTGIS_DEBUG_LEVEL > 3
1264  {
1265  char *wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, NULL);
1266  RASTER_DEBUGF(4, "LWGEOM wkt = %s", wkt);
1267  rtdealloc(wkt);
1268 
1269  size_t lwwkbsize = 0;
1270  uint8_t *lwwkb = lwgeom_to_wkb(lwgeom, WKB_ISO | WKB_NDR, &lwwkbsize);
1271  if (lwwkbsize) {
1272  d_print_binary_hex("LWGEOM wkb", lwwkb, lwwkbsize);
1273  rtdealloc(lwwkb);
1274  }
1275  }
1276 #endif
1277 
1278  /* set pixel value */
1279  pols[j].val = dValue;
1280  }
1281 
1282  *pnElements = nFeatureCount;
1283 
1284  RASTER_DEBUG(3, "destroying GDAL MEM raster");
1285  GDALClose(memdataset);
1286  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1287 
1288  RASTER_DEBUG(3, "destroying OGR MEM vector");
1289  OGR_Fld_Destroy(hFldDfn);
1290  OGR_DS_DeleteLayer(memdatasource, 0);
1291  if (NULL != pszQuery) rtdealloc(pszQuery);
1292  OGRReleaseDataSource(memdatasource);
1293 
1294  return pols;
1295 }
1296 
#define TRUE
Definition: dbfopen.c:169
#define FALSE
Definition: dbfopen.c:168
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, uint8_t want3d)
void lwgeom_geos_error(const char *fmt,...)
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:330
#define WKB_ISO
Definition: liblwgeom.h:2066
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:2005
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition: lwgeom.c:371
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:320
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:520
uint8_t * lwgeom_to_wkb(const LWGEOM *geom, uint8_t variant, size_t *size_out)
Convert LWGEOM to a char* in WKB format.
Definition: lwout_wkb.c:764
LWMPOLY * lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:251
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:62
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:335
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:1085
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
#define WKT_ISO
Definition: liblwgeom.h:2075
#define WKB_NDR
Definition: liblwgeom.h:2069
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:676
LWPOLY * lwpoly_construct(int srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points)
Definition: lwpoly.c:43
int ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int allow_duplicates)
Append a point to the end of an existing POINTARRAY If allow_duplicate is LW_FALSE,...
Definition: ptarray.c:156
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:175
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:76
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:206
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:188
LWGEOM * lwgeom_from_wkb(const uint8_t *wkb, const size_t wkb_size, const char check)
WKB inputs must have a declared size, to prevent malformed WKB from reading off the end of the memory...
Definition: lwin_wkb.c:789
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
void lwgeom_set_srid(LWGEOM *geom, int srid)
Set the SRID on an LWGEOM For collections, only the parent gets an SRID, all the children get SRID_UN...
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:435
LWGEOM * lwgeom_make_valid(LWGEOM *geom)
Attempts to make an invalid geometries valid w/out losing points.
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_band.c:640
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
rt_errorstate rt_raster_cell_to_geopoint(rt_raster raster, double xr, double yr, double *xw, double *yw, double *gt)
Convert an xr, yr raster point to an xw, yw point on map.
Definition: rt_raster.c:755
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:356
int rt_util_gdal_register_all(int force_register_all)
Definition: rt_util.c:336
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_raster.c:181
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition: rt_raster.c:213
void rtinfo(const char *fmt,...)
Definition: rt_context.c:211
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:674
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1221
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:714
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
rt_errorstate
Enum definitions.
Definition: librtcore.h:179
@ ES_NONE
Definition: librtcore.h:180
@ ES_ERROR
Definition: librtcore.h:181
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
struct rt_geomval_t * rt_geomval
Definition: librtcore.h:149
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1730
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
void rtdealloc(void *mem)
Definition: rt_context.c:186
GDALDatasetH rt_raster_to_gdal_mem(rt_raster raster, const char *srs, uint32_t *bandNums, int *excludeNodataValues, int count, GDALDriverH *rtn_drv, int *destroy_rtn_drv)
Return GDAL dataset using GDAL MEM driver from raster.
Definition: rt_raster.c:1826
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159
rt_errorstate rt_raster_get_envelope(rt_raster raster, rt_envelope *env)
Get raster's envelope.
Definition: rt_raster.c:873
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition: rt_raster.c:190
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:706
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_band.c:649
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition: rt_raster.c:222
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1334
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
This library is the generic raster handling section of PostGIS.
int value
Definition: genraster.py:61
band
Definition: ovdump.py:57
nband
Definition: pixval.py:52
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
gt
Definition: window.py:77
static char * trim(const char *input)
Definition: raster2pgsql.c:262
rt_errorstate rt_raster_surface(rt_raster raster, int nband, LWMPOLY **surface)
Get a raster as a surface (multipolygon).
Definition: rt_geometry.c:355
rt_geomval rt_raster_gdal_polygonize(rt_raster raster, int nband, int exclude_nodata_value, int *pnElements)
Returns a set of "geomval" value, one for each group of pixel sharing the same value for the provided...
Definition: rt_geometry.c:940
LWPOLY * rt_raster_pixel_as_polygon(rt_raster rast, int x, int y)
Get a raster pixel as a polygon.
Definition: rt_geometry.c:610
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition: rt_geometry.c:803
static rt_errorstate _rti_raster_get_band_perimeter(rt_band band, uint16_t *trim)
Definition: rt_geometry.c:39
rt_errorstate rt_raster_get_envelope_geom(rt_raster raster, LWGEOM **env)
Get raster's envelope as a geometry.
Definition: rt_geometry.c:670
rt_errorstate rt_raster_get_perimeter(rt_raster raster, int nband, LWGEOM **perimeter)
Get raster perimeter.
Definition: rt_geometry.c:182
double x
Definition: liblwgeom.h:355
double y
Definition: liblwgeom.h:355
double MinX
Definition: librtcore.h:165
double MaxX
Definition: librtcore.h:166
double MinY
Definition: librtcore.h:167
double MaxY
Definition: librtcore.h:168
double val
Definition: librtcore.h:2354
LWPOLY * geom
Definition: librtcore.h:2353
unsigned int uint32_t
Definition: uthash.h:78
unsigned char uint8_t
Definition: uthash.h:79