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 
test_gdal_configured()27 static void test_gdal_configured() {
28 	CU_ASSERT(rt_util_gdal_configured());
29 }
30 
test_gdal_drivers()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 
test_gdal_rasterize()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 
fillRasterToPolygonize(int hasnodata,double nodataval)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 
test_gdal_polygonize()154 static void test_gdal_polygonize() {
155 	int i;
156 	rt_raster rt;
157 	int nPols = 0;
158 	rt_geomval gv = NULL;
159 	LWGEOM *gexpected, *gobserved;
160 	//char *wkt = NULL;
161 	gexpected = lwgeom_from_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))",
162 				   LW_PARSER_CHECK_NONE);
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 
170 	CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
171 
172 	gobserved = (LWGEOM *)gv[0].geom;
173 
174 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
175 
176 	CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
177 	gobserved = (LWGEOM *)gv[1].geom;
178 	gexpected = lwgeom_from_wkt("POLYGON((3 3,3 6,6 6,6 3,3 3))", LW_PARSER_CHECK_NONE);
179 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON );
180 
181 	CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 2.8, FLT_EPSILON);
182 	gobserved = (LWGEOM *)gv[2].geom;
183 	gexpected = lwgeom_from_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))",
184 				    LW_PARSER_CHECK_NONE);
185 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
186 
187 	gobserved = (LWGEOM *)gv[3].geom;
188 	gexpected = lwgeom_from_wkt(
189 	    "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))",
190 	    LW_PARSER_CHECK_NONE);
191 	CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
192 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
193 
194 	for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
195 	rtdealloc(gv);
196 	cu_free_raster(rt);
197 
198 	/* Second test: NODATA value = 1.8 */
199 	rt = fillRasterToPolygonize(1, 1.8);
200 
201 	/* We can check rt_raster_has_band here too */
202 	CU_ASSERT(rt_raster_has_band(rt, 0));
203 
204 	nPols = 0;
205 	gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
206 
207 	/*
208 	for (i = 0; i < nPols; i++) {
209 		wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
210 		printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
211 		rtdealloc(wkt);
212 	}
213 	*/
214 
215 	CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
216 	gobserved = (LWGEOM *)gv[1].geom;
217 	gexpected = lwgeom_from_wkt("POLYGON((3 3,3 6,6 6,6 3,3 3))", LW_PARSER_CHECK_NONE);
218 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), 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 	gobserved = (LWGEOM *)gv[2].geom;
225 	gexpected = lwgeom_from_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 				    LW_PARSER_CHECK_NONE);
227 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
228 	//wkt = lwgeom_to_text((const LWGEOM *) gv[2].geom);
229 	//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))");
230 	//rtdealloc(wkt);
231 
232 	CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
233 	gobserved = (LWGEOM *)gv[3].geom;
234 	gexpected = lwgeom_from_wkt(
235 	    "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))",
236 	    LW_PARSER_CHECK_NONE);
237 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
238 
239 	for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
240 	rtdealloc(gv);
241 	cu_free_raster(rt);
242 
243 	/* Third test: NODATA value = 2.8 */
244 	rt = fillRasterToPolygonize(1, 2.8);
245 
246 	/* We can check rt_raster_has_band here too */
247 	CU_ASSERT(rt_raster_has_band(rt, 0));
248 
249 	nPols = 0;
250 	gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
251 
252 	/*
253 	for (i = 0; i < nPols; i++) {
254 		wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
255 		printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
256 		rtdealloc(wkt);
257 	}
258 	*/
259 
260 	CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
261 	gobserved = (LWGEOM *)gv[3].geom;
262 	gexpected = lwgeom_from_wkt(
263 	    "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))",
264 	    LW_PARSER_CHECK_NONE);
265 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
266 
267 	gobserved = (LWGEOM *)gv[0].geom;
268 	gexpected = lwgeom_from_wkt(
269 	    "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))", LW_PARSER_CHECK_NONE);
270 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
271 
272 	CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
273 	gobserved = (LWGEOM *)gv[1].geom;
274 	gexpected = lwgeom_from_wkt("POLYGON((3 3,3 6,6 6,6 3,3 3))", LW_PARSER_CHECK_NONE);
275 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
276 
277 	for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
278 	rtdealloc(gv);
279 	cu_free_raster(rt);
280 
281 	/* Fourth test: NODATA value = 0 */
282 	rt = fillRasterToPolygonize(1, 0.0);
283 	/* We can check rt_raster_has_band here too */
284 	CU_ASSERT(rt_raster_has_band(rt, 0));
285 
286 	nPols = 0;
287 	gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
288 
289 	/*
290 	for (i = 0; i < nPols; i++) {
291 		wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
292 		printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
293 		rtdealloc(wkt);
294 	}
295 	*/
296 
297 	CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
298 	gobserved = (LWGEOM *)gv[0].geom;
299 	gexpected = lwgeom_from_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))", LW_PARSER_CHECK_NONE);
300 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
301 
302 	CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 2.8, FLT_EPSILON);
303 	gobserved = (LWGEOM *)gv[1].geom;
304 	gexpected = lwgeom_from_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))",
305 				    LW_PARSER_CHECK_NONE);
306 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
307 
308 	for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
309 	rtdealloc(gv);
310 	cu_free_raster(rt);
311 
312 	/* Last test: There is no NODATA value (all values are valid) */
313 	rt = fillRasterToPolygonize(0, 0.0);
314 	/* We can check rt_raster_has_band here too */
315 	CU_ASSERT(rt_raster_has_band(rt, 0));
316 
317 	nPols = 0;
318 	gv = rt_raster_gdal_polygonize(rt, 0, TRUE, &nPols);
319 
320 	/*
321 	for (i = 0; i < nPols; i++) {
322 		wkt = lwgeom_to_text((const LWGEOM *) gv[i].geom);
323 		printf("(i, val, geom) = (%d, %f, %s)\n", i, gv[i].val, wkt);
324 		rtdealloc(wkt);
325 	}
326 	*/
327 
328 	CU_ASSERT_DOUBLE_EQUAL(gv[0].val, 1.8, FLT_EPSILON);
329 	gobserved = (LWGEOM *)gv[0].geom;
330 	gexpected = lwgeom_from_wkt(
331 	    "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))", LW_PARSER_CHECK_NONE);
332 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
333 
334 	CU_ASSERT_DOUBLE_EQUAL(gv[1].val, 0.0, FLT_EPSILON);
335 	gobserved = (LWGEOM *)gv[1].geom;
336 	gexpected = lwgeom_from_wkt("POLYGON((3 3,3 6,6 6,6 3,3 3))", LW_PARSER_CHECK_NONE);
337 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
338 
339 	CU_ASSERT_DOUBLE_EQUAL(gv[2].val, 2.8, FLT_EPSILON);
340 	gobserved = (LWGEOM *)gv[2].geom;
341 	gexpected = lwgeom_from_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))",
342 				    LW_PARSER_CHECK_NONE);
343 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
344 
345 	CU_ASSERT_DOUBLE_EQUAL(gv[3].val, 0.0, FLT_EPSILON);
346 	gobserved = (LWGEOM *)gv[3].geom;
347 	gexpected = lwgeom_from_wkt(
348 	    "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))",
349 	    LW_PARSER_CHECK_NONE);
350 	CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
351 
352 	for (i = 0; i < nPols; i++) lwgeom_free((LWGEOM *) gv[i].geom);
353 	rtdealloc(gv);
354 	cu_free_raster(rt);
355 }
356 
test_raster_to_gdal()357 static void test_raster_to_gdal() {
358 	rt_pixtype pixtype = PT_64BF;
359 	rt_raster raster = NULL;
360 	rt_band band = NULL;
361 	uint32_t x;
362 	uint32_t width = 100;
363 	uint32_t y;
364 	uint32_t height = 100;
365 	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\"]]";
366 
367 	uint64_t gdalSize;
368 	uint8_t *gdal = NULL;
369 
370 	raster = rt_raster_new(width, height);
371 	CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
372 
373 	band = cu_add_band(raster, pixtype, 1, 0);
374 	CU_ASSERT(band != NULL);
375 
376 	rt_raster_set_offsets(raster, -500000, 600000);
377 	rt_raster_set_scale(raster, 1000, -1000);
378 
379 	for (x = 0; x < width; x++) {
380 		for (y = 0; y < height; y++) {
381 			rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
382 		}
383 	}
384 
385 	gdal = rt_raster_to_gdal(raster, srs, "GTiff", NULL, &gdalSize);
386 	/*printf("gdalSize: %d\n", (int) gdalSize);*/
387 	CU_ASSERT(gdalSize);
388 
389 	/*
390 	FILE *fh = NULL;
391 	fh = fopen("/tmp/out.tif", "w");
392 	fwrite(gdal, sizeof(uint8_t), gdalSize, fh);
393 	fclose(fh);
394 	*/
395 
396 	if (gdal) CPLFree(gdal);
397 
398 	cu_free_raster(raster);
399 
400 	raster = rt_raster_new(width, height);
401 	CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
402 
403 	band = cu_add_band(raster, pixtype, 1, 0);
404 	CU_ASSERT(band != NULL);
405 
406 	rt_raster_set_offsets(raster, -500000, 600000);
407 	rt_raster_set_scale(raster, 1000, -1000);
408 
409 	for (x = 0; x < width; x++) {
410 		for (y = 0; y < height; y++) {
411 			rt_band_set_pixel(band, x, y, x, NULL);
412 		}
413 	}
414 
415 	/* add check that band isn't NODATA */
416 	CU_ASSERT_EQUAL(rt_band_check_is_nodata(band), FALSE);
417 
418 	gdal = rt_raster_to_gdal(raster, srs, "PNG", NULL, &gdalSize);
419 	/*printf("gdalSize: %d\n", (int) gdalSize);*/
420 	CU_ASSERT(gdalSize);
421 
422 	if (gdal) CPLFree(gdal);
423 
424 	gdal = rt_raster_to_gdal(raster, srs, "PCIDSK", NULL, &gdalSize);
425 	CU_ASSERT(gdal == NULL);
426 
427 	cu_free_raster(raster);
428 }
429 
test_gdal_to_raster()430 static void test_gdal_to_raster() {
431 	rt_pixtype pixtype = PT_64BF;
432 	rt_band band = NULL;
433 
434 	rt_raster raster;
435 	rt_raster rast;
436 	const uint32_t width = 100;
437 	const uint32_t height = 100;
438 	uint32_t x;
439 	uint32_t y;
440 	int v;
441 	double values[width][height];
442 	int rtn = 0;
443 	double value;
444 
445 	GDALDriverH gddrv = NULL;
446 	int destroy = 0;
447 	GDALDatasetH gdds = NULL;
448 
449 	raster = rt_raster_new(width, height);
450 	CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
451 
452 	band = cu_add_band(raster, pixtype, 1, 0);
453 	CU_ASSERT(band != NULL);
454 
455 	for (x = 0; x < width; x++) {
456 		for (y = 0; y < height; y++) {
457 			values[x][y] = (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1);
458 			rt_band_set_pixel(band, x, y, values[x][y], NULL);
459 		}
460 	}
461 
462 	gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
463 	CU_ASSERT(gddrv != NULL);
464 	CU_ASSERT_EQUAL(destroy, 0);
465 	CU_ASSERT(gdds != NULL);
466 	CU_ASSERT_EQUAL((uint32_t)GDALGetRasterXSize(gdds), width);
467 	CU_ASSERT_EQUAL((uint32_t)GDALGetRasterYSize(gdds), height);
468 
469 	rast = rt_raster_from_gdal_dataset(gdds);
470 	CU_ASSERT(rast != NULL);
471 	CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
472 
473 	band = rt_raster_get_band(rast, 0);
474 	CU_ASSERT(band != NULL);
475 
476 	for (x = 0; x < width; x++) {
477 		for (y = 0; y < height; y++) {
478 			rtn = rt_band_get_pixel(band, x, y, &value, NULL);
479  			CU_ASSERT_EQUAL(rtn, ES_NONE);
480 			CU_ASSERT_DOUBLE_EQUAL(value, values[x][y], DBL_EPSILON);
481 		}
482 	}
483 
484 	GDALClose(gdds);
485 	gdds = NULL;
486 	gddrv = NULL;
487 
488 	cu_free_raster(rast);
489 	cu_free_raster(raster);
490 
491 	raster = rt_raster_new(width, height);
492 	CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
493 
494 	pixtype = PT_8BSI;
495 	band = cu_add_band(raster, pixtype, 1, 0);
496 	CU_ASSERT(band != NULL);
497 
498 	v = -127;
499 	for (x = 0; x < width; x++) {
500 		for (y = 0; y < height; y++) {
501 			values[x][y] = v++;
502 			rt_band_set_pixel(band, x, y, values[x][y], NULL);
503 			if (v == 128)
504 				v = -127;
505 		}
506 	}
507 
508 	gdds = rt_raster_to_gdal_mem(raster, NULL, NULL, NULL, 0, &gddrv, &destroy);
509 	CU_ASSERT(gddrv != NULL);
510 	CU_ASSERT_EQUAL(destroy, 0);
511 	CU_ASSERT(gdds != NULL);
512 	CU_ASSERT_EQUAL((uint32_t)GDALGetRasterXSize(gdds), width);
513 	CU_ASSERT_EQUAL((uint32_t)GDALGetRasterYSize(gdds), height);
514 
515 	rast = rt_raster_from_gdal_dataset(gdds);
516 	CU_ASSERT(rast != NULL);
517 	CU_ASSERT_EQUAL(rt_raster_get_num_bands(rast), 1);
518 
519 	band = rt_raster_get_band(rast, 0);
520 	CU_ASSERT(band != NULL);
521 	CU_ASSERT_EQUAL(rt_band_get_pixtype(band), PT_16BSI);
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], 1.);
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 
test_gdal_warp()539 static void test_gdal_warp() {
540 	rt_pixtype pixtype = PT_64BF;
541 	rt_band band = NULL;
542 
543 	rt_raster raster;
544 	rt_raster rast;
545 	uint32_t x;
546 	uint32_t width = 100;
547 	uint32_t y;
548 	uint32_t height = 100;
549 	double value = 0;
550 
551 	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\"]]";
552 
553 	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]]";
554 
555 	raster = rt_raster_new(width, height);
556 	CU_ASSERT(raster != NULL); /* or we're out of virtual memory */
557 
558 	band = cu_add_band(raster, pixtype, 1, 0);
559 	CU_ASSERT(band != NULL);
560 
561 	rt_raster_set_offsets(raster, -500000, 600000);
562 	rt_raster_set_scale(raster, 1000, -1000);
563 
564 	for (x = 0; x < width; x++) {
565 		for (y = 0; y < height; y++) {
566 			rt_band_set_pixel(band, x, y, (((double) x * y) + (x + y) + (x + y * x)) / (x + y + 1), NULL);
567 		}
568 	}
569 
570 	rast = rt_raster_gdal_warp(
571 		raster,
572 		src_srs, dst_srs,
573 		NULL, NULL,
574 		NULL, NULL,
575 		NULL, NULL,
576 		NULL, NULL,
577 		NULL, NULL,
578 		GRA_NearestNeighbour, -1
579 	);
580 	CU_ASSERT(rast != NULL);
581 	CU_ASSERT_EQUAL(rt_raster_get_width(rast), 122);
582 	CU_ASSERT_EQUAL(rt_raster_get_height(rast), 116);
583 	CU_ASSERT_NOT_EQUAL(rt_raster_get_num_bands(rast), 0);
584 
585 	band = rt_raster_get_band(rast, 0);
586 	CU_ASSERT(band != NULL);
587 
588 	CU_ASSERT(rt_band_get_hasnodata_flag(band));
589 	rt_band_get_nodata(band, &value);
590 	CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
591 
592 	CU_ASSERT_EQUAL(rt_band_get_pixel(band, 0, 0, &value, NULL), ES_NONE);
593 	CU_ASSERT_DOUBLE_EQUAL(value, 0., DBL_EPSILON);
594 
595 	cu_free_raster(rast);
596 	cu_free_raster(raster);
597 }
598 
599 /* register tests */
600 void gdal_suite_setup(void);
gdal_suite_setup(void)601 void gdal_suite_setup(void)
602 {
603 	CU_pSuite suite = CU_add_suite("gdal", NULL, NULL);
604 	PG_ADD_TEST(suite, test_gdal_configured);
605 	PG_ADD_TEST(suite, test_gdal_drivers);
606 	PG_ADD_TEST(suite, test_gdal_rasterize);
607 	PG_ADD_TEST(suite, test_gdal_polygonize);
608 	PG_ADD_TEST(suite, test_raster_to_gdal);
609 	PG_ADD_TEST(suite, test_gdal_to_raster);
610 	PG_ADD_TEST(suite, test_gdal_warp);
611 }
612 
613