PostGIS  2.4.9dev-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() {
50  rt_raster raster;
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  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
176 
177  wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom);
178  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))");
179  rtdealloc(wkt);
180 
181  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
182  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
183  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))");
184  rtdealloc(wkt);
185 
186  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 2.8, FLT_EPSILON);
187 
188  wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom);
189  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))");
190  rtdealloc(wkt);
191 
192  CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
193  wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom);
194  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))");
195  rtdealloc(wkt);
196 
197  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
198  rtdealloc(gv);
199  cu_free_raster(rt);
200 
201  /* Second test: NODATA value = 1.8 */
202  rt = fillRasterToPolygonize(1, 1.8);
203 
204  /* We can check rt_raster_has_band here too */
205  CU_ASSERT(rt_raster_has_band(rt, 0));
206 
207  nPols = 0;
208  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
209 
210  /*
211  for (i = 0; i < nPols; i++) {
212  wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
213  printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
214  rtdealloc(wkt);
215  }
216  */
217 
218  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
219  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
220  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))");
221  rtdealloc(wkt);
222 
223  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 2.8, FLT_EPSILON);
224  wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom);
225  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))");
226  rtdealloc(wkt);
227 
228  CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
229  wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom);
230  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))");
231  rtdealloc(wkt);
232 
233  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
234  rtdealloc(gv);
235  cu_free_raster(rt);
236 
237  /* Third test: NODATA value = 2.8 */
238  rt = fillRasterToPolygonize(1, 2.8);
239 
240  /* We can check rt_raster_has_band here too */
241  CU_ASSERT(rt_raster_has_band(rt, 0));
242 
243  nPols = 0;
244  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
245 
246  /*
247  for (i = 0; i < nPols; i++) {
248  wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
249  printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
250  rtdealloc(wkt);
251  }
252  */
253 
254  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
255 
256  CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
257  wkt = lwgeom_to_text((const LWGEOM *) gv[3].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 
261  wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom);
262  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))");
263  rtdealloc(wkt);
264 
265  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
266  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
267  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))");
268  rtdealloc(wkt);
269 
270  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
271  rtdealloc(gv);
272  cu_free_raster(rt);
273 
274  /* Fourth test: NODATA value = 0 */
275  rt = fillRasterToPolygonize(1, 0.0);
276  /* We can check rt_raster_has_band here too */
277  CU_ASSERT(rt_raster_has_band(rt, 0));
278 
279  nPols = 0;
280  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
281 
282  /*
283  for (i = 0; i < nPols; i++) {
284  wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
285  printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
286  rtdealloc(wkt);
287  }
288  */
289 
290  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
291 
292  wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom);
293  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))");
294  rtdealloc(wkt);
295 
296  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 2.8, FLT_EPSILON);
297 
298  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
299  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))");
300  rtdealloc(wkt);
301 
302  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
303  rtdealloc(gv);
304  cu_free_raster(rt);
305 
306  /* Last test: There is no NODATA value (all values are valid) */
307  rt = fillRasterToPolygonize(0, 0.0);
308  /* We can check rt_raster_has_band here too */
309  CU_ASSERT(rt_raster_has_band(rt, 0));
310 
311  nPols = 0;
312  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
313 
314  /*
315  for (i = 0; i < nPols; i++) {
316  wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
317  printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
318  rtdealloc(wkt);
319  }
320  */
321 
322  CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
323 
324  wkt = lwgeom_to_text((const LWGEOM *) gv[0].geom);
325  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))");
326  rtdealloc(wkt);
327 
328  CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
329  wkt = lwgeom_to_text((const LWGEOM *) gv[1].geom);
330  CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((3 3,3 6,6 6,6 3,3 3))");
331  rtdealloc(wkt);
332 
333  CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 2.8, FLT_EPSILON);
334 
335  wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom);
336  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))");
337  rtdealloc(wkt);
338 
339  CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
340  wkt = lwgeom_to_text((const LWGEOM *) gv[3].geom);
341  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))");
342  rtdealloc(wkt);
343 
344  for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
345  rtdealloc(gv);
346  cu_free_raster(rt);
347 }
348 
349 static void test_raster_to_gdal() {
350  rt_pixtype pixtype = PT_64BF;
351  rt_raster raster = NULL;
352  rt_band band = NULL;
353  uint32_t x;
354  uint32_t width = 100;
355  uint32_t y;
356  uint32_t height = 100;
357  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\"]]";
358 
359  uint64_t gdalSize;
360  uint8_t *gdal = NULL;
361 
362  raster = rt_raster_new(width, height);
363  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
364 
365  band = cu_add_band(raster, pixtype, 1, 0);
366  CU_ASSERT(band != NULL);
367 
368  rt_raster_set_offsets(raster, -500000, 600000);
369  rt_raster_set_scale(raster, 1000, -1000);
370 
371  for (x = 0; x < width; x++) {
372  for (y = 0; y < height; y++) {
373  rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
374  }
375  }
376 
377  gdal = rt_raster_to_gdal(raster, srs, "GTiff", NULL, &gdalSize);
378  /*printf("gdalSize: %d\n", (int) gdalSize);*/
379  CU_ASSERT(gdalSize);
380 
381  /*
382  FILE *fh = NULL;
383  fh = fopen("/tmp/out.tif", "w");
384  fwrite(gdal, sizeof(uint8_t), gdalSize, fh);
385  fclose(fh);
386  */
387 
388  if (gdal) CPLFree(gdal);
389 
390  cu_free_raster(raster);
391 
392  raster = rt_raster_new(width, height);
393  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
394 
395  band = cu_add_band(raster, pixtype, 1, 0);
396  CU_ASSERT(band != NULL);
397 
398  rt_raster_set_offsets(raster, -500000, 600000);
399  rt_raster_set_scale(raster, 1000, -1000);
400 
401  for (x = 0; x < width; x++) {
402  for (y = 0; y < height; y++) {
403  rt_band_set_pixel(band, x, y, x, NULL);
404  }
405  }
406 
407  /* add check that band isn't NODATA */
408  CU_ASSERT_EQUAL(rt_band_check_is_nodata(band), FALSE);
409 
410  gdal = rt_raster_to_gdal(raster, srs, "PNG", NULL, &gdalSize);
411  /*printf("gdalSize: %d\n", (int) gdalSize);*/
412  CU_ASSERT(gdalSize);
413 
414  if (gdal) CPLFree(gdal);
415 
416  cu_free_raster(raster);
417 }
418 
419 static void test_gdal_to_raster() {
420  rt_pixtype pixtype = PT_64BF;
421  rt_band band = NULL;
422 
423  rt_raster raster;
424  rt_raster rast;
425  const uint32_t width = 100;
426  const uint32_t height = 100;
427  uint32_t x;
428  uint32_t y;
429  int v;
430  double values[width][height];
431  int rtn = 0;
432  double value;
433 
434  GDALDriverH gddrv = NULL;
435  int destroy = 0;
436  GDALDatasetH gdds = NULL;
437 
438  raster = rt_raster_new(width, height);
439  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
440 
441  band = cu_add_band(raster, pixtype, 1, 0);
442  CU_ASSERT(band != NULL);
443 
444  for (x = 0; x < width; x++) {
445  for (y = 0; y < height; y++) {
446  values[x][y] = (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1);
447  rt_band_set_pixel(band, x, y, values[x][y], NULL);
448  }
449  }
450 
451  gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
452  CU_ASSERT(gddrv != NULL);
453  CU_ASSERT_EQUAL(destroy, 0);
454  CU_ASSERT(gdds != NULL);
455  CU_ASSERT_EQUAL(GDALGetRasterXSize(gdds), width);
456  CU_ASSERT_EQUAL(GDALGetRasterYSize(gdds), height);
457 
458  rast = rt_raster_from_gdal_dataset(gdds);
459  CU_ASSERT(rast != NULL);
460  CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
461 
462  band = rt_raster_get_band(rast, 0);
463  CU_ASSERT(band != NULL);
464 
465  for (x = 0; x < width; x++) {
466  for (y = 0; y < height; y++) {
467  rtn = rt_band_get_pixel(band, x, y, &value, NULL);
468  CU_ASSERT_EQUAL(rtn, ES_NONE);
469  CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], DBL_EPSILON);
470  }
471  }
472 
473  GDALClose(gdds);
474  gdds = NULL;
475  gddrv = NULL;
476 
477  cu_free_raster(rast);
478  cu_free_raster(raster);
479 
480  raster = rt_raster_new(width, height);
481  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
482 
483  pixtype = PT_8BSI;
484  band = cu_add_band(raster, pixtype, 1, 0);
485  CU_ASSERT(band != NULL);
486 
487  v = -127;
488  for (x = 0; x < width; x++) {
489  for (y = 0; y < height; y++) {
490  values[x][y] = v++;
491  rt_band_set_pixel(band, x, y, values[x][y], NULL);
492  if (v == 128)
493  v = -127;
494  }
495  }
496 
497  gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
498  CU_ASSERT(gddrv != NULL);
499  CU_ASSERT_EQUAL(destroy, 0);
500  CU_ASSERT(gdds != NULL);
501  CU_ASSERT_EQUAL(GDALGetRasterXSize(gdds), width);
502  CU_ASSERT_EQUAL(GDALGetRasterYSize(gdds), height);
503 
504  rast = rt_raster_from_gdal_dataset(gdds);
505  CU_ASSERT(rast != NULL);
506  CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
507 
508  band = rt_raster_get_band(rast, 0);
509  CU_ASSERT(band != NULL);
510  CU_ASSERT_EQUAL(rt_band_get_pixtype(band), PT_16BSI);
511 
512  for (x = 0; x < width; x++) {
513  for (y = 0; y < height; y++) {
514  rtn = rt_band_get_pixel(band, x, y, &value, NULL);
515  CU_ASSERT_EQUAL(rtn, ES_NONE);
516  CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], 1.);
517  }
518  }
519 
520  GDALClose(gdds);
521  gdds = NULL;
522  gddrv = NULL;
523 
524  cu_free_raster(rast);
525  cu_free_raster(raster);
526 }
527 
528 static void test_gdal_warp() {
529  rt_pixtype pixtype = PT_64BF;
530  rt_band band = NULL;
531 
532  rt_raster raster;
533  rt_raster rast;
534  uint32_t x;
535  uint32_t width = 100;
536  uint32_t y;
537  uint32_t height = 100;
538  double value = 0;
539 
540  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\"]]";
541 
542  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]]";
543 
544  raster = rt_raster_new(width, height);
545  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
546 
547  band = cu_add_band(raster, pixtype, 1, 0);
548  CU_ASSERT(band != NULL);
549 
550  rt_raster_set_offsets(raster, -500000, 600000);
551  rt_raster_set_scale(raster, 1000, -1000);
552 
553  for (x = 0; x < width; x++) {
554  for (y = 0; y < height; y++) {
555  rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
556  }
557  }
558 
559  rast = rt_raster_gdal_warp(
560  raster,
561  src_srs, dst_srs,
562  NULL, NULL,
563  NULL, NULL,
564  NULL, NULL,
565  NULL, NULL,
566  NULL, NULL,
567  GRA_NearestNeighbour, -1
568  );
569  CU_ASSERT(rast != NULL);
570  CU_ASSERT_EQUAL(rt_raster_get_width(rast), 122);
571  CU_ASSERT_EQUAL(rt_raster_get_height(rast), 116);
572  CU_ASSERT_NOT_EQUAL(rt_raster_get_num_bands(rast), 0);
573 
574  band = rt_raster_get_band(rast, 0);
575  CU_ASSERT(band != NULL);
576 
577  CU_ASSERT(rt_band_get_hasnodata_flag(band));
578  rt_band_get_nodata(band, &value);
579  CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
580 
581  CU_ASSERT_EQUAL(rt_band_get_pixel(band, 0, 0, &value, NULL), ES_NONE);
582  CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
583 
584  cu_free_raster(rast);
585  cu_free_raster(raster);
586 }
587 
588 /* register tests */
589 void gdal_suite_setup(void);
590 void gdal_suite_setup(void)
591 {
592  CU_pSuite suite = CU_add_suite("gdal", NULL, NULL);
593  PG_ADD_TEST(suite, test_gdal_configured);
594  PG_ADD_TEST(suite, test_gdal_drivers);
595  PG_ADD_TEST(suite, test_gdal_rasterize);
596  PG_ADD_TEST(suite, test_gdal_polygonize);
597  PG_ADD_TEST(suite, test_raster_to_gdal);
598  PG_ADD_TEST(suite, test_gdal_to_raster);
599  PG_ADD_TEST(suite, test_gdal_warp);
600 }
601 
test_gdal_to_raster
static void test_gdal_to_raster()
Definition: cu_gdal.c:419
rt_band_t
Definition: librtcore.h:2264
rt_raster_get_x_offset
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
Definition: rt_raster.c:213
rt_raster_get_num_bands
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
test_gdal_drivers
static void test_gdal_drivers()
Definition: cu_gdal.c:31
lwgeom_to_wkt
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:669
rtrowdump.raster
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
test_gdal_polygonize
static void test_gdal_polygonize()
Definition: cu_gdal.c:162
rt_raster_gdal_polygonize
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
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1099
rt_gdaldriver
struct rt_gdaldriver_t * rt_gdaldriver
Definition: librtcore.h:154
ovdump.band
band
Definition: ovdump.py:57
rt_raster_to_gdal_mem
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
rt_raster_gdal_rasterize
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
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
rt_raster_gdal_warp
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
test_raster_to_gdal
static void test_raster_to_gdal()
Definition: cu_gdal.c:349
rt_pixtype
rt_pixtype
Definition: librtcore.h:185
rt_gdaldriver_t
Definition: librtcore.h:2422
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
uint32_t
unsigned int uint32_t
Definition: uthash.h:78
cu_free_raster
void cu_free_raster(rt_raster raster)
Definition: raster/test/cunit/cu_tester.c:209
cu_add_band
rt_band cu_add_band(rt_raster raster, rt_pixtype pixtype, int hasnodata, double nodataval)
Definition: raster/test/cunit/cu_tester.c:222
lwgeom_to_text
static char * lwgeom_to_text(const LWGEOM *lwgeom)
Definition: cu_gdal.c:100
WKT_ISO
#define WKT_ISO
Definition: liblwgeom.h:2083
rt_band_get_pixel
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
rt_raster_set_scale
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:137
rt_raster_set_offsets
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:199
PG_ADD_TEST
#define PG_ADD_TEST(suite, testfunc)
Definition: liblwgeom/cunit/cu_tester.h:17
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
rt_band_check_is_nodata
int rt_band_check_is_nodata(rt_band band)
Returns TRUE if the band is only nodata values.
Definition: rt_band.c:1619
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
rt_band_get_width
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_band.c:507
rt_raster_has_band
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
rt_band_get_height
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_band.c:516
rt_geomval_t
Definition: librtcore.h:2303
cu_tester.h
fillRasterToPolygonize
static rt_raster fillRasterToPolygonize(int hasnodata, double nodataval)
Definition: cu_gdal.c:109
rt_raster_new
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
test_gdal_warp
static void test_gdal_warp()
Definition: cu_gdal.c:528
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
rtdealloc
void rtdealloc(void *mem)
Definition: rt_context.c:186
gdal_suite_setup
void gdal_suite_setup(void)
Definition: cu_gdal.c:590
rt_raster_to_gdal
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
FALSE
#define FALSE
Definition: dbfopen.c:168
rt_raster_gdal_drivers
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_band_get_pixtype
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:498
rt_raster_t
Definition: librtcore.h:2235
rt_band_set_pixel
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel&#39;s value.
Definition: rt_band.c:841
rt_util_gdal_configured
int rt_util_gdal_configured(void)
Definition: rt_util.c:313
genraster.value
int value
Definition: genraster.py:61
rt_raster_from_gdal_dataset
rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds)
Return a raster from a GDAL dataset.
Definition: rt_raster.c:2165
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
test_gdal_rasterize
static void test_gdal_rasterize()
Definition: cu_gdal.c:49
uint8_t
unsigned char uint8_t
Definition: uthash.h:79
TRUE
#define TRUE
Definition: dbfopen.c:169
rt_raster_get_y_offset
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
Definition: rt_raster.c:222
test_gdal_configured
static void test_gdal_configured()
Definition: cu_gdal.c:27