PostGIS  2.2.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@keybit.net>
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 */
207  if (rt_raster_is_empty(raster))
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 */
253  if (rt_band_get_isnodata_flag(band) != 0)
254  continue;
255 
256  if (_rti_raster_get_band_perimeter(band, trim) != ES_NONE) {
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 */
374  if (rt_raster_is_empty(raster))
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 */
404  band = rt_raster_get_band(raster, nband);
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 */
411  if (!rt_band_get_hasnodata_flag(band)) {
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 #if POSTGIS_GEOS_VERSION >= 33
474  gc = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, geomscount);
475 #else
476  gc = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, geomscount);
477 #endif
478 
479  if (gc == NULL) {
480 #if POSTGIS_GEOS_VERSION >= 33
481  rterror("rt_raster_surface: Could not create GEOS GEOMETRYCOLLECTION from set of pixel polygons");
482 #else
483  rterror("rt_raster_surface: Could not create GEOS MULTIPOLYGON from set of pixel polygons");
484 #endif
485 
486  for (i = 0; i < geomscount; i++)
487  GEOSGeom_destroy(geoms[i]);
488  rtdealloc(geoms);
489  return ES_ERROR;
490  }
491 
492  /* run the union */
493 #if POSTGIS_GEOS_VERSION >= 33
494  gunion = GEOSUnaryUnion(gc);
495 #else
496  gunion = GEOSUnionCascaded(gc);
497 #endif
498  GEOSGeom_destroy(gc);
499  rtdealloc(geoms);
500 
501  if (gunion == NULL) {
502 #if POSTGIS_GEOS_VERSION >= 33
503  rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnaryUnion()");
504 #else
505  rterror("rt_raster_surface: Could not union the pixel polygons using GEOSUnionCascaded()");
506 #endif
507  return ES_ERROR;
508  }
509 
510  /* convert union result from GEOSGeometry to LWGEOM */
511  mpoly = GEOS2LWGEOM(gunion, 0);
512 
513  /*
514  is geometry valid?
515  if not, try to make valid
516  */
517  do {
518  LWGEOM *mpolyValid = NULL;
519 
520 #if POSTGIS_GEOS_VERSION < 33
521  break;
522 #endif
523 
524  if (GEOSisValid(gunion))
525  break;
526 
527  /* make geometry valid */
528  mpolyValid = lwgeom_make_valid(mpoly);
529  if (mpolyValid == NULL) {
530  rtwarn("Cannot fix invalid geometry");
531  break;
532  }
533 
534  lwgeom_free(mpoly);
535  mpoly = mpolyValid;
536  }
537  while (0);
538 
539  GEOSGeom_destroy(gunion);
540  }
541  else {
542  mpoly = lwpoly_as_lwgeom(gv[0].geom);
543  rtdealloc(gv);
544 
545 #if POSTGIS_DEBUG_LEVEL > 3
546  {
547  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
548  RASTER_DEBUGF(4, "geom 0 = %s", wkt);
549  rtdealloc(wkt);
550  }
551 #endif
552  }
553 
554  /* specify SRID */
555  lwgeom_set_srid(mpoly, rt_raster_get_srid(raster));
556 
557  if (mpoly != NULL) {
558  /* convert to multi */
559  if (!lwgeom_is_collection(mpoly)) {
560  tmp = mpoly;
561 
562 #if POSTGIS_DEBUG_LEVEL > 3
563  {
564  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
565  RASTER_DEBUGF(4, "before multi = %s", wkt);
566  rtdealloc(wkt);
567  }
568 #endif
569 
570  RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
571 
572  /*
573  lwgeom_as_multi() only does a shallow clone internally
574  so input and output geometries may share memory
575  hence the deep clone of the output geometry for returning
576  is the only way to guarentee the memory isn't shared
577  */
578  mpoly = lwgeom_as_multi(tmp);
579  clone = lwgeom_clone_deep(mpoly);
580  lwgeom_free(tmp);
581  lwgeom_free(mpoly);
582  mpoly = clone;
583 
584  RASTER_DEBUGF(4, "mpoly @ %p", mpoly);
585 
586 #if POSTGIS_DEBUG_LEVEL > 3
587  {
588  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
589  RASTER_DEBUGF(4, "after multi = %s", wkt);
590  rtdealloc(wkt);
591  }
592 #endif
593  }
594 
595 #if POSTGIS_DEBUG_LEVEL > 3
596  {
597  char *wkt = lwgeom_to_wkt(mpoly, WKT_ISO, DBL_DIG, NULL);
598  RASTER_DEBUGF(4, "returning geometry = %s", wkt);
599  rtdealloc(wkt);
600  }
601 #endif
602 
603  *surface = lwgeom_as_lwmpoly(mpoly);
604  return ES_NONE;
605  }
606 
607  return ES_NONE;
608 }
609 
610 /******************************************************************************
611 * rt_raster_pixel_as_polygon()
612 ******************************************************************************/
613 
628 LWPOLY*
630 {
631  double scale_x, scale_y;
632  double skew_x, skew_y;
633  double ul_x, ul_y;
634  int srid;
635  POINTARRAY **points;
636  POINT4D p, p0;
637  LWPOLY *poly;
638 
639  assert(rast != NULL);
640 
641  scale_x = rt_raster_get_x_scale(rast);
642  scale_y = rt_raster_get_y_scale(rast);
643  skew_x = rt_raster_get_x_skew(rast);
644  skew_y = rt_raster_get_y_skew(rast);
645  ul_x = rt_raster_get_x_offset(rast);
646  ul_y = rt_raster_get_y_offset(rast);
647  srid = rt_raster_get_srid(rast);
648 
649  points = rtalloc(sizeof(POINTARRAY *)*1);
650  points[0] = ptarray_construct(0, 0, 5);
651 
652  p0.x = scale_x * x + skew_x * y + ul_x;
653  p0.y = scale_y * y + skew_y * x + ul_y;
654  ptarray_set_point4d(points[0], 0, &p0);
655 
656  p.x = p0.x + scale_x;
657  p.y = p0.y + skew_y;
658  ptarray_set_point4d(points[0], 1, &p);
659 
660  p.x = p0.x + scale_x + skew_x;
661  p.y = p0.y + scale_y + skew_y;
662  ptarray_set_point4d(points[0], 2, &p);
663 
664  p.x = p0.x + skew_x;
665  p.y = p0.y + scale_y;
666  ptarray_set_point4d(points[0], 3, &p);
667 
668  /* close it */
669  ptarray_set_point4d(points[0], 4, &p0);
670 
671  poly = lwpoly_construct(srid, NULL, 1, points);
672 
673  return poly;
674 }
675 
676 /******************************************************************************
677 * rt_raster_get_envelope_geom()
678 ******************************************************************************/
679 
690  double gt[6] = {0.0};
691  int srid = SRID_UNKNOWN;
692 
693  POINTARRAY *pts = NULL;
694  POINT4D p4d;
695 
696  assert(env != NULL);
697  *env = NULL;
698 
699  /* raster is NULL, envelope is NULL */
700  if (raster == NULL)
701  return ES_NONE;
702 
703  /* raster metadata */
704  srid = rt_raster_get_srid(raster);
706 
708  3,
709  "rt_raster_get_envelope: raster is %dx%d",
710  raster->width,
711  raster->height
712  );
713 
714  /* return point or line since at least one of the two dimensions is 0 */
715  if ((!raster->width) || (!raster->height)) {
716  p4d.x = gt[0];
717  p4d.y = gt[3];
718 
719  /* return point */
720  if (!raster->width && !raster->height) {
721  LWPOINT *point = lwpoint_make2d(srid, p4d.x, p4d.y);
722  *env = lwpoint_as_lwgeom(point);
723  }
724  /* return linestring */
725  else {
726  LWLINE *line = NULL;
727  pts = ptarray_construct_empty(0, 0, 2);
728 
729  /* first point of line */
730  ptarray_append_point(pts, &p4d, LW_TRUE);
731 
732  /* second point of line */
734  raster,
735  rt_raster_get_width(raster), rt_raster_get_height(raster),
736  &p4d.x, &p4d.y,
737  gt
738  ) != ES_NONE) {
739  rterror("rt_raster_get_envelope: Could not get second point for linestring");
740  return ES_ERROR;
741  }
742  ptarray_append_point(pts, &p4d, LW_TRUE);
743  line = lwline_construct(srid, NULL, pts);
744 
745  *env = lwline_as_lwgeom(line);
746  }
747 
748  return ES_NONE;
749  }
750  else {
751  rt_envelope rtenv;
752  int err = ES_NONE;
753  POINTARRAY **rings = NULL;
754  LWPOLY* poly = NULL;
755 
756  /* only one ring */
757  rings = (POINTARRAY **) rtalloc(sizeof (POINTARRAY*));
758  if (!rings) {
759  rterror("rt_raster_get_envelope_geom: Could not allocate memory for polygon ring");
760  return ES_ERROR;
761  }
762  rings[0] = ptarray_construct(0, 0, 5);
763  if (!rings[0]) {
764  rterror("rt_raster_get_envelope_geom: Could not construct point array");
765  return ES_ERROR;
766  }
767  pts = rings[0];
768 
769  err = rt_raster_get_envelope(raster, &rtenv);
770  if (err != ES_NONE) {
771  rterror("rt_raster_get_envelope_geom: Could not get raster envelope");
772  return err;
773  }
774 
775  /* build ring */
776 
777  /* minx, maxy */
778  p4d.x = rtenv.MinX;
779  p4d.y = rtenv.MaxY;
780  ptarray_set_point4d(pts, 0, &p4d);
781  ptarray_set_point4d(pts, 4, &p4d);
782 
783  /* maxx, maxy */
784  p4d.x = rtenv.MaxX;
785  p4d.y = rtenv.MaxY;
786  ptarray_set_point4d(pts, 1, &p4d);
787 
788  /* maxx, miny */
789  p4d.x = rtenv.MaxX;
790  p4d.y = rtenv.MinY;
791  ptarray_set_point4d(pts, 2, &p4d);
792 
793  /* minx, miny */
794  p4d.x = rtenv.MinX;
795  p4d.y = rtenv.MinY;
796  ptarray_set_point4d(pts, 3, &p4d);
797 
798  poly = lwpoly_construct(srid, 0, 1, rings);
799  *env = lwpoly_as_lwgeom(poly);
800  }
801 
802  return ES_NONE;
803 }
804 
805 /******************************************************************************
806 * rt_raster_get_convex_hull()
807 ******************************************************************************/
808 
823  double gt[6] = {0.0};
824  int srid = SRID_UNKNOWN;
825 
826  POINTARRAY *pts = NULL;
827  POINT4D p4d;
828 
829  assert(hull != NULL);
830  *hull = NULL;
831 
832  /* raster is NULL, convex hull is NULL */
833  if (raster == NULL)
834  return ES_NONE;
835 
836  /* raster metadata */
837  srid = rt_raster_get_srid(raster);
839 
840  RASTER_DEBUGF(3, "rt_raster_get_convex_hull: raster is %dx%d", raster->width, raster->height);
841 
842  /* return point or line since at least one of the two dimensions is 0 */
843  if ((!raster->width) || (!raster->height)) {
844  p4d.x = gt[0];
845  p4d.y = gt[3];
846 
847  /* return point */
848  if (!raster->width && !raster->height) {
849  LWPOINT *point = lwpoint_make2d(srid, p4d.x, p4d.y);
850  *hull = lwpoint_as_lwgeom(point);
851  }
852  /* return linestring */
853  else {
854  LWLINE *line = NULL;
855  pts = ptarray_construct_empty(0, 0, 2);
856 
857  /* first point of line */
858  ptarray_append_point(pts, &p4d, LW_TRUE);
859 
860  /* second point of line */
862  raster,
863  rt_raster_get_width(raster), rt_raster_get_height(raster),
864  &p4d.x, &p4d.y,
865  gt
866  ) != ES_NONE) {
867  rterror("rt_raster_get_convex_hull: Could not get second point for linestring");
868  return ES_ERROR;
869  }
870  ptarray_append_point(pts, &p4d, LW_TRUE);
871  line = lwline_construct(srid, NULL, pts);
872 
873  *hull = lwline_as_lwgeom(line);
874  }
875 
876  return ES_NONE;
877  }
878  else {
879  POINTARRAY **rings = NULL;
880  LWPOLY* poly = NULL;
881 
882  /* only one ring */
883  rings = (POINTARRAY **) rtalloc(sizeof (POINTARRAY*));
884  if (!rings) {
885  rterror("rt_raster_get_convex_hull: Could not allocate memory for polygon ring");
886  return ES_ERROR;
887  }
888  rings[0] = ptarray_construct(0, 0, 5);
889  /* TODO: handle error on ptarray construction */
890  /* XXX jorgearevalo: the error conditions aren't managed in ptarray_construct */
891  if (!rings[0]) {
892  rterror("rt_raster_get_convex_hull: Could not construct point array");
893  return ES_ERROR;
894  }
895  pts = rings[0];
896 
897  /* Upper-left corner (first and last points) */
898  p4d.x = gt[0];
899  p4d.y = gt[3];
900  ptarray_set_point4d(pts, 0, &p4d);
901  ptarray_set_point4d(pts, 4, &p4d);
902 
903  /* Upper-right corner (we go clockwise) */
905  raster,
906  raster->width, 0,
907  &p4d.x, &p4d.y,
908  gt
909  );
910  ptarray_set_point4d(pts, 1, &p4d);
911 
912  /* Lower-right corner */
914  raster,
915  raster->width, raster->height,
916  &p4d.x, &p4d.y,
917  gt
918  );
919  ptarray_set_point4d(pts, 2, &p4d);
920 
921  /* Lower-left corner */
923  raster,
924  0, raster->height,
925  &p4d.x, &p4d.y,
926  gt
927  );
928  ptarray_set_point4d(pts, 3, &p4d);
929 
930  poly = lwpoly_construct(srid, 0, 1, rings);
931  *hull = lwpoly_as_lwgeom(poly);
932  }
933 
934  return ES_NONE;
935 }
936 
937 /******************************************************************************
938 * rt_raster_gdal_polygonize()
939 ******************************************************************************/
940 
960  rt_raster raster, int nband,
961  int exclude_nodata_value,
962  int *pnElements
963 ) {
964  CPLErr cplerr = CE_None;
965  char *pszQuery;
966  long j;
967  OGRSFDriverH ogr_drv = NULL;
968  GDALDriverH gdal_drv = NULL;
969  int destroy_gdal_drv = 0;
970  GDALDatasetH memdataset = NULL;
971  GDALRasterBandH gdal_band = NULL;
972  OGRDataSourceH memdatasource = NULL;
973  rt_geomval pols = NULL;
974  OGRLayerH hLayer = NULL;
975  OGRFeatureH hFeature = NULL;
976  OGRGeometryH hGeom = NULL;
977  OGRFieldDefnH hFldDfn = NULL;
978  unsigned char *wkb = NULL;
979  int wkbsize = 0;
980  LWGEOM *lwgeom = NULL;
981  int nFeatureCount = 0;
982  rt_band band = NULL;
983  int iPixVal = -1;
984  double dValue = 0.0;
985  int iBandHasNodataValue = FALSE;
986  double dBandNoData = 0.0;
987 
988  /* for checking that a geometry is valid */
989  GEOSGeometry *ggeom = NULL;
990  int isValid;
991  LWGEOM *lwgeomValid = NULL;
992 
993  uint32_t bandNums[1] = {nband};
994  int excludeNodataValues[1] = {exclude_nodata_value};
995 
996  /* checks */
997  assert(NULL != raster);
998  assert(NULL != pnElements);
999 
1000  RASTER_DEBUG(2, "In rt_raster_gdal_polygonize");
1001 
1002  *pnElements = 0;
1003 
1004  /*******************************
1005  * Get band
1006  *******************************/
1007  band = rt_raster_get_band(raster, nband);
1008  if (NULL == band) {
1009  rterror("rt_raster_gdal_polygonize: Error getting band %d from raster", nband);
1010  return NULL;
1011  }
1012 
1013  if (exclude_nodata_value) {
1014 
1015  /* band is NODATA */
1016  if (rt_band_get_isnodata_flag(band)) {
1017  RASTER_DEBUG(3, "Band is NODATA. Returning null");
1018  *pnElements = 0;
1019  return NULL;
1020  }
1021 
1022  iBandHasNodataValue = rt_band_get_hasnodata_flag(band);
1023  if (iBandHasNodataValue)
1024  rt_band_get_nodata(band, &dBandNoData);
1025  else
1026  exclude_nodata_value = FALSE;
1027  }
1028 
1029  /*****************************************************
1030  * Convert raster to GDAL MEM dataset
1031  *****************************************************/
1032  memdataset = rt_raster_to_gdal_mem(raster, NULL, bandNums, excludeNodataValues, 1, &gdal_drv, &destroy_gdal_drv);
1033  if (NULL == memdataset) {
1034  rterror("rt_raster_gdal_polygonize: Couldn't convert raster to GDAL MEM dataset");
1035  return NULL;
1036  }
1037 
1038  /*****************************
1039  * Register ogr mem driver
1040  *****************************/
1041 #ifdef GDAL_DCAP_RASTER
1042  /* in GDAL 2.0, OGRRegisterAll() is an alias to GDALAllRegister() */
1044 #else
1045  OGRRegisterAll();
1046 #endif
1047 
1048  RASTER_DEBUG(3, "creating OGR MEM vector");
1049 
1050  /*****************************************************
1051  * Create an OGR in-memory vector for layers
1052  *****************************************************/
1053  ogr_drv = OGRGetDriverByName("Memory");
1054  memdatasource = OGR_Dr_CreateDataSource(ogr_drv, "", NULL);
1055  if (NULL == memdatasource) {
1056  rterror("rt_raster_gdal_polygonize: Couldn't create a OGR Datasource to store pols");
1057  GDALClose(memdataset);
1058  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1059  return NULL;
1060  }
1061 
1062  /* Can MEM driver create new layers? */
1063  if (!OGR_DS_TestCapability(memdatasource, ODsCCreateLayer)) {
1064  rterror("rt_raster_gdal_polygonize: MEM driver can't create new layers, aborting");
1065 
1066  /* xxx jorgearevalo: what should we do now? */
1067  GDALClose(memdataset);
1068  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1069  OGRReleaseDataSource(memdatasource);
1070 
1071  return NULL;
1072  }
1073 
1074  RASTER_DEBUG(3, "polygonizying GDAL MEM raster band");
1075 
1076  /*****************************
1077  * Polygonize the raster band
1078  *****************************/
1079 
1085  hLayer = OGR_DS_CreateLayer(memdatasource, "PolygonizedLayer", NULL, wkbPolygon, NULL);
1086 
1087  if (NULL == hLayer) {
1088  rterror("rt_raster_gdal_polygonize: Couldn't create layer to store polygons");
1089 
1090  GDALClose(memdataset);
1091  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1092  OGRReleaseDataSource(memdatasource);
1093 
1094  return NULL;
1095  }
1096 
1101  /* First, create a field definition to create the field */
1102  hFldDfn = OGR_Fld_Create("PixelValue", OFTReal);
1103 
1104  /* Second, create the field */
1105  if (OGR_L_CreateField(hLayer, hFldDfn, TRUE) != OGRERR_NONE) {
1106  rtwarn("Couldn't create a field in OGR Layer. The polygons generated won't be able to store the pixel value");
1107  iPixVal = -1;
1108  }
1109  else {
1110  /* Index to the new field created in the layer */
1111  iPixVal = 0;
1112  }
1113 
1114  /* Get GDAL raster band */
1115  gdal_band = GDALGetRasterBand(memdataset, 1);
1116  if (NULL == gdal_band) {
1117  rterror("rt_raster_gdal_polygonize: Couldn't get GDAL band to polygonize");
1118 
1119  GDALClose(memdataset);
1120  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1121  OGR_Fld_Destroy(hFldDfn);
1122  OGR_DS_DeleteLayer(memdatasource, 0);
1123  OGRReleaseDataSource(memdatasource);
1124 
1125  return NULL;
1126  }
1127 
1131 #ifdef GDALFPOLYGONIZE
1132  cplerr = GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1133 #else
1134  cplerr = GDALPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1135 #endif
1136 
1137  if (cplerr != CE_None) {
1138  rterror("rt_raster_gdal_polygonize: Could not polygonize GDAL band");
1139 
1140  GDALClose(memdataset);
1141  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1142  OGR_Fld_Destroy(hFldDfn);
1143  OGR_DS_DeleteLayer(memdatasource, 0);
1144  OGRReleaseDataSource(memdatasource);
1145 
1146  return NULL;
1147  }
1148 
1155  if (iBandHasNodataValue) {
1156  pszQuery = (char *) rtalloc(50 * sizeof (char));
1157  sprintf(pszQuery, "PixelValue != %f", dBandNoData );
1158  OGRErr e = OGR_L_SetAttributeFilter(hLayer, pszQuery);
1159  if (e != OGRERR_NONE) {
1160  rtwarn("Error filtering NODATA values for band. All values will be treated as data values");
1161  }
1162  }
1163  else {
1164  pszQuery = NULL;
1165  }
1166 
1167  /*********************************************************************
1168  * Transform OGR layers to WKB polygons
1169  * XXX jorgearevalo: GDALPolygonize does not set the coordinate system
1170  * on the output layer. Application code should do this when the layer
1171  * is created, presumably matching the raster coordinate system.
1172  *********************************************************************/
1173  nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);
1174 
1175  /* Allocate memory for pols */
1176  pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t));
1177 
1178  if (NULL == pols) {
1179  rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set");
1180 
1181  GDALClose(memdataset);
1182  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1183  OGR_Fld_Destroy(hFldDfn);
1184  OGR_DS_DeleteLayer(memdatasource, 0);
1185  if (NULL != pszQuery)
1186  rtdealloc(pszQuery);
1187  OGRReleaseDataSource(memdatasource);
1188 
1189  return NULL;
1190  }
1191 
1192  /* initialize GEOS */
1193  initGEOS(rtinfo, lwgeom_geos_error);
1194 
1195  RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount);
1196 
1197  /* Reset feature reading to start in the first feature */
1198  OGR_L_ResetReading(hLayer);
1199 
1200  for (j = 0; j < nFeatureCount; j++) {
1201  hFeature = OGR_L_GetNextFeature(hLayer);
1202  dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
1203 
1204  hGeom = OGR_F_GetGeometryRef(hFeature);
1205  wkbsize = OGR_G_WkbSize(hGeom);
1206 
1207  /* allocate wkb buffer */
1208  wkb = rtalloc(sizeof(unsigned char) * wkbsize);
1209  if (wkb == NULL) {
1210  rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer");
1211 
1212  OGR_F_Destroy(hFeature);
1213  GDALClose(memdataset);
1214  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1215  OGR_Fld_Destroy(hFldDfn);
1216  OGR_DS_DeleteLayer(memdatasource, 0);
1217  if (NULL != pszQuery)
1218  rtdealloc(pszQuery);
1219  OGRReleaseDataSource(memdatasource);
1220 
1221  return NULL;
1222  }
1223 
1224  /* export WKB with LSB byte order */
1225  OGR_G_ExportToWkb(hGeom, wkbNDR, wkb);
1226 
1227  /* convert WKB to LWGEOM */
1228  lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE);
1229 
1230 #if POSTGIS_DEBUG_LEVEL > 3
1231  {
1232  char *wkt = NULL;
1233  OGR_G_ExportToWkt(hGeom, &wkt);
1234  RASTER_DEBUGF(4, "GDAL wkt = %s", wkt);
1235  CPLFree(wkt);
1236 
1237  d_print_binary_hex("GDAL wkb", wkb, wkbsize);
1238  }
1239 #endif
1240 
1241  /* cleanup unnecessary stuff */
1242  rtdealloc(wkb);
1243  wkb = NULL;
1244  wkbsize = 0;
1245 
1246  OGR_F_Destroy(hFeature);
1247 
1248  /* specify SRID */
1249  lwgeom_set_srid(lwgeom, rt_raster_get_srid(raster));
1250 
1251  /*
1252  is geometry valid?
1253  if not, try to make valid
1254  */
1255  do {
1256 #if POSTGIS_GEOS_VERSION < 33
1257  /* nothing can be done if the geometry was invalid if GEOS < 3.3 */
1258  break;
1259 #endif
1260 
1261  ggeom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
1262  if (ggeom == NULL) {
1263  rtwarn("Cannot test geometry for validity");
1264  break;
1265  }
1266 
1267  isValid = GEOSisValid(ggeom);
1268 
1269  GEOSGeom_destroy(ggeom);
1270  ggeom = NULL;
1271 
1272  /* geometry is valid */
1273  if (isValid)
1274  break;
1275 
1276  RASTER_DEBUG(3, "fixing invalid geometry");
1277 
1278  /* make geometry valid */
1279  lwgeomValid = lwgeom_make_valid(lwgeom);
1280  if (lwgeomValid == NULL) {
1281  rtwarn("Cannot fix invalid geometry");
1282  break;
1283  }
1284 
1285  lwgeom_free(lwgeom);
1286  lwgeom = lwgeomValid;
1287  }
1288  while (0);
1289 
1290  /* save lwgeom */
1291  pols[j].geom = lwgeom_as_lwpoly(lwgeom);
1292 
1293 #if POSTGIS_DEBUG_LEVEL > 3
1294  {
1295  char *wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, NULL);
1296  RASTER_DEBUGF(4, "LWGEOM wkt = %s", wkt);
1297  rtdealloc(wkt);
1298 
1299  size_t lwwkbsize = 0;
1300  uint8_t *lwwkb = lwgeom_to_wkb(lwgeom, WKB_ISO | WKB_NDR, &lwwkbsize);
1301  if (lwwkbsize) {
1302  d_print_binary_hex("LWGEOM wkb", lwwkb, lwwkbsize);
1303  rtdealloc(lwwkb);
1304  }
1305  }
1306 #endif
1307 
1308  /* set pixel value */
1309  pols[j].val = dValue;
1310  }
1311 
1312  *pnElements = nFeatureCount;
1313 
1314  RASTER_DEBUG(3, "destroying GDAL MEM raster");
1315  GDALClose(memdataset);
1316  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1317 
1318  RASTER_DEBUG(3, "destroying OGR MEM vector");
1319  OGR_Fld_Destroy(hFldDfn);
1320  OGR_DS_DeleteLayer(memdatasource, 0);
1321  if (NULL != pszQuery) rtdealloc(pszQuery);
1322  OGRReleaseDataSource(memdatasource);
1323 
1324  return pols;
1325 }
1326 
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:536
double x
Definition: liblwgeom.h:336
double MinY
Definition: librtcore.h:179
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
int rt_util_gdal_register_all(int force_register_all)
Definition: rt_util.c:334
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition: rt_raster.c:213
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
#define WKB_NDR
Definition: liblwgeom.h:1933
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
Definition: rt_raster.c:190
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:991
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:655
double MaxY
Definition: librtcore.h:180
rt_raster raster
Definition: librtcore.h:2287
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:706
double MaxX
Definition: librtcore.h:178
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:132
static rt_errorstate _rti_raster_get_band_perimeter(rt_band band, uint16_t *trim)
Definition: rt_geometry.c:39
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1050
tuple gt
Definition: window.py:77
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:1809
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
tuple band
Definition: ovdump.py:57
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
tuple rast
Definition: rtpixdump.py:61
rt_errorstate
Enum definitions.
Definition: librtcore.h:191
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:433
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:125
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition: lwgeom.c:284
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:239
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition: rt_raster.c:1338
uint16_t height
Definition: librtcore.h:2265
LWMPOLY * lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:170
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
LWPOLY * geom
Definition: librtcore.h:2316
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:1869
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
LWPOLY * rt_raster_pixel_as_polygon(rt_raster rast, int x, int y)
Get a raster pixel as a polygon.
Definition: rt_geometry.c:629
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition: rt_geometry.c:822
double val
Definition: librtcore.h:2317
#define WKT_ISO
Definition: liblwgeom.h:1939
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:750
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:249
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_TRUE, then a duplicate point will not be added.
Definition: ptarray.c:156
LWGEOM * lwgeom_make_valid(LWGEOM *geom)
Attempts to make an invalid geometries valid w/out losing points.
void lwgeom_geos_error(const char *fmt,...)
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
LWPOLY * lwpoly_construct(int srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points)
Definition: lwpoly.c:29
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:79
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1088
tuple nband
Definition: pixval.py:52
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:61
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
uint16_t width
Definition: librtcore.h:2264
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:172
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:959
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
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:541
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
Definition: rt_raster.c:150
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_band.c:507
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_band.c:516
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:311
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:356
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
Definition: rt_raster.c:159
#define WKB_ISO
Definition: liblwgeom.h:1930
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_raster.c:181
tuple x
Definition: pixval.py:53
rt_errorstate rt_raster_get_envelope_geom(rt_raster raster, LWGEOM **env)
Get raster's envelope as a geometry.
Definition: rt_geometry.c:689
struct rt_geomval_t * rt_geomval
Definition: librtcore.h:161
rt_message_handler err
Definition: rt_context.c:112
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
void rtdealloc(void *mem)
Definition: rt_context.c:186
#define FALSE
Definition: dbfopen.c:168
double MinX
Definition: librtcore.h:177
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:307
This library is the generic raster handling section of PostGIS.
rt_errorstate rt_raster_get_perimeter(rt_raster raster, int nband, LWGEOM **perimeter)
Get raster perimeter.
Definition: rt_geometry.c:182
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
static char * trim(const char *input)
Definition: raster2pgsql.c:262
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:254
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition: rt_band.c:581
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...
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
double y
Definition: liblwgeom.h:336
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:754
rt_errorstate rt_raster_surface(rt_raster raster, int nband, LWMPOLY **surface)
Get a raster as a surface (multipolygon).
Definition: rt_geometry.c:355
#define TRUE
Definition: dbfopen.c:169
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_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition: rt_raster.c:222
tuple y
Definition: pixval.py:54