PostGIS  2.3.8dev-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 */
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  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 */
536  lwgeom_set_srid(mpoly, rt_raster_get_srid(raster));
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);
626  ul_x = rt_raster_get_x_offset(rast);
627  ul_y = rt_raster_get_y_offset(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,
716  rt_raster_get_width(raster), rt_raster_get_height(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,
844  rt_raster_get_width(raster), rt_raster_get_height(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  *******************************/
988  band = rt_raster_get_band(raster, nband);
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 */
997  if (rt_band_get_isnodata_flag(band)) {
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 
1112 #ifdef GDALFPOLYGONIZE
1113  cplerr = GDALFPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1114 #else
1115  cplerr = GDALPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1116 #endif
1117 
1118  if (cplerr != CE_None) {
1119  rterror("rt_raster_gdal_polygonize: Could not polygonize GDAL band");
1120 
1121  GDALClose(memdataset);
1122  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1123  OGR_Fld_Destroy(hFldDfn);
1124  OGR_DS_DeleteLayer(memdatasource, 0);
1125  OGRReleaseDataSource(memdatasource);
1126 
1127  return NULL;
1128  }
1129 
1136  if (iBandHasNodataValue) {
1137  pszQuery = (char *) rtalloc(50 * sizeof (char));
1138  sprintf(pszQuery, "PixelValue != %f", dBandNoData );
1139  OGRErr e = OGR_L_SetAttributeFilter(hLayer, pszQuery);
1140  if (e != OGRERR_NONE) {
1141  rtwarn("Error filtering NODATA values for band. All values will be treated as data values");
1142  }
1143  }
1144  else {
1145  pszQuery = NULL;
1146  }
1147 
1148  /*********************************************************************
1149  * Transform OGR layers to WKB polygons
1150  * XXX jorgearevalo: GDALPolygonize does not set the coordinate system
1151  * on the output layer. Application code should do this when the layer
1152  * is created, presumably matching the raster coordinate system.
1153  *********************************************************************/
1154  nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);
1155 
1156  /* Allocate memory for pols */
1157  pols = (rt_geomval) rtalloc(nFeatureCount * sizeof(struct rt_geomval_t));
1158 
1159  if (NULL == pols) {
1160  rterror("rt_raster_gdal_polygonize: Could not allocate memory for geomval set");
1161 
1162  GDALClose(memdataset);
1163  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1164  OGR_Fld_Destroy(hFldDfn);
1165  OGR_DS_DeleteLayer(memdatasource, 0);
1166  if (NULL != pszQuery)
1167  rtdealloc(pszQuery);
1168  OGRReleaseDataSource(memdatasource);
1169 
1170  return NULL;
1171  }
1172 
1173  /* initialize GEOS */
1174  initGEOS(rtinfo, lwgeom_geos_error);
1175 
1176  RASTER_DEBUGF(3, "storing polygons (%d)", nFeatureCount);
1177 
1178  /* Reset feature reading to start in the first feature */
1179  OGR_L_ResetReading(hLayer);
1180 
1181  for (j = 0; j < nFeatureCount; j++) {
1182  hFeature = OGR_L_GetNextFeature(hLayer);
1183  dValue = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
1184 
1185  hGeom = OGR_F_GetGeometryRef(hFeature);
1186  wkbsize = OGR_G_WkbSize(hGeom);
1187 
1188  /* allocate wkb buffer */
1189  wkb = rtalloc(sizeof(unsigned char) * wkbsize);
1190  if (wkb == NULL) {
1191  rterror("rt_raster_gdal_polygonize: Could not allocate memory for WKB buffer");
1192 
1193  OGR_F_Destroy(hFeature);
1194  GDALClose(memdataset);
1195  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1196  OGR_Fld_Destroy(hFldDfn);
1197  OGR_DS_DeleteLayer(memdatasource, 0);
1198  if (NULL != pszQuery)
1199  rtdealloc(pszQuery);
1200  OGRReleaseDataSource(memdatasource);
1201 
1202  return NULL;
1203  }
1204 
1205  /* export WKB with LSB byte order */
1206  OGR_G_ExportToWkb(hGeom, wkbNDR, wkb);
1207 
1208  /* convert WKB to LWGEOM */
1209  lwgeom = lwgeom_from_wkb(wkb, wkbsize, LW_PARSER_CHECK_NONE);
1210 
1211 #if POSTGIS_DEBUG_LEVEL > 3
1212  {
1213  char *wkt = NULL;
1214  OGR_G_ExportToWkt(hGeom, &wkt);
1215  RASTER_DEBUGF(4, "GDAL wkt = %s", wkt);
1216  CPLFree(wkt);
1217 
1218  d_print_binary_hex("GDAL wkb", wkb, wkbsize);
1219  }
1220 #endif
1221 
1222  /* cleanup unnecessary stuff */
1223  rtdealloc(wkb);
1224  wkb = NULL;
1225  wkbsize = 0;
1226 
1227  OGR_F_Destroy(hFeature);
1228 
1229  /* specify SRID */
1230  lwgeom_set_srid(lwgeom, rt_raster_get_srid(raster));
1231 
1232  /*
1233  is geometry valid?
1234  if not, try to make valid
1235  */
1236  do {
1237  ggeom = (GEOSGeometry *) LWGEOM2GEOS(lwgeom, 0);
1238  if (ggeom == NULL) {
1239  rtwarn("Cannot test geometry for validity");
1240  break;
1241  }
1242 
1243  isValid = GEOSisValid(ggeom);
1244 
1245  GEOSGeom_destroy(ggeom);
1246  ggeom = NULL;
1247 
1248  /* geometry is valid */
1249  if (isValid)
1250  break;
1251 
1252  RASTER_DEBUG(3, "fixing invalid geometry");
1253 
1254  /* make geometry valid */
1255  lwgeomValid = lwgeom_make_valid(lwgeom);
1256  if (lwgeomValid == NULL) {
1257  rtwarn("Cannot fix invalid geometry");
1258  break;
1259  }
1260 
1261  lwgeom_free(lwgeom);
1262  lwgeom = lwgeomValid;
1263  }
1264  while (0);
1265 
1266  /* save lwgeom */
1267  pols[j].geom = lwgeom_as_lwpoly(lwgeom);
1268 
1269 #if POSTGIS_DEBUG_LEVEL > 3
1270  {
1271  char *wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, NULL);
1272  RASTER_DEBUGF(4, "LWGEOM wkt = %s", wkt);
1273  rtdealloc(wkt);
1274 
1275  size_t lwwkbsize = 0;
1276  uint8_t *lwwkb = lwgeom_to_wkb(lwgeom, WKB_ISO | WKB_NDR, &lwwkbsize);
1277  if (lwwkbsize) {
1278  d_print_binary_hex("LWGEOM wkb", lwwkb, lwwkbsize);
1279  rtdealloc(lwwkb);
1280  }
1281  }
1282 #endif
1283 
1284  /* set pixel value */
1285  pols[j].val = dValue;
1286  }
1287 
1288  *pnElements = nFeatureCount;
1289 
1290  RASTER_DEBUG(3, "destroying GDAL MEM raster");
1291  GDALClose(memdataset);
1292  if (destroy_gdal_drv) GDALDestroyDriver(gdal_drv);
1293 
1294  RASTER_DEBUG(3, "destroying OGR MEM vector");
1295  OGR_Fld_Destroy(hFldDfn);
1296  OGR_DS_DeleteLayer(memdatasource, 0);
1297  if (NULL != pszQuery) rtdealloc(pszQuery);
1298  OGRReleaseDataSource(memdatasource);
1299 
1300  return pols;
1301 }
1302 
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:549
double x
Definition: liblwgeom.h:351
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:2049
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:1004
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:669
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
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:145
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:1063
band
Definition: ovdump.py:57
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
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
rt_errorstate
Enum definitions.
Definition: librtcore.h:191
gt
Definition: window.py:77
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:446
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:138
LWGEOM * lwgeom_as_multi(const LWGEOM *lwgeom)
Create a new LWGEOM of the appropriate MULTI* type.
Definition: lwgeom.c:297
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:252
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:183
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:1985
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
nband
Definition: pixval.py:52
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&#39;s convex hull.
Definition: rt_geometry.c:803
double val
Definition: librtcore.h:2317
#define WKT_ISO
Definition: liblwgeom.h:2055
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
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:262
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, 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:43
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:152
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
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
uint16_t width
Definition: librtcore.h:2264
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:187
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
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&#39;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:2046
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
Definition: rt_raster.c:181
rt_errorstate rt_raster_get_envelope_geom(rt_raster raster, LWGEOM **env)
Get raster&#39;s envelope as a geometry.
Definition: rt_geometry.c:670
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:267
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...
int value
Definition: genraster.py:61
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
double y
Definition: liblwgeom.h:351
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:765
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&#39;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