PostGIS  2.3.7dev-r@@SVN_REVISION@@
cu_gdal.c
Go to the documentation of this file.
1 /*
2  * PostGIS Raster - Raster Types for PostGIS
3  * http://trac.osgeo.org/postgis/wiki/WKTRaster
4  *
5  * Copyright (C) 2012 Regents of the University of California
6  * <bkpark@ucdavis.edu>
7  *
8  * This program is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License
10  * as published by the Free Software Foundation; either version 2
11  * of the License, or (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software Foundation,
20  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
21  *
22  */
23 
24 #include "CUnit/Basic.h"
25 #include "cu_tester.h"
26 
27 static void test_gdal_configured() {
28  CU_ASSERT(rt_util_gdal_configured());
29 }
30 
31 static void test_gdal_drivers() {
32  int i;
33  uint32_t size;
34  rt_gdaldriver drv = NULL;
35 
36  drv = (rt_gdaldriver) rt_raster_gdal_drivers(&size, 1);
37  CU_ASSERT(drv != NULL);
38 
39  for (i = 0; i < size; i++) {
40  CU_ASSERT(strlen(drv[i].short_name));
41  rtdealloc(drv[i].short_name);
42  rtdealloc(drv[i].long_name);
43  rtdealloc(drv[i].create_options);
44  }
45 
46  rtdealloc(drv);
47 }
48 
49 static void test_gdal_rasterize() {
51  char srs[] = "PROJCS[\"unnamed\",GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6370997,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center\",45],PARAMETER[\"longitude_of_center\",-100],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1],AUTHORITY[\"EPSG\",\"2163\"]]";
52  const char wkb_hex[] = "010300000001000000050000000000000080841ec100000000600122410000000080841ec100000000804f22410000000040e81dc100000000804f22410000000040e81dc100000000600122410000000080841ec10000000060012241";
53  const char *pos = wkb_hex;
54  unsigned char *wkb = NULL;
55  int wkb_len = 0;
56  int i;
57  double scale_x = 100;
58  double scale_y = -100;
59 
60  rt_pixtype pixtype[] = {PT_8BUI};
61  double init[] = {0};
62  double value[] = {1};
63  double nodata[] = {0};
64  uint8_t nodata_mask[] = {1};
65 
66  /* hex to byte */
67  wkb_len = (int) ceil(((double) strlen(wkb_hex)) / 2);
68  wkb = (unsigned char *) rtalloc(sizeof(unsigned char) * wkb_len);
69  for (i = 0; i < wkb_len; i++) {
70  sscanf(pos, "%2hhx", &wkb[i]);
71  pos += 2;
72  }
73 
74  raster = rt_raster_gdal_rasterize(
75  wkb,
76  wkb_len, srs,
77  1, pixtype,
78  init, value,
79  nodata, nodata_mask,
80  NULL, NULL,
81  &scale_x, &scale_y,
82  NULL, NULL,
83  NULL, NULL,
84  NULL, NULL,
85  NULL
86  );
87 
88  CU_ASSERT(raster != NULL);
89  CU_ASSERT_EQUAL(rt_raster_get_width(raster), 100);
90  CU_ASSERT_EQUAL(rt_raster_get_height(raster), 100);
91  CU_ASSERT_NOT_EQUAL(rt_raster_get_num_bands(raster), 0);
92  CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_x_offset(raster), -500000, DBL_EPSILON);
93  CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_y_offset(raster), 600000, DBL_EPSILON);
94 
95  rtdealloc(wkb);
96  cu_free_raster(raster);
97 }
98 
99 static char *
100 lwgeom_to_text(const LWGEOM *lwgeom) {
101  char *wkt;
102  size_t wkt_size;
103 
104  wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, &wkt_size);
105 
106  return wkt;
107 }
108 
109 static rt_raster fillRasterToPolygonize(int hasnodata, double nodataval) {
110  rt_band band = NULL;
111  rt_pixtype pixtype = PT_32BF;
112 
113  /* Create raster */
114  uint16_t width = 9;
115  uint16_t height = 9;
116 
117  rt_raster raster = rt_raster_new(width, height);
118  rt_raster_set_scale(raster, 1, 1);
119 
120  band = cu_add_band(raster, pixtype, hasnodata, nodataval);
121  CU_ASSERT(band != NULL);
122 
123  {
124  int x, y;
125  for (x = 0; x < rt_band_get_width(band); ++x)
126  for (y = 0; y < rt_band_get_height(band); ++y)
127  rt_band_set_pixel(band, x, y, 0.0, NULL);
128  }
129 
130  rt_band_set_pixel(band, 3, 1, 1.8, NULL);
131  rt_band_set_pixel(band, 4, 1, 1.8, NULL);
132  rt_band_set_pixel(band, 5, 1, 2.8, NULL);
133  rt_band_set_pixel(band, 2, 2, 1.8, NULL);
134  rt_band_set_pixel(band, 3, 2, 1.8, NULL);
135  rt_band_set_pixel(band, 4, 2, 1.8, NULL);
136  rt_band_set_pixel(band, 5, 2, 2.8, NULL);
137  rt_band_set_pixel(band, 6, 2, 2.8, NULL);
138  rt_band_set_pixel(band, 1, 3, 1.8, NULL);
139  rt_band_set_pixel(band, 2, 3, 1.8, NULL);
140  rt_band_set_pixel(band, 6, 3, 2.8, NULL);
141  rt_band_set_pixel(band, 7, 3, 2.8, NULL);
142  rt_band_set_pixel(band, 1, 4, 1.8, NULL);
143  rt_band_set_pixel(band, 2, 4, 1.8, NULL);
144  rt_band_set_pixel(band, 6, 4, 2.8, NULL);
145  rt_band_set_pixel(band, 7, 4, 2.8, NULL);
146  rt_band_set_pixel(band, 1, 5, 1.8, NULL);
147  rt_band_set_pixel(band, 2, 5, 1.8, NULL);
148  rt_band_set_pixel(band, 6, 5, 2.8, NULL);
149  rt_band_set_pixel(band, 7, 5, 2.8, NULL);
150  rt_band_set_pixel(band, 2, 6, 1.8, NULL);
151  rt_band_set_pixel(band, 3, 6, 1.8, NULL);
152  rt_band_set_pixel(band, 4, 6, 1.8, NULL);
153  rt_band_set_pixel(band, 5, 6, 2.8, NULL);
154  rt_band_set_pixel(band, 6, 6, 2.8, NULL);
155  rt_band_set_pixel(band, 3, 7, 1.8, NULL);
156  rt_band_set_pixel(band, 4, 7, 1.8, NULL);
157  rt_band_set_pixel(band, 5, 7, 2.8, NULL);
158 
159  return raster;
160 }
161 
162 static void test_gdal_polygonize() {
163  int i;
164  rt_raster rt;
165  int nPols = 0;
166  rt_geomval gv = NULL;
167  char *wkt = NULL;
168 
169  rt = fillRasterToPolygonize(1, -1.0);
170  CU_ASSERT(rt_raster_has_band(rt, 0));
171 
172  nPols = 0;
173  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
174 
175 #ifdef GDALFPOLYGONIZE
176  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
177 #else
178  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 2.0, 1.);
179 #endif
180 
181  wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom);
182  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))");
183  rtdealloc(wkt);
184 
185  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
186  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
187  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))");
188  rtdealloc(wkt);
189 
190 #ifdef GDALFPOLYGONIZE
191  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 2.8, FLT_EPSILON);
192 #else
193  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 3.0, 1.);
194 #endif
195 
196  wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom);
197  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))");
198  rtdealloc(wkt);
199 
200  CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
201  wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom);
202  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))");
203  rtdealloc(wkt);
204 
205  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
206  rtdealloc(gv);
207  cu_free_raster(rt);
208 
209  /* Second test: NODATA value = 1.8 */
210 #ifdef GDALFPOLYGONIZE
211  rt = fillRasterToPolygonize(1, 1.8);
212 #else
213  rt = fillRasterToPolygonize(1, 2.0);
214 #endif
215 
216  /* We can check rt_raster_has_band here too */
217  CU_ASSERT(rt_raster_has_band(rt, 0));
218 
219  nPols = 0;
220  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
221 
222  /*
223  for (i = 0; i < nPols; i++) {
224  wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
225  printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
226  rtdealloc(wkt);
227  }
228  */
229 
230 #ifdef GDALFPOLYGONIZE
231  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
232  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
233  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))");
234  rtdealloc(wkt);
235 
236  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 2.8, FLT_EPSILON);
237  wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom);
238  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))");
239  rtdealloc(wkt);
240 
241  CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
242  wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom);
243  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))");
244  rtdealloc(wkt);
245 #else
246  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 0.0, 1.);
247  wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom);
248  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))");
249  rtdealloc(wkt);
250 
251  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 3.0, 1.);
252  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
253  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))");
254  rtdealloc(wkt);
255 
256  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 0.0, 1.);
257  wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom);
258  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))");
259  rtdealloc(wkt);
260 #endif
261 
262  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
263  rtdealloc(gv);
264  cu_free_raster(rt);
265 
266  /* Third test: NODATA value = 2.8 */
267 #ifdef GDALFPOLYGONIZE
268  rt = fillRasterToPolygonize(1, 2.8);
269 #else
270  rt = fillRasterToPolygonize(1, 3.0);
271 #endif
272 
273  /* We can check rt_raster_has_band here too */
274  CU_ASSERT(rt_raster_has_band(rt, 0));
275 
276  nPols = 0;
277  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
278 
279  /*
280  for (i = 0; i < nPols; i++) {
281  wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
282  printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
283  rtdealloc(wkt);
284  }
285  */
286 
287 #ifdef GDALFPOLYGONIZE
288  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
289 
290  CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
291  wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom);
292  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))");
293  rtdealloc(wkt);
294 #else
295  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 2.0, 1.);
296 
297  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 0.0, 1.);
298  wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom);
299  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))");
300  rtdealloc(wkt);
301 #endif
302 
303  wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom);
304  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))");
305  rtdealloc(wkt);
306 
307  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
308  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
309  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))");
310  rtdealloc(wkt);
311 
312  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
313  rtdealloc(gv);
314  cu_free_raster(rt);
315 
316  /* Fourth test: NODATA value = 0 */
317  rt = fillRasterToPolygonize(1, 0.0);
318  /* We can check rt_raster_has_band here too */
319  CU_ASSERT(rt_raster_has_band(rt, 0));
320 
321  nPols = 0;
322  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
323 
324  /*
325  for (i = 0; i < nPols; i++) {
326  wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
327  printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
328  rtdealloc(wkt);
329  }
330  */
331 
332 #ifdef GDALFPOLYGONIZE
333  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
334 #else
335  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 2.0, 1.);
336 #endif
337 
338  wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom);
339  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))");
340  rtdealloc(wkt);
341 
342 #ifdef GDALFPOLYGONIZE
343  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 2.8, FLT_EPSILON);
344 #else
345  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 3.0, 1.);
346 #endif
347 
348  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
349  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))");
350  rtdealloc(wkt);
351 
352  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
353  rtdealloc(gv);
354  cu_free_raster(rt);
355 
356  /* Last test: There is no NODATA value (all values are valid) */
357  rt = fillRasterToPolygonize(0, 0.0);
358  /* We can check rt_raster_has_band here too */
359  CU_ASSERT(rt_raster_has_band(rt, 0));
360 
361  nPols = 0;
362  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
363 
364  /*
365  for (i = 0; i < nPols; i++) {
366  wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
367  printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
368  rtdealloc(wkt);
369  }
370  */
371 
372 #ifdef GDALFPOLYGONIZE
373  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
374 #else
375  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 2.0, 1.);
376 #endif
377 
378  wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom);
379  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5 1,3 1))");
380  rtdealloc(wkt);
381 
382  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
383  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
384  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))");
385  rtdealloc(wkt);
386 
387 #ifdef GDALFPOLYGONIZE
388  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 2.8, FLT_EPSILON);
389 #else
390  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 3.0, 1.);
391 #endif
392 
393  wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom);
394  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5 1))");
395  rtdealloc(wkt);
396 
397  CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
398  wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom);
399  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3 2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))");
400  rtdealloc(wkt);
401 
402  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
403  rtdealloc(gv);
404  cu_free_raster(rt);
405 }
406 
407 static void test_raster_to_gdal() {
408  rt_pixtype pixtype = PT_64BF;
409  rt_raster raster = NULL;
410  rt_band band = NULL;
411  uint32_t x;
412  uint32_t width = 100;
413  uint32_t y;
414  uint32_t height = 100;
415  char srs[] = "PROJCS[\"unnamed\",GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6370997,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center\",45],PARAMETER[\"longitude_of_center\",-100],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1],AUTHORITY[\"EPSG\",\"2163\"]]";
416 
417  uint64_t gdalSize;
418  uint8_t *gdal = NULL;
419 
420  raster = rt_raster_new(width, height);
421  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
422 
423  band = cu_add_band(raster, pixtype, 1, 0);
424  CU_ASSERT(band != NULL);
425 
426  rt_raster_set_offsets(raster, -500000, 600000);
427  rt_raster_set_scale(raster, 1000, -1000);
428 
429  for (x = 0; x < width; x++) {
430  for (y = 0; y < height; y++) {
431  rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
432  }
433  }
434 
435  gdal = rt_raster_to_gdal(raster, srs, "GTiff", NULL, &gdalSize);
436  /*printf("gdalSize: %d\n", (int) gdalSize);*/
437  CU_ASSERT(gdalSize);
438 
439  /*
440  FILE *fh = NULL;
441  fh = fopen("/tmp/out.tif", "w");
442  fwrite(gdal, sizeof(uint8_t), gdalSize, fh);
443  fclose(fh);
444  */
445 
446  if (gdal) CPLFree(gdal);
447 
448  cu_free_raster(raster);
449 
450  raster = rt_raster_new(width, height);
451  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
452 
453  band = cu_add_band(raster, pixtype, 1, 0);
454  CU_ASSERT(band != NULL);
455 
456  rt_raster_set_offsets(raster, -500000, 600000);
457  rt_raster_set_scale(raster, 1000, -1000);
458 
459  for (x = 0; x < width; x++) {
460  for (y = 0; y < height; y++) {
461  rt_band_set_pixel(band, x, y, x, NULL);
462  }
463  }
464 
465  /* add check that band isn't NODATA */
466  CU_ASSERT_EQUAL(rt_band_check_is_nodata(band), FALSE);
467 
468  gdal = rt_raster_to_gdal(raster, srs, "PNG", NULL, &gdalSize);
469  /*printf("gdalSize: %d\n", (int) gdalSize);*/
470  CU_ASSERT(gdalSize);
471 
472  if (gdal) CPLFree(gdal);
473 
474  cu_free_raster(raster);
475 }
476 
477 static void test_gdal_to_raster() {
478  rt_pixtype pixtype = PT_64BF;
479  rt_band band = NULL;
480 
482  rt_raster rast;
483  const uint32_t width = 100;
484  const uint32_t height = 100;
485  uint32_t x;
486  uint32_t y;
487  int v;
488  double values[width][height];
489  int rtn = 0;
490  double value;
491 
492  GDALDriverH gddrv = NULL;
493  int destroy = 0;
494  GDALDatasetH gdds = NULL;
495 
496  raster = rt_raster_new(width, height);
497  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
498 
499  band = cu_add_band(raster, pixtype, 1, 0);
500  CU_ASSERT(band != NULL);
501 
502  for (x = 0; x < width; x++) {
503  for (y = 0; y < height; y++) {
504  values[x][y] = (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1);
505  rt_band_set_pixel(band, x, y, values[x][y], NULL);
506  }
507  }
508 
509  gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
510  CU_ASSERT(gddrv != NULL);
511  CU_ASSERT_EQUAL(destroy, 0);
512  CU_ASSERT(gdds != NULL);
513  CU_ASSERT_EQUAL(GDALGetRasterXSize(gdds), width);
514  CU_ASSERT_EQUAL(GDALGetRasterYSize(gdds), height);
515 
516  rast = rt_raster_from_gdal_dataset(gdds);
517  CU_ASSERT(rast != NULL);
518  CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
519 
520  band = rt_raster_get_band(rast, 0);
521  CU_ASSERT(band != NULL);
522 
523  for (x = 0; x < width; x++) {
524  for (y = 0; y < height; y++) {
525  rtn = rt_band_get_pixel(band, x, y, &value, NULL);
526  CU_ASSERT_EQUAL(rtn, ES_NONE);
527  CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], DBL_EPSILON);
528  }
529  }
530 
531  GDALClose(gdds);
532  gdds = NULL;
533  gddrv = NULL;
534 
535  cu_free_raster(rast);
536  cu_free_raster(raster);
537 
538  raster = rt_raster_new(width, height);
539  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
540 
541  pixtype = PT_8BSI;
542  band = cu_add_band(raster, pixtype, 1, 0);
543  CU_ASSERT(band != NULL);
544 
545  v = -127;
546  for (x = 0; x < width; x++) {
547  for (y = 0; y < height; y++) {
548  values[x][y] = v++;
549  rt_band_set_pixel(band, x, y, values[x][y], NULL);
550  if (v == 128)
551  v = -127;
552  }
553  }
554 
555  gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
556  CU_ASSERT(gddrv != NULL);
557  CU_ASSERT_EQUAL(destroy, 0);
558  CU_ASSERT(gdds != NULL);
559  CU_ASSERT_EQUAL(GDALGetRasterXSize(gdds), width);
560  CU_ASSERT_EQUAL(GDALGetRasterYSize(gdds), height);
561 
562  rast = rt_raster_from_gdal_dataset(gdds);
563  CU_ASSERT(rast != NULL);
564  CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
565 
566  band = rt_raster_get_band(rast, 0);
567  CU_ASSERT(band != NULL);
568  CU_ASSERT_EQUAL(rt_band_get_pixtype(band), PT_16BSI);
569 
570  for (x = 0; x < width; x++) {
571  for (y = 0; y < height; y++) {
572  rtn = rt_band_get_pixel(band, x, y, &value, NULL);
573  CU_ASSERT_EQUAL(rtn, ES_NONE);
574  CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], 1.);
575  }
576  }
577 
578  GDALClose(gdds);
579  gdds = NULL;
580  gddrv = NULL;
581 
582  cu_free_raster(rast);
583  cu_free_raster(raster);
584 }
585 
586 static void test_gdal_warp() {
587  rt_pixtype pixtype = PT_64BF;
588  rt_band band = NULL;
589 
591  rt_raster rast;
592  uint32_t x;
593  uint32_t width = 100;
594  uint32_t y;
595  uint32_t height = 100;
596  double value = 0;
597 
598  char src_srs[] = "PROJCS[\"unnamed\",GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6370997,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center\",45],PARAMETER[\"longitude_of_center\",-100],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1],AUTHORITY[\"EPSG\",\"2163\"]]";
599 
600  char dst_srs[] = "PROJCS[\"NAD83 / California Albers\",GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4269\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Albers_Conic_Equal_Area\"],PARAMETER[\"standard_parallel_1\",34],PARAMETER[\"standard_parallel_2\",40.5],PARAMETER[\"latitude_of_center\",0],PARAMETER[\"longitude_of_center\",-120],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",-4000000],AUTHORITY[\"EPSG\",\"3310\"],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH]]";
601 
602  raster = rt_raster_new(width, height);
603  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
604 
605  band = cu_add_band(raster, pixtype, 1, 0);
606  CU_ASSERT(band != NULL);
607 
608  rt_raster_set_offsets(raster, -500000, 600000);
609  rt_raster_set_scale(raster, 1000, -1000);
610 
611  for (x = 0; x < width; x++) {
612  for (y = 0; y < height; y++) {
613  rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
614  }
615  }
616 
617  rast = rt_raster_gdal_warp(
618  raster,
619  src_srs, dst_srs,
620  NULL, NULL,
621  NULL, NULL,
622  NULL, NULL,
623  NULL, NULL,
624  NULL, NULL,
625  GRA_NearestNeighbour, -1
626  );
627  CU_ASSERT(rast != NULL);
628  CU_ASSERT_EQUAL(rt_raster_get_width(rast), 122);
629  CU_ASSERT_EQUAL(rt_raster_get_height(rast), 116);
630  CU_ASSERT_NOT_EQUAL(rt_raster_get_num_bands(rast), 0);
631 
632  band = rt_raster_get_band(rast, 0);
633  CU_ASSERT(band != NULL);
634 
635  CU_ASSERT(rt_band_get_hasnodata_flag(band));
636  rt_band_get_nodata(band, &value);
637  CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
638 
639  CU_ASSERT_EQUAL(rt_band_get_pixel(band, 0, 0, &value, NULL), ES_NONE);
640  CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
641 
642  cu_free_raster(rast);
643  cu_free_raster(raster);
644 }
645 
646 /* register tests */
647 void gdal_suite_setup(void);
649 {
650  CU_pSuite suite = CU_add_suite("gdal", NULL, NULL);
657  PG_ADD_TEST(suite, test_gdal_warp);
658 }
659 
static void test_gdal_to_raster()
Definition: cu_gdal.c:477
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
static void test_gdal_drivers()
Definition: cu_gdal.c:31
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:669
static void test_gdal_polygonize()
Definition: cu_gdal.c:162
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
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1063
struct rt_gdaldriver_t * rt_gdaldriver
Definition: librtcore.h:166
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
tuple rt
Definition: rtrowdump.py:98
tuple band
Definition: ovdump.py:57
rt_raster rt_raster_gdal_rasterize(const unsigned char *wkb, uint32_t wkb_len, const char *srs, uint32_t num_bands, rt_pixtype *pixtype, double *init, double *value, double *nodata, uint8_t *hasnodata, int *width, int *height, double *scale_x, double *scale_y, double *ul_xw, double *ul_yw, double *grid_xw, double *grid_yw, double *skew_x, double *skew_y, char **options)
Return a raster of the provided geometry.
Definition: rt_raster.c:2492
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
tuple rast
Definition: rtpixdump.py:61
rt_raster rt_raster_gdal_warp(rt_raster raster, const char *src_srs, const char *dst_srs, double *scale_x, double *scale_y, int *width, int *height, double *ul_xw, double *ul_yw, double *grid_xw, double *grid_yw, double *skew_x, double *skew_y, GDALResampleAlg resample_alg, double max_err)
Return a warped raster using GDAL Warp API.
Definition: rt_warp.c:178
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
static void test_raster_to_gdal()
Definition: cu_gdal.c:407
rt_pixtype
Definition: librtcore.h:197
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
void cu_free_raster(rt_raster raster)
rt_band cu_add_band(rt_raster raster, rt_pixtype pixtype, int hasnodata, double nodataval)
static char * lwgeom_to_text(const LWGEOM *lwgeom)
Definition: cu_gdal.c:100
#define WKT_ISO
Definition: liblwgeom.h:2055
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
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:137
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:199
#define PG_ADD_TEST(suite, testfunc)
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
int rt_band_check_is_nodata(rt_band band)
Returns TRUE if the band is only nodata values.
Definition: rt_band.c:1619
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_band.c:507
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
Definition: rt_raster.c:1351
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_band.c:516
static rt_raster fillRasterToPolygonize(int hasnodata, double nodataval)
Definition: cu_gdal.c:109
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
static void test_gdal_warp()
Definition: cu_gdal.c:586
tuple x
Definition: pixval.py:53
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
void rtdealloc(void *mem)
Definition: rt_context.c:186
void gdal_suite_setup(void)
Definition: cu_gdal.c:648
uint8_t * rt_raster_to_gdal(rt_raster raster, const char *srs, char *format, char **options, uint64_t *gdalsize)
Return formatted GDAL raster from raster.
Definition: rt_raster.c:1602
#define FALSE
Definition: dbfopen.c:168
rt_gdaldriver rt_raster_gdal_drivers(uint32_t *drv_count, uint8_t cancc)
Returns a set of available GDAL drivers.
Definition: rt_raster.c:1705
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:498
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
Definition: rt_band.c:841
int rt_util_gdal_configured(void)
Definition: rt_util.c:313
rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds)
Return a raster from a GDAL dataset.
Definition: rt_raster.c:2165
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
static void test_gdal_rasterize()
Definition: cu_gdal.c:49
#define TRUE
Definition: dbfopen.c:169
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
static void test_gdal_configured()
Definition: cu_gdal.c:27