PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches
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 * Copyright (C) 2025 Darafei Praliaskouski <me@komzpa.net>
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License
11 * as published by the Free Software Foundation; either version 2
12 * of the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software Foundation,
21 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 *
23 */
24
25#include "CUnit/Basic.h"
26#include "cu_tester.h"
27
28static void test_gdal_configured() {
29 CU_ASSERT(rt_util_gdal_configured());
30}
31
32static void test_gdal_drivers() {
33 uint32_t i;
34 uint32_t size;
35 rt_gdaldriver drv = NULL;
36
37 drv = (rt_gdaldriver) rt_raster_gdal_drivers(&size, 1);
38 CU_ASSERT(drv != NULL);
39
40 for (i = 0; i < size; i++) {
41 CU_ASSERT(strlen(drv[i].short_name));
42 rtdealloc(drv[i].short_name);
43 rtdealloc(drv[i].long_name);
44 rtdealloc(drv[i].create_options);
45 }
46
47 rtdealloc(drv);
48}
49
50static void test_gdal_rasterize() {
51 rt_raster raster;
52 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\"]]";
53 const char wkb_hex[] = "010300000001000000050000000000000080841ec100000000600122410000000080841ec100000000804f22410000000040e81dc100000000804f22410000000040e81dc100000000600122410000000080841ec10000000060012241";
54 const char *pos = wkb_hex;
55 unsigned char *wkb = NULL;
56 int wkb_len = 0;
57 int i;
58 double scale_x = 100;
59 double scale_y = -100;
60
61 rt_pixtype pixtype[] = {PT_8BUI};
62 double init[] = {0};
63 double value[] = {1};
64 double nodata[] = {0};
65 uint8_t nodata_mask[] = {1};
66
67 /* hex to byte */
68 wkb_len = (int) ceil(((double) strlen(wkb_hex)) / 2);
69 wkb = (unsigned char *) rtalloc(sizeof(unsigned char) * wkb_len);
70 for (i = 0; i < wkb_len; i++) {
71 int b = 0;
72 sscanf(pos, "%2x", &b);
73 wkb[i] = (unsigned char)b;
74 pos += 2;
75 }
76
77 raster = rt_raster_gdal_rasterize(
78 wkb,
79 wkb_len, srs,
80 1, pixtype,
81 init, value,
82 nodata, nodata_mask,
83 NULL, NULL,
84 &scale_x, &scale_y,
85 NULL, NULL,
86 NULL, NULL,
87 NULL, NULL,
88 NULL
89 );
90
91 CU_ASSERT(raster != NULL);
92 CU_ASSERT_EQUAL(rt_raster_get_width(raster), 100);
93 CU_ASSERT_EQUAL(rt_raster_get_height(raster), 100);
94 CU_ASSERT_NOT_EQUAL(rt_raster_get_num_bands(raster), 0);
95 CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_x_offset(raster), -500000, DBL_EPSILON);
96 CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_y_offset(raster), 600000, DBL_EPSILON);
97
98 rtdealloc(wkb);
99 cu_free_raster(raster);
100}
101
102static rt_raster fillRasterToPolygonize(int hasnodata, double nodataval) {
103 rt_band band = NULL;
104 rt_pixtype pixtype = PT_32BF;
105
106 /* Create raster */
107 uint16_t width = 9;
108 uint16_t height = 9;
109
110 rt_raster raster = rt_raster_new(width, height);
111 rt_raster_set_scale(raster, 1, 1);
112
113 band = cu_add_band(raster, pixtype, hasnodata, nodataval);
114 CU_ASSERT(band != NULL);
115
116 {
117 int x, y;
118 for (x = 0; x < rt_band_get_width(band); ++x)
119 for (y = 0; y < rt_band_get_height(band); ++y)
120 rt_band_set_pixel(band, x, y, 0.0, NULL);
121 }
122
123 rt_band_set_pixel(band, 3, 1, 1.8, NULL);
124 rt_band_set_pixel(band, 4, 1, 1.8, NULL);
125 rt_band_set_pixel(band, 5, 1, 2.8, NULL);
126 rt_band_set_pixel(band, 2, 2, 1.8, NULL);
127 rt_band_set_pixel(band, 3, 2, 1.8, NULL);
128 rt_band_set_pixel(band, 4, 2, 1.8, NULL);
129 rt_band_set_pixel(band, 5, 2, 2.8, NULL);
130 rt_band_set_pixel(band, 6, 2, 2.8, NULL);
131 rt_band_set_pixel(band, 1, 3, 1.8, NULL);
132 rt_band_set_pixel(band, 2, 3, 1.8, NULL);
133 rt_band_set_pixel(band, 6, 3, 2.8, NULL);
134 rt_band_set_pixel(band, 7, 3, 2.8, NULL);
135 rt_band_set_pixel(band, 1, 4, 1.8, NULL);
136 rt_band_set_pixel(band, 2, 4, 1.8, NULL);
137 rt_band_set_pixel(band, 6, 4, 2.8, NULL);
138 rt_band_set_pixel(band, 7, 4, 2.8, NULL);
139 rt_band_set_pixel(band, 1, 5, 1.8, NULL);
140 rt_band_set_pixel(band, 2, 5, 1.8, NULL);
141 rt_band_set_pixel(band, 6, 5, 2.8, NULL);
142 rt_band_set_pixel(band, 7, 5, 2.8, NULL);
143 rt_band_set_pixel(band, 2, 6, 1.8, NULL);
144 rt_band_set_pixel(band, 3, 6, 1.8, NULL);
145 rt_band_set_pixel(band, 4, 6, 1.8, NULL);
146 rt_band_set_pixel(band, 5, 6, 2.8, NULL);
147 rt_band_set_pixel(band, 6, 6, 2.8, NULL);
148 rt_band_set_pixel(band, 3, 7, 1.8, NULL);
149 rt_band_set_pixel(band, 4, 7, 1.8, NULL);
150 rt_band_set_pixel(band, 5, 7, 2.8, NULL);
151
152 return raster;
153}
154
155static void test_gdal_polygonize() {
156 int i;
157 rt_raster rt;
158 int nPols = 0;
159 double total_area = 0;
160 double total_val = 0;
161 rt_geomval gv = NULL;
162 LWGEOM *gobserved;
163 //char *wkt = NULL;
164
165 rt = fillRasterToPolygonize(1, -1.0);
166 CU_ASSERT(rt_raster_has_band(rt, 0));
167
168 nPols = 0;
169 gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
170 CU_ASSERT_DOUBLE_EQUAL(nPols, 4, FLT_EPSILON);
171 total_area = 0; total_val = 0;
172 for (i = 0; i < nPols; i++) {
173 total_val += gv[i].val;
174 gobserved = (LWGEOM *) gv[i].geom;
175 total_area += lwgeom_area(gobserved);
176 lwgeom_free((LWGEOM *) gv[i].geom);
177 }
178 printf("total area, total val, nPols = %f, %f, %i\n", total_area, total_val, nPols);
179 CU_ASSERT_DOUBLE_EQUAL(total_val, 1.8 + 0.0 + 2.8 + 0, FLT_EPSILON);
180 CU_ASSERT_DOUBLE_EQUAL(total_area, 81, FLT_EPSILON);
181
182 rtdealloc(gv);
183 cu_free_raster(rt);
184
185 /* Second test: NODATA value = 1.8 */
186 rt = fillRasterToPolygonize(1, 1.8);
187
188 /* We can check rt_raster_has_band here too */
189 CU_ASSERT(rt_raster_has_band(rt, 0));
190
191 nPols = 0;
192 gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
193 CU_ASSERT_DOUBLE_EQUAL(nPols, 4, FLT_EPSILON);
194 total_area = 0; total_val = 0;
195 for (i = 0; i < nPols; i++) {
196 total_val += gv[i].val;
197 gobserved = (LWGEOM *) gv[i].geom;
198 total_area += lwgeom_area(gobserved);
199 lwgeom_free((LWGEOM *) gv[i].geom);
200 }
201 printf("total area, total_val, polys = %f, %f, %i\n", total_area, total_val, nPols);
202 CU_ASSERT_DOUBLE_EQUAL(total_val, 4.6, FLT_EPSILON);
203 CU_ASSERT_DOUBLE_EQUAL(total_area, 81, FLT_EPSILON);
204
205
206 rtdealloc(gv);
207 cu_free_raster(rt);
208
209 /* Third test: NODATA value = 2.8 */
210 rt = fillRasterToPolygonize(1, 2.8);
211
212 /* We can check rt_raster_has_band here too */
213 CU_ASSERT(rt_raster_has_band(rt, 0));
214
215 nPols = 0;
216 gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
217 CU_ASSERT_DOUBLE_EQUAL(nPols, 4, FLT_EPSILON);
218 total_area = 0; total_val = 0;
219 for (i = 0; i < nPols; i++) {
220 total_val += gv[i].val;
221 gobserved = (LWGEOM *) gv[i].geom;
222 total_area += lwgeom_area(gobserved);
223 lwgeom_free((LWGEOM *) gv[i].geom);
224 }
225
226 printf("total area, total_val, polys = %f, %f, %i\n", total_area, total_val, nPols);
227 CU_ASSERT_DOUBLE_EQUAL(total_val, 4.6, FLT_EPSILON);
228 CU_ASSERT_DOUBLE_EQUAL(total_area, 81, FLT_EPSILON);
229
230 rtdealloc(gv);
231 cu_free_raster(rt);
232
233 /* Fourth test: NODATA value = 0 */
234 rt = fillRasterToPolygonize(1, 0.0);
235 /* We can check rt_raster_has_band here too */
236 CU_ASSERT(rt_raster_has_band(rt, 0));
237
238 nPols = 0;
239 gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
240
241 CU_ASSERT_DOUBLE_EQUAL(nPols, 2, FLT_EPSILON);
242 total_area = 0; total_val = 0;
243 for (i = 0; i < nPols; i++) {
244 total_val += gv[i].val;
245 gobserved = (LWGEOM *) gv[i].geom;
246 total_area += lwgeom_area(gobserved);
247 lwgeom_free((LWGEOM *) gv[i].geom);
248 }
249
250 printf("total area, total_val, polys = %f, %f, %i\n", total_area, total_val, nPols);
251 CU_ASSERT_DOUBLE_EQUAL(total_val, 4.6, FLT_EPSILON);
252 CU_ASSERT_DOUBLE_EQUAL(total_area, 28, FLT_EPSILON);
253
254 rtdealloc(gv);
255 cu_free_raster(rt);
256
257 /* Last test: There is no NODATA value (all values are valid) */
258 rt = fillRasterToPolygonize(0, 0.0);
259 /* We can check rt_raster_has_band here too */
260 CU_ASSERT(rt_raster_has_band(rt, 0));
261
262 nPols = 0;
263 gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
264
265 CU_ASSERT_DOUBLE_EQUAL(nPols, 4, FLT_EPSILON);
266 total_area = 0; total_val = 0;
267 for (i = 0; i < nPols; i++) {
268 total_val += gv[i].val;
269 gobserved = (LWGEOM *) gv[i].geom;
270 total_area += lwgeom_area(gobserved);
271 lwgeom_free((LWGEOM *) gv[i].geom);
272 }
273
274 printf("total area, total_val, polys = %f, %f, %i\n", total_area, total_val, nPols);
275 CU_ASSERT_DOUBLE_EQUAL(total_val, 1.8 + 0.0 + 2.8 + 0.0, FLT_EPSILON);
276 CU_ASSERT_DOUBLE_EQUAL(total_area, 81, FLT_EPSILON);
277 rtdealloc(gv);
278 cu_free_raster(rt);
279}
280
281static void
282test_gdal_polygonize_interrupt(void)
283{
284 rt_raster rt;
285 int nPols = 0;
286 rt_geomval gv = NULL;
287
288 rt = fillRasterToPolygonize(0, 0.0);
289 CU_ASSERT(rt != NULL);
290
291 /* Why: confirm that the GDAL callback honours liblwgeom interrupts (#4222). */
292 lwgeom_request_interrupt();
293 gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
294 lwgeom_cancel_interrupt();
295
296 CU_ASSERT_PTR_NULL(gv);
297 CU_ASSERT_EQUAL(nPols, 0);
298
299 if (gv)
300 rtdealloc(gv);
301 cu_free_raster(rt);
302}
303
304static void test_raster_to_gdal() {
305 rt_pixtype pixtype = PT_64BF;
306 rt_raster raster = NULL;
307 rt_band band = NULL;
308 uint32_t x;
309 uint32_t width = 100;
310 uint32_t y;
311 uint32_t height = 100;
312 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\"]]";
313
314 uint64_t gdalSize;
315 uint8_t *gdal = NULL;
316
317 raster = rt_raster_new(width, height);
318 CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
319
320 band = cu_add_band(raster, pixtype, 1, 0);
321 CU_ASSERT(band != NULL);
322
323 rt_raster_set_offsets(raster, -500000, 600000);
324 rt_raster_set_scale(raster, 1000, -1000);
325
326 for (x = 0; x < width; x++) {
327 for (y = 0; y < height; y++) {
328 rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
329 }
330 }
331
332 gdal = rt_raster_to_gdal(raster, srs, "GTiff", NULL, &gdalSize);
333 /*printf("gdalSize: %d\n", (int) gdalSize);*/
334 CU_ASSERT(gdalSize);
335
336 /*
337 FILE *fh = NULL;
338 fh = fopen("/tmp/out.tif", "w");
339 fwrite(gdal, sizeof(uint8_t), gdalSize, fh);
340 fclose(fh);
341 */
342
343 if (gdal) CPLFree(gdal);
344
345 cu_free_raster(raster);
346
347 raster = rt_raster_new(width, height);
348 CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
349
350 band = cu_add_band(raster, pixtype, 1, 0);
351 CU_ASSERT(band != NULL);
352
353 rt_raster_set_offsets(raster, -500000, 600000);
354 rt_raster_set_scale(raster, 1000, -1000);
355
356 for (x = 0; x < width; x++) {
357 for (y = 0; y < height; y++) {
358 rt_band_set_pixel(band, x, y, x, NULL);
359 }
360 }
361
362 /* add check that band isn't NODATA */
363 CU_ASSERT_EQUAL(rt_band_check_is_nodata(band), FALSE);
364
365 gdal = rt_raster_to_gdal(raster, srs, "PNG", NULL, &gdalSize);
366 /*printf("gdalSize: %d\n", (int) gdalSize);*/
367 CU_ASSERT(gdalSize);
368
369 if (gdal) CPLFree(gdal);
370
371 gdal = rt_raster_to_gdal(raster, srs, "PCIDSK", NULL, &gdalSize);
372 CU_ASSERT(gdal == NULL);
373
374 cu_free_raster(raster);
375}
376
377static void test_gdal_to_raster() {
378 rt_pixtype pixtype = PT_64BF;
379 rt_band band = NULL;
380
381 rt_raster raster;
382 rt_raster rast;
383 const uint32_t width = 100;
384 const uint32_t height = 100;
385 uint32_t x;
386 uint32_t y;
387 int v;
388 double values[width][height];
389 int rtn = 0;
390 double value;
391
392 GDALDriverH gddrv = NULL;
393 int destroy = 0;
394 GDALDatasetH gdds = NULL;
395
396 raster = rt_raster_new(width, height);
397 CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
398
399 band = cu_add_band(raster, pixtype, 1, 0);
400 CU_ASSERT(band != NULL);
401
402 for (x = 0; x < width; x++) {
403 for (y = 0; y < height; y++) {
404 values[x][y] = (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1);
405 rt_band_set_pixel(band, x, y, values[x][y], NULL);
406 }
407 }
408
409 gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
410 CU_ASSERT(gddrv != NULL);
411 CU_ASSERT(gdds != NULL);
412 CU_ASSERT_EQUAL((uint32_t)GDALGetRasterXSize(gdds), width);
413 CU_ASSERT_EQUAL((uint32_t)GDALGetRasterYSize(gdds), height);
414
415 rast = rt_raster_from_gdal_dataset(gdds);
416 CU_ASSERT(rast != NULL);
417 CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
418
419 band = rt_raster_get_band(rast, 0);
420 CU_ASSERT(band != NULL);
421
422 for (x = 0; x < width; x++) {
423 for (y = 0; y < height; y++) {
424 rtn = rt_band_get_pixel(band, x, y, &value, NULL);
425 CU_ASSERT_EQUAL(rtn, ES_NONE);
426 CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], DBL_EPSILON);
427 }
428 }
429
430 if (destroy && gddrv) {
431 GDALDeregisterDriver(gddrv);
432 GDALDestroyDriver(gddrv);
433 }
434 GDALClose(gdds);
435 gdds = NULL;
436 gddrv = NULL;
437
438 cu_free_raster(rast);
439 cu_free_raster(raster);
440
441 raster = rt_raster_new(width, height);
442 CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
443
444 pixtype = PT_8BSI;
445 band = cu_add_band(raster, pixtype, 1, 0);
446 CU_ASSERT(band != NULL);
447
448 v = -127;
449 for (x = 0; x < width; x++) {
450 for (y = 0; y < height; y++) {
451 values[x][y] = v++;
452 rt_band_set_pixel(band, x, y, values[x][y], NULL);
453 if (v == 128)
454 v = -127;
455 }
456 }
457
458 gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
459 CU_ASSERT(gddrv != NULL);
460 CU_ASSERT(gdds != NULL);
461 CU_ASSERT_EQUAL((uint32_t)GDALGetRasterXSize(gdds), width);
462 CU_ASSERT_EQUAL((uint32_t)GDALGetRasterYSize(gdds), height);
463
464 rast = rt_raster_from_gdal_dataset(gdds);
465 CU_ASSERT(rast != NULL);
466 CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
467
468 band = rt_raster_get_band(rast, 0);
469 CU_ASSERT(band != NULL);
470#if POSTGIS_GDAL_VERSION < 30700
471 CU_ASSERT_EQUAL(rt_band_get_pixtype(band), PT_16BSI);
472#else
473 CU_ASSERT_EQUAL(rt_band_get_pixtype(band), PT_8BSI);
474#endif
475 for (x = 0; x < width; x++) {
476 for (y = 0; y < height; y++) {
477 rtn = rt_band_get_pixel(band, x, y, &value, NULL);
478#if POSTGIS_GDAL_VERSION < 30700
479 CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], 1.);
480 CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], DBL_EPSILON);
481#else
482 CU_ASSERT_EQUAL(rtn, ES_NONE);
483 CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], DBL_EPSILON);
484#endif
485 }
486 }
487 if (destroy && gddrv) {
488 GDALDeregisterDriver(gddrv);
489 GDALDestroyDriver(gddrv);
490 }
491
492 GDALClose(gdds);
493 gdds = NULL;
494 gddrv = NULL;
495
496 cu_free_raster(rast);
497 cu_free_raster(raster);
498}
499
500static void test_gdal_warp() {
501 rt_pixtype pixtype = PT_64BF;
502 rt_band band = NULL;
503
504 rt_raster raster;
505 rt_raster rast;
506 uint32_t x;
507 uint32_t width = 100;
508 uint32_t y;
509 uint32_t height = 100;
510 double value = 0;
511
512 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\"]]";
513
514 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]]";
515
516 raster = rt_raster_new(width, height);
517 CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
518
519 band = cu_add_band(raster, pixtype, 1, 0);
520 CU_ASSERT(band != NULL);
521
522 rt_raster_set_offsets(raster, -500000, 600000);
523 rt_raster_set_scale(raster, 1000, -1000);
524
525 for (x = 0; x < width; x++) {
526 for (y = 0; y < height; y++) {
527 rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
528 }
529 }
530
531 rast = rt_raster_gdal_warp(
532 raster,
533 src_srs, dst_srs,
534 NULL, NULL,
535 NULL, NULL,
536 NULL, NULL,
537 NULL, NULL,
538 NULL, NULL,
539 GRA_NearestNeighbour, -1
540 );
541 CU_ASSERT(rast != NULL);
542 CU_ASSERT_EQUAL(rt_raster_get_width(rast), 122);
543 CU_ASSERT_EQUAL(rt_raster_get_height(rast), 116);
544 CU_ASSERT_NOT_EQUAL(rt_raster_get_num_bands(rast), 0);
545
546 band = rt_raster_get_band(rast, 0);
547 CU_ASSERT(band != NULL);
548
549 CU_ASSERT(rt_band_get_hasnodata_flag(band));
550 rt_band_get_nodata(band, &value);
551 CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
552
553 CU_ASSERT_EQUAL(rt_band_get_pixel(band, 0, 0, &value, NULL), ES_NONE);
554 CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
555
556 cu_free_raster(rast);
557 cu_free_raster(raster);
558}
559
560static void test_gdal_warp_preserves_data(void) {
561 const char *filename = POSTGIS_TOP_SRC_DIR "/raster/test/regress/loader/Projected.tif";
562
563 GDALDatasetH hDS_in = NULL;
564 rt_raster rast_in = NULL;
565 rt_raster rast_out = NULL;
566 int band_count_in, band_count_out, i;
567
568 // double scale_x = 0.0, scale_y = 0.0;
569 // double dim_x = 0.0, dim_y = 0.0;
570 // int width = 0, height = 0;
571 // double grid_xw = 0.0, grid_yw = 0.0;
572 // double skew_x = 0.0, skew_y = 0.0;
573
574 double max_err = 0.125;
575 GDALResampleAlg alg = GRA_NearestNeighbour;
576
577 const char *src_srs = "EPSG:4326";
578 const char *dst_srs = "EPSG:3857";
579
580 /* Handle to TIFF */
581 GDALAllRegister();
582 hDS_in = GDALOpen(filename, GA_ReadOnly);
583 CU_ASSERT(hDS_in != NULL);
584
585 /* Read TIFF into memory as rt_raster */
586 rast_in = rt_raster_from_gdal_dataset(hDS_in);
587 CU_ASSERT(rast_in != NULL);
588
589 /* Warp raster using default options */
590 rast_out = rt_raster_gdal_warp(rast_in,
591 src_srs, dst_srs,
592 NULL, NULL, // &scale_x, &scale_y,
593 NULL, NULL, // &dim_x, &dim_y,
594 NULL, NULL, // &width, &height,
595 NULL, NULL, // &grid_xw, &grid_yw,
596 NULL, NULL, // &skew_x, &skew_y,
597 alg, max_err);
598 CU_ASSERT(rast_out != NULL);
599
600 band_count_in = rt_raster_get_num_bands(rast_in);
601 band_count_out = rt_raster_get_num_bands(rast_out);
602 CU_ASSERT_EQUAL(band_count_in, band_count_out);
603
604 for (i = 0; i < band_count_in; i++) {
605 double tolerance = 0.1;
606 rt_bandstats stats_in, stats_out;
607 rt_band band_in = rt_raster_get_band(rast_in, i);
608 rt_band band_out = rt_raster_get_band(rast_out, i);
609
610 CU_ASSERT(band_in != NULL);
611 CU_ASSERT(band_out != NULL);
612
613 stats_in = rt_band_get_summary_stats(band_in, 1, 1, 0, NULL, NULL, NULL);
614 stats_out = rt_band_get_summary_stats(band_out, 1, 1, 0, NULL, NULL, NULL);
615
616 CU_ASSERT_DOUBLE_EQUAL(stats_in->min, stats_out->min, fabs(stats_in->min) * tolerance);
617 CU_ASSERT_DOUBLE_EQUAL(stats_in->max, stats_out->max, fabs(stats_in->max) * tolerance);
618 CU_ASSERT_DOUBLE_EQUAL(stats_in->mean, stats_out->mean, fabs(stats_in->mean) * tolerance);
619 }
620
621 rt_raster_destroy(rast_in);
622 rt_raster_destroy(rast_out);
623 GDALClose(hDS_in);
624}
625
626/* register tests */
627void gdal_suite_setup(void);
628void gdal_suite_setup(void)
629{
630 CU_pSuite suite = CU_add_suite("gdal", NULL, NULL);
631 PG_ADD_TEST(suite, test_gdal_configured);
632 PG_ADD_TEST(suite, test_gdal_drivers);
633 PG_ADD_TEST(suite, test_gdal_rasterize);
634 PG_ADD_TEST(suite, test_gdal_polygonize);
635 PG_ADD_TEST(suite, test_gdal_polygonize_interrupt);
636 PG_ADD_TEST(suite, test_raster_to_gdal);
637 PG_ADD_TEST(suite, test_gdal_to_raster);
638 PG_ADD_TEST(suite, test_gdal_warp);
639 PG_ADD_TEST(suite, test_gdal_warp_preserves_data);
640}
641
fillRasterToPolygonize
static rt_raster fillRasterToPolygonize(int hasnodata, double nodataval)
Definition cu_gdal.c:102
test_gdal_warp_preserves_data
static void test_gdal_warp_preserves_data(void)
Definition cu_gdal.c:560
test_gdal_warp
static void test_gdal_warp()
Definition cu_gdal.c:500
test_gdal_rasterize
static void test_gdal_rasterize()
Definition cu_gdal.c:50
test_gdal_to_raster
static void test_gdal_to_raster()
Definition cu_gdal.c:377
test_raster_to_gdal
static void test_raster_to_gdal()
Definition cu_gdal.c:304
gdal_suite_setup
void gdal_suite_setup(void)
Definition cu_gdal.c:628
test_gdal_drivers
static void test_gdal_drivers()
Definition cu_gdal.c:32
test_gdal_configured
static void test_gdal_configured()
Definition cu_gdal.c:28
test_gdal_polygonize
static void test_gdal_polygonize()
Definition cu_gdal.c:155
test_gdal_polygonize_interrupt
static void test_gdal_polygonize_interrupt(void)
Definition cu_gdal.c:282
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_request_interrupt
void lwgeom_request_interrupt(void)
Request interruption of any running code.
Definition lwgeom_api.c:658
lwgeom_cancel_interrupt
void lwgeom_cancel_interrupt(void)
Cancel any interruption request.
Definition lwgeom_api.c:662
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition lwgeom.c:1246
lwgeom_area
double lwgeom_area(const LWGEOM *geom)
Definition lwgeom.c:1999
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:799
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:425
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_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:833
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_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:1551
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:188
PT_32BF
@ PT_32BF
Definition librtcore.h:199
PT_8BSI
@ PT_8BSI
Definition librtcore.h:192
PT_16BSI
@ PT_16BSI
Definition librtcore.h:194
PT_64BF
@ PT_64BF
Definition librtcore.h:200
PT_8BUI
@ PT_8BUI
Definition librtcore.h:193
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:976
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:2089
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:1140
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:2067
rt_band_get_pixtype
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition rt_band.c:790
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:111
rt_band_get_height
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition rt_band.c:808
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
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:2515
rt_bandstats_t::mean
double mean
Definition librtcore.h:2567
rt_bandstats_t::max
double max
Definition librtcore.h:2565
rt_bandstats_t::min
double min
Definition librtcore.h:2564
rt_bandstats_t
Definition librtcore.h:2560
rt_gdaldriver_t
Definition librtcore.h:2687
rt_geomval_t::val
double val
Definition librtcore.h:2556
rt_geomval_t
Definition librtcore.h:2554
rt_raster_t
Definition librtcore.h:2486