PostGIS  3.7.0dev-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  uint32_t 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  int b = 0;
71  sscanf(pos, "%2x", &b);
72  wkb[i] = (unsigned char)b;
73  pos += 2;
74  }
75 
76  raster = rt_raster_gdal_rasterize(
77  wkb,
78  wkb_len, srs,
79  1, pixtype,
80  init, value,
81  nodata, nodata_mask,
82  NULL, NULL,
83  &scale_x, &scale_y,
84  NULL, NULL,
85  NULL, NULL,
86  NULL, NULL,
87  NULL
88  );
89 
90  CU_ASSERT(raster != NULL);
91  CU_ASSERT_EQUAL(rt_raster_get_width(raster), 100);
92  CU_ASSERT_EQUAL(rt_raster_get_height(raster), 100);
93  CU_ASSERT_NOT_EQUAL(rt_raster_get_num_bands(raster), 0);
94  CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_x_offset(raster), -500000, DBL_EPSILON);
95  CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_y_offset(raster), 600000, DBL_EPSILON);
96 
97  rtdealloc(wkb);
98  cu_free_raster(raster);
99 }
100 
101 static rt_raster fillRasterToPolygonize(int hasnodata, double nodataval) {
102  rt_band band = NULL;
103  rt_pixtype pixtype = PT_32BF;
104 
105  /* Create raster */
106  uint16_t width = 9;
107  uint16_t height = 9;
108 
109  rt_raster raster = rt_raster_new(width, height);
110  rt_raster_set_scale(raster, 1, 1);
111 
112  band = cu_add_band(raster, pixtype, hasnodata, nodataval);
113  CU_ASSERT(band != NULL);
114 
115  {
116  int x, y;
117  for (x = 0; x < rt_band_get_width(band); ++x)
118  for (y = 0; y < rt_band_get_height(band); ++y)
119  rt_band_set_pixel(band, x, y, 0.0, NULL);
120  }
121 
122  rt_band_set_pixel(band, 3, 1, 1.8, NULL);
123  rt_band_set_pixel(band, 4, 1, 1.8, NULL);
124  rt_band_set_pixel(band, 5, 1, 2.8, NULL);
125  rt_band_set_pixel(band, 2, 2, 1.8, NULL);
126  rt_band_set_pixel(band, 3, 2, 1.8, NULL);
127  rt_band_set_pixel(band, 4, 2, 1.8, NULL);
128  rt_band_set_pixel(band, 5, 2, 2.8, NULL);
129  rt_band_set_pixel(band, 6, 2, 2.8, NULL);
130  rt_band_set_pixel(band, 1, 3, 1.8, NULL);
131  rt_band_set_pixel(band, 2, 3, 1.8, NULL);
132  rt_band_set_pixel(band, 6, 3, 2.8, NULL);
133  rt_band_set_pixel(band, 7, 3, 2.8, NULL);
134  rt_band_set_pixel(band, 1, 4, 1.8, NULL);
135  rt_band_set_pixel(band, 2, 4, 1.8, NULL);
136  rt_band_set_pixel(band, 6, 4, 2.8, NULL);
137  rt_band_set_pixel(band, 7, 4, 2.8, NULL);
138  rt_band_set_pixel(band, 1, 5, 1.8, NULL);
139  rt_band_set_pixel(band, 2, 5, 1.8, NULL);
140  rt_band_set_pixel(band, 6, 5, 2.8, NULL);
141  rt_band_set_pixel(band, 7, 5, 2.8, NULL);
142  rt_band_set_pixel(band, 2, 6, 1.8, NULL);
143  rt_band_set_pixel(band, 3, 6, 1.8, NULL);
144  rt_band_set_pixel(band, 4, 6, 1.8, NULL);
145  rt_band_set_pixel(band, 5, 6, 2.8, NULL);
146  rt_band_set_pixel(band, 6, 6, 2.8, NULL);
147  rt_band_set_pixel(band, 3, 7, 1.8, NULL);
148  rt_band_set_pixel(band, 4, 7, 1.8, NULL);
149  rt_band_set_pixel(band, 5, 7, 2.8, NULL);
150 
151  return raster;
152 }
153 
154 static void test_gdal_polygonize() {
155  int i;
156  rt_raster rt;
157  int nPols = 0;
158  double total_area = 0;
159  double total_val = 0;
160  rt_geomval gv = NULL;
161  LWGEOM *gobserved;
162  //char *wkt = NULL;
163 
164  rt = fillRasterToPolygonize(1, -1.0);
165  CU_ASSERT(rt_raster_has_band(rt, 0));
166 
167  nPols = 0;
168  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
169  CU_ASSERT_DOUBLE_EQUAL(nPols, 4, FLT_EPSILON);
170  total_area = 0; total_val = 0;
171  for (i = 0; i < nPols; i++) {
172  total_val += gv[i].val;
173  gobserved = (LWGEOM *) gv[i].geom;
174  total_area += lwgeom_area(gobserved);
175  lwgeom_free((LWGEOM *) gv[i].geom);
176  }
177  printf("total area, total val, nPols = %f, %f, %i\n", total_area, total_val, nPols);
178  CU_ASSERT_DOUBLE_EQUAL(total_val, 1.8 + 0.0 + 2.8 + 0, FLT_EPSILON);
179  CU_ASSERT_DOUBLE_EQUAL(total_area, 81, FLT_EPSILON);
180 
181  rtdealloc(gv);
182  cu_free_raster(rt);
183 
184  /* Second test: NODATA value = 1.8 */
185  rt = fillRasterToPolygonize(1, 1.8);
186 
187  /* We can check rt_raster_has_band here too */
188  CU_ASSERT(rt_raster_has_band(rt, 0));
189 
190  nPols = 0;
191  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
192  CU_ASSERT_DOUBLE_EQUAL(nPols, 4, FLT_EPSILON);
193  total_area = 0; total_val = 0;
194  for (i = 0; i < nPols; i++) {
195  total_val += gv[i].val;
196  gobserved = (LWGEOM *) gv[i].geom;
197  total_area += lwgeom_area(gobserved);
198  lwgeom_free((LWGEOM *) gv[i].geom);
199  }
200  printf("total area, total_val, polys = %f, %f, %i\n", total_area, total_val, nPols);
201  CU_ASSERT_DOUBLE_EQUAL(total_val, 4.6, FLT_EPSILON);
202  CU_ASSERT_DOUBLE_EQUAL(total_area, 81, FLT_EPSILON);
203 
204 
205  rtdealloc(gv);
206  cu_free_raster(rt);
207 
208  /* Third test: NODATA value = 2.8 */
209  rt = fillRasterToPolygonize(1, 2.8);
210 
211  /* We can check rt_raster_has_band here too */
212  CU_ASSERT(rt_raster_has_band(rt, 0));
213 
214  nPols = 0;
215  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
216  CU_ASSERT_DOUBLE_EQUAL(nPols, 4, FLT_EPSILON);
217  total_area = 0; total_val = 0;
218  for (i = 0; i < nPols; i++) {
219  total_val += gv[i].val;
220  gobserved = (LWGEOM *) gv[i].geom;
221  total_area += lwgeom_area(gobserved);
222  lwgeom_free((LWGEOM *) gv[i].geom);
223  }
224 
225  printf("total area, total_val, polys = %f, %f, %i\n", total_area, total_val, nPols);
226  CU_ASSERT_DOUBLE_EQUAL(total_val, 4.6, FLT_EPSILON);
227  CU_ASSERT_DOUBLE_EQUAL(total_area, 81, FLT_EPSILON);
228 
229  rtdealloc(gv);
230  cu_free_raster(rt);
231 
232  /* Fourth test: NODATA value = 0 */
233  rt = fillRasterToPolygonize(1, 0.0);
234  /* We can check rt_raster_has_band here too */
235  CU_ASSERT(rt_raster_has_band(rt, 0));
236 
237  nPols = 0;
238  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
239 
240  CU_ASSERT_DOUBLE_EQUAL(nPols, 2, FLT_EPSILON);
241  total_area = 0; total_val = 0;
242  for (i = 0; i < nPols; i++) {
243  total_val += gv[i].val;
244  gobserved = (LWGEOM *) gv[i].geom;
245  total_area += lwgeom_area(gobserved);
246  lwgeom_free((LWGEOM *) gv[i].geom);
247  }
248 
249  printf("total area, total_val, polys = %f, %f, %i\n", total_area, total_val, nPols);
250  CU_ASSERT_DOUBLE_EQUAL(total_val, 4.6, FLT_EPSILON);
251  CU_ASSERT_DOUBLE_EQUAL(total_area, 28, FLT_EPSILON);
252 
253  rtdealloc(gv);
254  cu_free_raster(rt);
255 
256  /* Last test: There is no NODATA value (all values are valid) */
257  rt = fillRasterToPolygonize(0, 0.0);
258  /* We can check rt_raster_has_band here too */
259  CU_ASSERT(rt_raster_has_band(rt, 0));
260 
261  nPols = 0;
262  gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
263 
264  CU_ASSERT_DOUBLE_EQUAL(nPols, 4, FLT_EPSILON);
265  total_area = 0; total_val = 0;
266  for (i = 0; i < nPols; i++) {
267  total_val += gv[i].val;
268  gobserved = (LWGEOM *) gv[i].geom;
269  total_area += lwgeom_area(gobserved);
270  lwgeom_free((LWGEOM *) gv[i].geom);
271  }
272 
273  printf("total area, total_val, polys = %f, %f, %i\n", total_area, total_val, nPols);
274  CU_ASSERT_DOUBLE_EQUAL(total_val, 1.8 + 0.0 + 2.8 + 0.0, FLT_EPSILON);
275  CU_ASSERT_DOUBLE_EQUAL(total_area, 81, FLT_EPSILON);
276  rtdealloc(gv);
277  cu_free_raster(rt);
278 }
279 
280 static void test_raster_to_gdal() {
281  rt_pixtype pixtype = PT_64BF;
282  rt_raster raster = NULL;
283  rt_band band = NULL;
284  uint32_t x;
285  uint32_t width = 100;
286  uint32_t y;
287  uint32_t height = 100;
288  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\"]]";
289 
290  uint64_t gdalSize;
291  uint8_t *gdal = NULL;
292 
293  raster = rt_raster_new(width, height);
294  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
295 
296  band = cu_add_band(raster, pixtype, 1, 0);
297  CU_ASSERT(band != NULL);
298 
299  rt_raster_set_offsets(raster, -500000, 600000);
300  rt_raster_set_scale(raster, 1000, -1000);
301 
302  for (x = 0; x < width; x++) {
303  for (y = 0; y < height; y++) {
304  rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
305  }
306  }
307 
308  gdal = rt_raster_to_gdal(raster, srs, "GTiff", NULL, &gdalSize);
309  /*printf("gdalSize: %d\n", (int) gdalSize);*/
310  CU_ASSERT(gdalSize);
311 
312  /*
313  FILE *fh = NULL;
314  fh = fopen("/tmp/out.tif", "w");
315  fwrite(gdal, sizeof(uint8_t), gdalSize, fh);
316  fclose(fh);
317  */
318 
319  if (gdal) CPLFree(gdal);
320 
321  cu_free_raster(raster);
322 
323  raster = rt_raster_new(width, height);
324  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
325 
326  band = cu_add_band(raster, pixtype, 1, 0);
327  CU_ASSERT(band != NULL);
328 
329  rt_raster_set_offsets(raster, -500000, 600000);
330  rt_raster_set_scale(raster, 1000, -1000);
331 
332  for (x = 0; x < width; x++) {
333  for (y = 0; y < height; y++) {
334  rt_band_set_pixel(band, x, y, x, NULL);
335  }
336  }
337 
338  /* add check that band isn't NODATA */
339  CU_ASSERT_EQUAL(rt_band_check_is_nodata(band), FALSE);
340 
341  gdal = rt_raster_to_gdal(raster, srs, "PNG", NULL, &gdalSize);
342  /*printf("gdalSize: %d\n", (int) gdalSize);*/
343  CU_ASSERT(gdalSize);
344 
345  if (gdal) CPLFree(gdal);
346 
347  gdal = rt_raster_to_gdal(raster, srs, "PCIDSK", NULL, &gdalSize);
348  CU_ASSERT(gdal == NULL);
349 
350  cu_free_raster(raster);
351 }
352 
353 static void test_gdal_to_raster() {
354  rt_pixtype pixtype = PT_64BF;
355  rt_band band = NULL;
356 
357  rt_raster raster;
358  rt_raster rast;
359  const uint32_t width = 100;
360  const uint32_t height = 100;
361  uint32_t x;
362  uint32_t y;
363  int v;
364  double values[width][height];
365  int rtn = 0;
366  double value;
367 
368  GDALDriverH gddrv = NULL;
369  int destroy = 0;
370  GDALDatasetH gdds = NULL;
371 
372  raster = rt_raster_new(width, height);
373  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
374 
375  band = cu_add_band(raster, pixtype, 1, 0);
376  CU_ASSERT(band != NULL);
377 
378  for (x = 0; x < width; x++) {
379  for (y = 0; y < height; y++) {
380  values[x][y] = (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1);
381  rt_band_set_pixel(band, x, y, values[x][y], NULL);
382  }
383  }
384 
385  gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
386  CU_ASSERT(gddrv != NULL);
387  CU_ASSERT(gdds != NULL);
388  CU_ASSERT_EQUAL((uint32_t)GDALGetRasterXSize(gdds), width);
389  CU_ASSERT_EQUAL((uint32_t)GDALGetRasterYSize(gdds), height);
390 
391  rast = rt_raster_from_gdal_dataset(gdds);
392  CU_ASSERT(rast != NULL);
393  CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
394 
395  band = rt_raster_get_band(rast, 0);
396  CU_ASSERT(band != NULL);
397 
398  for (x = 0; x < width; x++) {
399  for (y = 0; y < height; y++) {
400  rtn = rt_band_get_pixel(band, x, y, &value, NULL);
401  CU_ASSERT_EQUAL(rtn, ES_NONE);
402  CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], DBL_EPSILON);
403  }
404  }
405 
406  if (destroy && gddrv) {
407  GDALDeregisterDriver(gddrv);
408  GDALDestroyDriver(gddrv);
409  }
410  GDALClose(gdds);
411  gdds = NULL;
412  gddrv = NULL;
413 
414  cu_free_raster(rast);
415  cu_free_raster(raster);
416 
417  raster = rt_raster_new(width, height);
418  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
419 
420  pixtype = PT_8BSI;
421  band = cu_add_band(raster, pixtype, 1, 0);
422  CU_ASSERT(band != NULL);
423 
424  v = -127;
425  for (x = 0; x < width; x++) {
426  for (y = 0; y < height; y++) {
427  values[x][y] = v++;
428  rt_band_set_pixel(band, x, y, values[x][y], NULL);
429  if (v == 128)
430  v = -127;
431  }
432  }
433 
434  gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
435  CU_ASSERT(gddrv != NULL);
436  CU_ASSERT(gdds != NULL);
437  CU_ASSERT_EQUAL((uint32_t)GDALGetRasterXSize(gdds), width);
438  CU_ASSERT_EQUAL((uint32_t)GDALGetRasterYSize(gdds), height);
439 
440  rast = rt_raster_from_gdal_dataset(gdds);
441  CU_ASSERT(rast != NULL);
442  CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
443 
444  band = rt_raster_get_band(rast, 0);
445  CU_ASSERT(band != NULL);
446 #if POSTGIS_GDAL_VERSION < 30700
447  CU_ASSERT_EQUAL(rt_band_get_pixtype(band), PT_16BSI);
448 #else
449  CU_ASSERT_EQUAL(rt_band_get_pixtype(band), PT_8BSI);
450 #endif
451  for (x = 0; x < width; x++) {
452  for (y = 0; y < height; y++) {
453  rtn = rt_band_get_pixel(band, x, y, &value, NULL);
454 #if POSTGIS_GDAL_VERSION < 30700
455  CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], 1.);
456  CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], DBL_EPSILON);
457 #else
458  CU_ASSERT_EQUAL(rtn, ES_NONE);
459  CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], DBL_EPSILON);
460 #endif
461  }
462  }
463  if (destroy && gddrv) {
464  GDALDeregisterDriver(gddrv);
465  GDALDestroyDriver(gddrv);
466  }
467 
468  GDALClose(gdds);
469  gdds = NULL;
470  gddrv = NULL;
471 
472  cu_free_raster(rast);
473  cu_free_raster(raster);
474 }
475 
476 static void test_gdal_warp() {
477  rt_pixtype pixtype = PT_64BF;
478  rt_band band = NULL;
479 
480  rt_raster raster;
481  rt_raster rast;
482  uint32_t x;
483  uint32_t width = 100;
484  uint32_t y;
485  uint32_t height = 100;
486  double value = 0;
487 
488  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\"]]";
489 
490  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]]";
491 
492  raster = rt_raster_new(width, height);
493  CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
494 
495  band = cu_add_band(raster, pixtype, 1, 0);
496  CU_ASSERT(band != NULL);
497 
498  rt_raster_set_offsets(raster, -500000, 600000);
499  rt_raster_set_scale(raster, 1000, -1000);
500 
501  for (x = 0; x < width; x++) {
502  for (y = 0; y < height; y++) {
503  rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
504  }
505  }
506 
507  rast = rt_raster_gdal_warp(
508  raster,
509  src_srs, dst_srs,
510  NULL, NULL,
511  NULL, NULL,
512  NULL, NULL,
513  NULL, NULL,
514  NULL, NULL,
515  GRA_NearestNeighbour, -1
516  );
517  CU_ASSERT(rast != NULL);
518  CU_ASSERT_EQUAL(rt_raster_get_width(rast), 122);
519  CU_ASSERT_EQUAL(rt_raster_get_height(rast), 116);
520  CU_ASSERT_NOT_EQUAL(rt_raster_get_num_bands(rast), 0);
521 
522  band = rt_raster_get_band(rast, 0);
523  CU_ASSERT(band != NULL);
524 
525  CU_ASSERT(rt_band_get_hasnodata_flag(band));
526  rt_band_get_nodata(band, &value);
527  CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
528 
529  CU_ASSERT_EQUAL(rt_band_get_pixel(band, 0, 0, &value, NULL), ES_NONE);
530  CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
531 
532  cu_free_raster(rast);
533  cu_free_raster(raster);
534 }
535 
536 static void test_gdal_warp_preserves_data(void) {
537  const char *filename = POSTGIS_TOP_SRC_DIR "/raster/test/regress/loader/Projected.tif";
538 
539  GDALDatasetH hDS_in = NULL;
540  rt_raster rast_in = NULL;
541  rt_raster rast_out = NULL;
542  int band_count_in, band_count_out, i;
543 
544  // double scale_x = 0.0, scale_y = 0.0;
545  // double dim_x = 0.0, dim_y = 0.0;
546  // int width = 0, height = 0;
547  // double grid_xw = 0.0, grid_yw = 0.0;
548  // double skew_x = 0.0, skew_y = 0.0;
549 
550  double max_err = 0.125;
551  GDALResampleAlg alg = GRA_NearestNeighbour;
552 
553  const char *src_srs = "EPSG:4326";
554  const char *dst_srs = "EPSG:3857";
555 
556  /* Handle to TIFF */
557  GDALAllRegister();
558  hDS_in = GDALOpen(filename, GA_ReadOnly);
559  CU_ASSERT(hDS_in != NULL);
560 
561  /* Read TIFF into memory as rt_raster */
562  rast_in = rt_raster_from_gdal_dataset(hDS_in);
563  CU_ASSERT(rast_in != NULL);
564 
565  /* Warp raster using default options */
566  rast_out = rt_raster_gdal_warp(rast_in,
567  src_srs, dst_srs,
568  NULL, NULL, // &scale_x, &scale_y,
569  NULL, NULL, // &dim_x, &dim_y,
570  NULL, NULL, // &width, &height,
571  NULL, NULL, // &grid_xw, &grid_yw,
572  NULL, NULL, // &skew_x, &skew_y,
573  alg, max_err);
574  CU_ASSERT(rast_out != NULL);
575 
576  band_count_in = rt_raster_get_num_bands(rast_in);
577  band_count_out = rt_raster_get_num_bands(rast_out);
578  CU_ASSERT_EQUAL(band_count_in, band_count_out);
579 
580  for (i = 0; i < band_count_in; i++) {
581  double tolerance = 0.1;
582  rt_bandstats stats_in, stats_out;
583  rt_band band_in = rt_raster_get_band(rast_in, i);
584  rt_band band_out = rt_raster_get_band(rast_out, i);
585 
586  CU_ASSERT(band_in != NULL);
587  CU_ASSERT(band_out != NULL);
588 
589  stats_in = rt_band_get_summary_stats(band_in, 1, 1, 0, NULL, NULL, NULL);
590  stats_out = rt_band_get_summary_stats(band_out, 1, 1, 0, NULL, NULL, NULL);
591 
592  CU_ASSERT_DOUBLE_EQUAL(stats_in->min, stats_out->min, fabs(stats_in->min) * tolerance);
593  CU_ASSERT_DOUBLE_EQUAL(stats_in->max, stats_out->max, fabs(stats_in->max) * tolerance);
594  CU_ASSERT_DOUBLE_EQUAL(stats_in->mean, stats_out->mean, fabs(stats_in->mean) * tolerance);
595  }
596 
597  rt_raster_destroy(rast_in);
598  rt_raster_destroy(rast_out);
599  GDALClose(hDS_in);
600 }
601 
602 /* register tests */
603 void gdal_suite_setup(void);
604 void gdal_suite_setup(void)
605 {
606  CU_pSuite suite = CU_add_suite("gdal", NULL, NULL);
607  PG_ADD_TEST(suite, test_gdal_configured);
608  PG_ADD_TEST(suite, test_gdal_drivers);
609  PG_ADD_TEST(suite, test_gdal_rasterize);
610  PG_ADD_TEST(suite, test_gdal_polygonize);
611  PG_ADD_TEST(suite, test_raster_to_gdal);
612  PG_ADD_TEST(suite, test_gdal_to_raster);
613  PG_ADD_TEST(suite, test_gdal_warp);
614  PG_ADD_TEST(suite, test_gdal_warp_preserves_data);
615 }
616 
fillRasterToPolygonize
static rt_raster fillRasterToPolygonize(int hasnodata, double nodataval)
Definition: cu_gdal.c:101
test_gdal_warp_preserves_data
static void test_gdal_warp_preserves_data(void)
Definition: cu_gdal.c:536
test_gdal_warp
static void test_gdal_warp()
Definition: cu_gdal.c:476
test_gdal_rasterize
static void test_gdal_rasterize()
Definition: cu_gdal.c:49
test_gdal_to_raster
static void test_gdal_to_raster()
Definition: cu_gdal.c:353
test_raster_to_gdal
static void test_raster_to_gdal()
Definition: cu_gdal.c:280
gdal_suite_setup
void gdal_suite_setup(void)
Definition: cu_gdal.c:604
test_gdal_drivers
static void test_gdal_drivers()
Definition: cu_gdal.c:31
test_gdal_configured
static void test_gdal_configured()
Definition: cu_gdal.c:27
test_gdal_polygonize
static void test_gdal_polygonize()
Definition: cu_gdal.c:154
TRUE
#define TRUE
Definition: dbfopen.c:73
FALSE
#define FALSE
Definition: dbfopen.c:72
PG_ADD_TEST
#define PG_ADD_TEST(suite, testfunc)
Definition: liblwgeom/cunit/cu_tester.h:21
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1218
lwgeom_area
double lwgeom_area(const LWGEOM *geom)
Definition: lwgeom.c:1971
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:191
rt_band_get_width
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_band.c:791
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:217
rt_util_gdal_configured
int rt_util_gdal_configured(void)
Definition: rt_util.c:328
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:1716
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:141
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:1598
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:825
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:1527
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:86
rt_pixtype
rt_pixtype
Definition: librtcore.h:187
PT_32BF
@ PT_32BF
Definition: librtcore.h:197
PT_8BSI
@ PT_8BSI
Definition: librtcore.h:191
PT_16BSI
@ PT_16BSI
Definition: librtcore.h:193
PT_64BF
@ PT_64BF
Definition: librtcore.h:198
PT_8BUI
@ PT_8BUI
Definition: librtcore.h:192
rt_raster_new
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:52
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:1253
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:2175
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:975
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:2060
rt_band_set_pixel
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:1125
ES_NONE
@ ES_NONE
Definition: librtcore.h:182
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:2502
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:376
rt_gdaldriver
struct rt_gdaldriver_t * rt_gdaldriver
Definition: librtcore.h:155
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:133
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:2038
rt_band_get_pixtype
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:782
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:125
rtdealloc
void rtdealloc(void *mem)
Definition: rt_context.c:206
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:177
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:1813
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:203
rt_band_get_summary_stats
rt_bandstats rt_band_get_summary_stats(rt_band band, int exclude_nodata_value, double sample, int inc_vals, uint64_t *cK, double *cM, double *cQ)
Compute summary statistics for a band.
Definition: rt_statistics.c:110
rt_band_get_height
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_band.c:800
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:226
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:385
genraster.value
int value
Definition: genraster.py:62
ovdump.band
band
Definition: ovdump.py:58
rtrowdump.raster
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
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:224
cu_free_raster
void cu_free_raster(rt_raster raster)
Definition: raster/test/cunit/cu_tester.c:211
cu_tester.h
rt_band_t
Definition: librtcore.h:2503
rt_bandstats_t::mean
double mean
Definition: librtcore.h:2555
rt_bandstats_t::max
double max
Definition: librtcore.h:2553
rt_bandstats_t::min
double min
Definition: librtcore.h:2552
rt_bandstats_t
Definition: librtcore.h:2548
rt_gdaldriver_t
Definition: librtcore.h:2675
rt_geomval_t::val
double val
Definition: librtcore.h:2544
rt_geomval_t
Definition: librtcore.h:2542
rt_raster_t
Definition: librtcore.h:2474