PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches
cu_raster_geometry.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-2013 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
27static void test_raster_envelope() {
28 rt_raster raster = NULL;
29 rt_envelope rtenv;
30
31 /* width = 0, height = 0 */
32 raster = rt_raster_new(0, 0);
33 CU_ASSERT(raster != NULL);
34
35 rt_raster_set_offsets(raster, 0.5, 0.5);
36 rt_raster_set_scale(raster, 1, -1);
37
38 CU_ASSERT_EQUAL(rt_raster_get_envelope(raster, &rtenv), ES_NONE);
39 CU_ASSERT_DOUBLE_EQUAL(rtenv.MinX, 0.5, DBL_EPSILON);
40 CU_ASSERT_DOUBLE_EQUAL(rtenv.MaxX, 0.5, DBL_EPSILON);
41 CU_ASSERT_DOUBLE_EQUAL(rtenv.MinY, 0.5, DBL_EPSILON);
42 CU_ASSERT_DOUBLE_EQUAL(rtenv.MaxY, 0.5, DBL_EPSILON);
43 cu_free_raster(raster);
44
45 /* width = 0 */
46 raster = rt_raster_new(0, 5);
47 CU_ASSERT(raster != NULL);
48
49 rt_raster_set_offsets(raster, 0.5, 0.5);
50 rt_raster_set_scale(raster, 1, -1);
51
52 CU_ASSERT_EQUAL(rt_raster_get_envelope(raster, &rtenv), ES_NONE);
53 CU_ASSERT_DOUBLE_EQUAL(rtenv.MinX, 0.5, DBL_EPSILON);
54 CU_ASSERT_DOUBLE_EQUAL(rtenv.MaxX, 0.5, DBL_EPSILON);
55 CU_ASSERT_DOUBLE_EQUAL(rtenv.MinY, -4.5, DBL_EPSILON);
56 CU_ASSERT_DOUBLE_EQUAL(rtenv.MaxY, 0.5, DBL_EPSILON);
57 cu_free_raster(raster);
58
59 /* height = 0 */
60 raster = rt_raster_new(5, 0);
61 CU_ASSERT(raster != NULL);
62
63 rt_raster_set_offsets(raster, 0.5, 0.5);
64 rt_raster_set_scale(raster, 1, -1);
65
66 CU_ASSERT_EQUAL(rt_raster_get_envelope(raster, &rtenv), ES_NONE);
67 CU_ASSERT_DOUBLE_EQUAL(rtenv.MinX, 0.5, DBL_EPSILON);
68 CU_ASSERT_DOUBLE_EQUAL(rtenv.MaxX, 5.5, DBL_EPSILON);
69 CU_ASSERT_DOUBLE_EQUAL(rtenv.MinY, 0.5, DBL_EPSILON);
70 CU_ASSERT_DOUBLE_EQUAL(rtenv.MaxY, 0.5, DBL_EPSILON);
71 cu_free_raster(raster);
72
73 /* normal raster */
74 raster = rt_raster_new(5, 5);
75 CU_ASSERT(raster != NULL);
76
77 rt_raster_set_offsets(raster, 0.5, 0.5);
78 rt_raster_set_scale(raster, 1, -1);
79
80 CU_ASSERT_EQUAL(rt_raster_get_envelope(raster, &rtenv), ES_NONE);
81 CU_ASSERT_DOUBLE_EQUAL(rtenv.MinX, 0.5, DBL_EPSILON);
82 CU_ASSERT_DOUBLE_EQUAL(rtenv.MaxX, 5.5, DBL_EPSILON);
83 CU_ASSERT_DOUBLE_EQUAL(rtenv.MinY, -4.5, DBL_EPSILON);
84 CU_ASSERT_DOUBLE_EQUAL(rtenv.MaxY, 0.5, DBL_EPSILON);
85 cu_free_raster(raster);
86}
87
88static void test_raster_envelope_geom() {
89 rt_raster raster = NULL;
90 LWGEOM *env = NULL;
91 LWPOLY *poly = NULL;
92 POINTARRAY *ring = NULL;
93 POINT4D pt;
94
95 /* NULL raster */
96 CU_ASSERT_EQUAL(rt_raster_get_envelope_geom(NULL, &env), ES_NONE);
97 CU_ASSERT(env == NULL);
98
99 /* width = 0, height = 0 */
100 raster = rt_raster_new(0, 0);
101 CU_ASSERT(raster != NULL);
102
103 CU_ASSERT_EQUAL(rt_raster_get_envelope_geom(raster, &env), ES_NONE);
104 CU_ASSERT_EQUAL(env->type, POINTTYPE);
105 lwgeom_free(env);
106 cu_free_raster(raster);
107
108 /* width = 0 */
109 raster = rt_raster_new(0, 256);
110 CU_ASSERT(raster != NULL);
111
112 CU_ASSERT_EQUAL(rt_raster_get_envelope_geom(raster, &env), ES_NONE);
113 CU_ASSERT_EQUAL(env->type, LINETYPE);
114 lwgeom_free(env);
115 cu_free_raster(raster);
116
117 /* height = 0 */
118 raster = rt_raster_new(256, 0);
119 CU_ASSERT(raster != NULL);
120
121 CU_ASSERT_EQUAL(rt_raster_get_envelope_geom(raster, &env), ES_NONE);
122 CU_ASSERT_EQUAL(env->type, LINETYPE);
123 lwgeom_free(env);
124 cu_free_raster(raster);
125
126 /* normal raster */
127 raster = rt_raster_new(5, 5);
128 CU_ASSERT(raster != NULL);
129
130 rt_raster_set_offsets(raster, 0.5, 0.5);
131 rt_raster_set_scale(raster, 1, -1);
132
133 CU_ASSERT_EQUAL(rt_raster_get_envelope_geom(raster, &env), ES_NONE);
134 poly = lwgeom_as_lwpoly(env);
135 CU_ASSERT_EQUAL(poly->srid, rt_raster_get_srid(raster));
136 CU_ASSERT_EQUAL(poly->nrings, 1);
137
138 ring = poly->rings[0];
139 CU_ASSERT(ring != NULL);
140 CU_ASSERT_EQUAL(ring->npoints, 5);
141
142 getPoint4d_p(ring, 0, &pt);
143 CU_ASSERT_DOUBLE_EQUAL(pt.x, 0.5, DBL_EPSILON);
144 CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.5, DBL_EPSILON);
145
146 getPoint4d_p(ring, 1, &pt);
147 CU_ASSERT_DOUBLE_EQUAL(pt.x, 5.5, DBL_EPSILON);
148 CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.5, DBL_EPSILON);
149
150 getPoint4d_p(ring, 2, &pt);
151 CU_ASSERT_DOUBLE_EQUAL(pt.x, 5.5, DBL_EPSILON);
152 CU_ASSERT_DOUBLE_EQUAL(pt.y, -4.5, DBL_EPSILON);
153
154 getPoint4d_p(ring, 3, &pt);
155 CU_ASSERT_DOUBLE_EQUAL(pt.x, 0.5, DBL_EPSILON);
156 CU_ASSERT_DOUBLE_EQUAL(pt.y, -4.5, DBL_EPSILON);
157
158 getPoint4d_p(ring, 4, &pt);
159 CU_ASSERT_DOUBLE_EQUAL(pt.x, 0.5, DBL_EPSILON);
160 CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.5, DBL_EPSILON);
161
162 lwgeom_free(env);
163 cu_free_raster(raster);
164}
165
166static void test_raster_convex_hull() {
167 rt_raster raster = NULL;
168 LWGEOM *hull = NULL;
169 LWPOLY *poly = NULL;
170 POINTARRAY *ring = NULL;
171 POINT4D pt;
172
173 /* NULL raster */
174 CU_ASSERT_EQUAL(rt_raster_get_convex_hull(NULL, &hull), ES_NONE);
175 CU_ASSERT(hull == NULL);
176
177 /* width = 0, height = 0 */
178 raster = rt_raster_new(0, 0);
179 CU_ASSERT(raster != NULL);
180
181 CU_ASSERT_EQUAL(rt_raster_get_convex_hull(raster, &hull), ES_NONE);
182 CU_ASSERT_EQUAL(hull->type, POINTTYPE);
183 lwgeom_free(hull);
184 cu_free_raster(raster);
185
186 /* width = 0 */
187 raster = rt_raster_new(0, 256);
188 CU_ASSERT(raster != NULL);
189
190 CU_ASSERT_EQUAL(rt_raster_get_convex_hull(raster, &hull), ES_NONE);
191 CU_ASSERT_EQUAL(hull->type, LINETYPE);
192 lwgeom_free(hull);
193 cu_free_raster(raster);
194
195 /* height = 0 */
196 raster = rt_raster_new(256, 0);
197 CU_ASSERT(raster != NULL);
198
199 CU_ASSERT_EQUAL(rt_raster_get_convex_hull(raster, &hull), ES_NONE);
200 CU_ASSERT_EQUAL(hull->type, LINETYPE);
201 lwgeom_free(hull);
202 cu_free_raster(raster);
203
204 /* normal raster */
205 raster = rt_raster_new(256, 256);
206 CU_ASSERT(raster != NULL);
207
208 rt_raster_set_offsets(raster, 0.5, 0.5);
209 rt_raster_set_scale(raster, 1, 1);
210 rt_raster_set_skews(raster, 4, 5);
211
212 CU_ASSERT_EQUAL(rt_raster_get_convex_hull(raster, &hull), ES_NONE);
213 poly = lwgeom_as_lwpoly(hull);
214 CU_ASSERT_EQUAL(poly->srid, rt_raster_get_srid(raster));
215 CU_ASSERT_EQUAL(poly->nrings, 1);
216
217 ring = poly->rings[0];
218 CU_ASSERT(ring != NULL);
219 CU_ASSERT_EQUAL(ring->npoints, 5);
220
221 getPoint4d_p(ring, 0, &pt);
222 CU_ASSERT_DOUBLE_EQUAL(pt.x, 0.5, DBL_EPSILON);
223 CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.5, DBL_EPSILON);
224
225 getPoint4d_p(ring, 1, &pt);
226 CU_ASSERT_DOUBLE_EQUAL(pt.x, 256.5, DBL_EPSILON);
227 CU_ASSERT_DOUBLE_EQUAL(pt.y, 1280.5, DBL_EPSILON);
228
229 getPoint4d_p(ring, 2, &pt);
230 CU_ASSERT_DOUBLE_EQUAL(pt.x, 1280.5, DBL_EPSILON);
231 CU_ASSERT_DOUBLE_EQUAL(pt.y, 1536.5, DBL_EPSILON);
232
233 getPoint4d_p(ring, 3, &pt);
234 CU_ASSERT_DOUBLE_EQUAL(pt.x, 1024.5, DBL_EPSILON);
235 CU_ASSERT_DOUBLE_EQUAL(pt.y, 256.5, DBL_EPSILON);
236
237 getPoint4d_p(ring, 4, &pt);
238 CU_ASSERT_DOUBLE_EQUAL(pt.x, 0.5, DBL_EPSILON);
239 CU_ASSERT_DOUBLE_EQUAL(pt.y, 0.5, DBL_EPSILON);
240
241 lwgeom_free(hull);
242 cu_free_raster(raster);
243}
244
245static char *
246lwgeom_to_text(const LWGEOM *lwgeom) {
247 char *wkt;
248 size_t wkt_size;
249
250 wkt = lwgeom_to_wkt(lwgeom, WKT_ISO, DBL_DIG, &wkt_size);
251
252 return wkt;
253}
254
255static void test_raster_surface() {
256 rt_raster rast;
257 rt_band band;
258 const int maxX = 5;
259 const int maxY = 5;
260 int x, y;
261 LWMPOLY *mpoly = NULL;
262 int err;
263 LWGEOM *gexpected, *gobserved;
264
265 rast = rt_raster_new(maxX, maxY);
266 CU_ASSERT(rast != NULL);
267
268 rt_raster_set_offsets(rast, 0, 0);
269 rt_raster_set_scale(rast, 1, -1);
270
271 band = cu_add_band(rast, PT_32BUI, 1, 0);
272 CU_ASSERT(band != NULL);
273
274 for (y = 0; y < maxY; y++) {
275 for (x = 0; x < maxX; x++) {
276 rt_band_set_pixel(band, x, y, 1, NULL);
277 }
278 }
279
280 err = rt_raster_surface(rast, 0, &mpoly);
281 CU_ASSERT_EQUAL(err, ES_NONE);
282 CU_ASSERT(mpoly != NULL);
283 //wkt = lwgeom_to_text(lwmpoly_as_lwgeom(mpoly));
284 gobserved = (LWGEOM *)lwmpoly_as_lwgeom(mpoly);
285 gexpected = lwgeom_from_wkt("MULTIPOLYGON(((0 0,0 -5,5 -5,5 0,0 0)))", LW_PARSER_CHECK_NONE);
286 CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
287 lwmpoly_free(mpoly);
288 mpoly = NULL;
289
290 /* 0,0 NODATA */
291 rt_band_set_pixel(band, 0, 0, 0, NULL);
292
293 err = rt_raster_surface(rast, 0, &mpoly);
294 CU_ASSERT_EQUAL(err, ES_NONE);
295 CU_ASSERT(mpoly != NULL);
296 gobserved = (LWGEOM *)lwmpoly_as_lwgeom(mpoly);
297 gexpected = lwgeom_from_wkt("MULTIPOLYGON(((1 0,1 -1,0 -1,0 -5,4 -5,5 -5,5 0,1 0)))", LW_PARSER_CHECK_NONE);
298 CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
299 lwmpoly_free(mpoly);
300 mpoly = NULL;
301
302 /* plus 1,1 NODATA */
303 rt_band_set_pixel(band, 1, 1, 0, NULL);
304
305 err = rt_raster_surface(rast, 0, &mpoly);
306 CU_ASSERT_EQUAL(err, ES_NONE);
307 CU_ASSERT(mpoly != NULL);
308 gobserved = (LWGEOM *)lwmpoly_as_lwgeom(mpoly);
309 gexpected = lwgeom_from_wkt("MULTIPOLYGON(((1 0,1 -1,0 -1,0 -5,4 -5,5 -5,5 0,1 0),(1 -1,1 -2,2 -2,2 -1,1 -1)))",
310 LW_PARSER_CHECK_NONE);
311 CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
312 lwmpoly_free(mpoly);
313 mpoly = NULL;
314
315 /* plus 2,2 NODATA */
316 rt_band_set_pixel(band, 2, 2, 0, NULL);
317
318 err = rt_raster_surface(rast, 0, &mpoly);
319 CU_ASSERT_EQUAL(err, ES_NONE);
320 CU_ASSERT(mpoly != NULL);
321 gobserved = (LWGEOM *)lwmpoly_as_lwgeom(mpoly);
322 gexpected = lwgeom_from_wkt(
323 "MULTIPOLYGON(((1 -1,1 0,5 0,5 -5,4 -5,0 -5,0 -1,1 -1),(1 -1,1 -2,2 -2,2 -1,1 -1),(2 -2,2 -3,3 -3,3 -2,2 -2)))",
324 LW_PARSER_CHECK_NONE);
325 CU_ASSERT_DOUBLE_EQUAL(lwgeom_perimeter(gobserved), lwgeom_perimeter(gexpected), FLT_EPSILON);
326 lwmpoly_free(mpoly);
327 mpoly = NULL;
328
329 /* plus 3,3 NODATA */
330 rt_band_set_pixel(band, 3, 3, 0, NULL);
331
332 err = rt_raster_surface(rast, 0, &mpoly);
333 CU_ASSERT_EQUAL(err, ES_NONE);
334 CU_ASSERT(mpoly != NULL);
335 gobserved = (LWGEOM *)lwmpoly_as_lwgeom(mpoly);
336 gexpected = lwgeom_from_wkt(
337 "MULTIPOLYGON(((1 -1,1 0,5 0,5 -5,4 -5,0 -5,0 -1,1 -1),(1 -1,1 -2,2 -2,2 -1,1 -1),(2 -2,2 -3,3 -3,3 -2,2 -2),(3 -3,3 -4,4 -4,4 -3,3 -3)))",
338 LW_PARSER_CHECK_NONE);
339 CU_ASSERT_DOUBLE_EQUAL(lwgeom_perimeter(gobserved), lwgeom_perimeter(gexpected), FLT_EPSILON);
340 lwmpoly_free(mpoly);
341 mpoly = NULL;
342
343 /* plus 4,4 NODATA */
344 rt_band_set_pixel(band, 4, 4, 0, NULL);
345
346 err = rt_raster_surface(rast, 0, &mpoly);
347 CU_ASSERT_EQUAL(err, ES_NONE);
348 CU_ASSERT(mpoly != NULL);
349 gobserved = (LWGEOM *)lwmpoly_as_lwgeom(mpoly);
350 gexpected = lwgeom_from_wkt(
351 "MULTIPOLYGON(((4 -4,4 -5,0 -5,0 -1,1 -1,1 -2,2 -2,2 -3,3 -3,3 -4,4 -4)),((1 -1,1 0,5 0,5 -4,4 -4,4 -3,3 -3,3 -2,2 -2,2 -1,1 -1)))",
352 LW_PARSER_CHECK_NONE);
353 CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
354 lwmpoly_free(mpoly);
355 mpoly = NULL;
356
357 /* a whole lot more NODATA */
358 rt_band_set_pixel(band, 4, 0, 0, NULL);
359 rt_band_set_pixel(band, 3, 1, 0, NULL);
360 rt_band_set_pixel(band, 1, 3, 0, NULL);
361 rt_band_set_pixel(band, 0, 4, 0, NULL);
362
363 err = rt_raster_surface(rast, 0, &mpoly);
364 CU_ASSERT_EQUAL(err, ES_NONE);
365 CU_ASSERT(mpoly != NULL);
366 gobserved = (LWGEOM *)lwmpoly_as_lwgeom(mpoly);
367 gexpected = lwgeom_from_wkt(
368 "MULTIPOLYGON(((1 -4,2 -4,2 -3,3 -3,3 -4,4 -4,4 -5,3 -5,1 -5,1 -4)),((1 -4,0 -4,0 -1,1 -1,1 -2,2 -2,2 -3,1 -3,1 -4)),((3 -2,4 -2,4 -1,5 -1,5 -4,4 -4,4 -3,3 -3,3 -2)),((3 -2,2 -2,2 -1,1 -1,1 0,4 0,4 -1,3 -1,3 -2)))",
369 LW_PARSER_CHECK_NONE);
370 CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(gobserved), lwgeom_area(gexpected), FLT_EPSILON);
371 lwmpoly_free(mpoly);
372 mpoly = NULL;
373
374 cu_free_raster(rast);
375}
376
377static void test_raster_perimeter() {
378 rt_raster rast;
379 rt_band band;
380 const int maxX = 5;
381 const int maxY = 5;
382 int x, y;
383 char *wkt = NULL;
384 LWGEOM *geom = NULL;
385 int err;
386
387 rast = rt_raster_new(maxX, maxY);
388 CU_ASSERT(rast != NULL);
389
390 rt_raster_set_offsets(rast, 0, 0);
391 rt_raster_set_scale(rast, 1, -1);
392
393 band = cu_add_band(rast, PT_32BUI, 1, 0);
394 CU_ASSERT(band != NULL);
395
396 for (y = 0; y < maxY; y++) {
397 for (x = 0; x < maxX; x++) {
398 rt_band_set_pixel(band, x, y, 1, NULL);
399 }
400 }
401
402 err = rt_raster_get_perimeter(rast, -1, &geom);
403 CU_ASSERT_EQUAL(err, ES_NONE);
404 CU_ASSERT(geom != NULL);
405 wkt = lwgeom_to_text(geom);
406 CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,5 0,5 -5,0 -5,0 0))");
407 rtdealloc(wkt);
408 lwgeom_free(geom);
409 geom = NULL;
410
411 /* row 0 is NODATA */
412 rt_band_set_pixel(band, 0, 0, 0, NULL);
413 rt_band_set_pixel(band, 1, 0, 0, NULL);
414 rt_band_set_pixel(band, 2, 0, 0, NULL);
415 rt_band_set_pixel(band, 3, 0, 0, NULL);
416 rt_band_set_pixel(band, 4, 0, 0, NULL);
417
418 err = rt_raster_get_perimeter(rast, 0, &geom);
419 CU_ASSERT_EQUAL(err, ES_NONE);
420 CU_ASSERT(geom != NULL);
421 wkt = lwgeom_to_text(geom);
422 CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 -1,5 -1,5 -5,0 -5,0 -1))");
423 rtdealloc(wkt);
424 lwgeom_free(geom);
425 geom = NULL;
426
427 /* column 4 is NODATA */
428 /* pixel 4, 0 already set to NODATA */
429 rt_band_set_pixel(band, 4, 1, 0, NULL);
430 rt_band_set_pixel(band, 4, 2, 0, NULL);
431 rt_band_set_pixel(band, 4, 3, 0, NULL);
432 rt_band_set_pixel(band, 4, 4, 0, NULL);
433
434 err = rt_raster_get_perimeter(rast, 0, &geom);
435 CU_ASSERT_EQUAL(err, ES_NONE);
436 CU_ASSERT(geom != NULL);
437 wkt = lwgeom_to_text(geom);
438 CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 -1,4 -1,4 -5,0 -5,0 -1))");
439 rtdealloc(wkt);
440 lwgeom_free(geom);
441 geom = NULL;
442
443 /* row 4 is NODATA */
444 rt_band_set_pixel(band, 0, 4, 0, NULL);
445 rt_band_set_pixel(band, 1, 4, 0, NULL);
446 rt_band_set_pixel(band, 2, 4, 0, NULL);
447 rt_band_set_pixel(band, 3, 4, 0, NULL);
448 /* pixel 4, 4 already set to NODATA */
449
450 err = rt_raster_get_perimeter(rast, 0, &geom);
451 CU_ASSERT_EQUAL(err, ES_NONE);
452 CU_ASSERT(geom != NULL);
453 wkt = lwgeom_to_text(geom);
454 CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 -1,4 -1,4 -4,0 -4,0 -1))");
455 rtdealloc(wkt);
456 lwgeom_free(geom);
457 geom = NULL;
458
459 /* column 0 is NODATA */
460 /* pixel 0, 0 already set to NODATA*/
461 rt_band_set_pixel(band, 0, 1, 0, NULL);
462 rt_band_set_pixel(band, 0, 2, 0, NULL);
463 rt_band_set_pixel(band, 0, 3, 0, NULL);
464 /* pixel 0, 4 already set to NODATA*/
465
466 err = rt_raster_get_perimeter(rast, 0, &geom);
467 CU_ASSERT_EQUAL(err, ES_NONE);
468 CU_ASSERT(geom != NULL);
469 wkt = lwgeom_to_text(geom);
470 CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((1 -1,4 -1,4 -4,1 -4,1 -1))");
471 rtdealloc(wkt);
472 lwgeom_free(geom);
473 geom = NULL;
474
475 /* columns 1 and 3 are NODATA */
476 /* pixel 1, 0 already set to NODATA */
477 rt_band_set_pixel(band, 1, 1, 0, NULL);
478 rt_band_set_pixel(band, 1, 2, 0, NULL);
479 rt_band_set_pixel(band, 1, 3, 0, NULL);
480 /* pixel 1, 4 already set to NODATA */
481 /* pixel 3, 0 already set to NODATA */
482 rt_band_set_pixel(band, 3, 1, 0, NULL);
483 rt_band_set_pixel(band, 3, 2, 0, NULL);
484 rt_band_set_pixel(band, 3, 3, 0, NULL);
485 /* pixel 3, 4 already set to NODATA */
486
487 err = rt_raster_get_perimeter(rast, 0, &geom);
488 CU_ASSERT_EQUAL(err, ES_NONE);
489 CU_ASSERT(geom != NULL);
490 wkt = lwgeom_to_text(geom);
491 CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((2 -1,3 -1,3 -4,2 -4,2 -1))");
492 rtdealloc(wkt);
493 lwgeom_free(geom);
494 geom = NULL;
495
496 /* more pixels are NODATA */
497 rt_band_set_pixel(band, 2, 1, 0, NULL);
498 rt_band_set_pixel(band, 2, 3, 0, NULL);
499
500 err = rt_raster_get_perimeter(rast, 0, &geom);
501 CU_ASSERT_EQUAL(err, ES_NONE);
502 CU_ASSERT(geom != NULL);
503 wkt = lwgeom_to_text(geom);
504 CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((2 -2,3 -2,3 -3,2 -3,2 -2))");
505 rtdealloc(wkt);
506 lwgeom_free(geom);
507 geom = NULL;
508
509 /* second band */
510 band = cu_add_band(rast, PT_32BUI, 1, 0);
511 CU_ASSERT(band != NULL);
512
513 for (y = 0; y < maxY; y++) {
514 for (x = 0; x < maxX; x++) {
515 rt_band_set_pixel(band, x, y, 1, NULL);
516 }
517 }
518
519 err = rt_raster_get_perimeter(rast, -1, &geom);
520 CU_ASSERT_EQUAL(err, ES_NONE);
521 CU_ASSERT(geom != NULL);
522 wkt = lwgeom_to_text(geom);
523 CU_ASSERT_STRING_EQUAL(wkt, "POLYGON((0 0,5 0,5 -5,0 -5,0 0))");
524 rtdealloc(wkt);
525 lwgeom_free(geom);
526 geom = NULL;
527
528 cu_free_raster(rast);
529}
530
531static void test_raster_pixel_as_polygon() {
532 rt_raster rast;
533 rt_band band;
534 uint32_t x, y;
535 const uint32_t maxX = 10;
536 const uint32_t maxY = 10;
537 LWPOLY *poly = NULL;
538
539 rast = rt_raster_new(maxX, maxY);
540 CU_ASSERT(rast != NULL);
541
542 band = cu_add_band(rast, PT_32BUI, 1, 0);
543 CU_ASSERT(band != NULL);
544
545 for (x = 0; x < maxX; x++) {
546 for (y = 0; y < maxY; y++) {
547 rt_band_set_pixel(band, x, y, 1, NULL);
548 }
549 }
550
551 rt_band_set_pixel(band, 0, 0, 0, NULL);
552 rt_band_set_pixel(band, 3, 0, 0, NULL);
553 rt_band_set_pixel(band, 6, 0, 0, NULL);
554 rt_band_set_pixel(band, 9, 0, 0, NULL);
555 rt_band_set_pixel(band, 1, 2, 0, NULL);
556 rt_band_set_pixel(band, 4, 2, 0, NULL);
557 rt_band_set_pixel(band, 7, 2, 0, NULL);
558 rt_band_set_pixel(band, 2, 4, 0, NULL);
559 rt_band_set_pixel(band, 5, 4, 0, NULL);
560 rt_band_set_pixel(band, 8, 4, 0, NULL);
561 rt_band_set_pixel(band, 0, 6, 0, NULL);
562 rt_band_set_pixel(band, 3, 6, 0, NULL);
563 rt_band_set_pixel(band, 6, 6, 0, NULL);
564 rt_band_set_pixel(band, 9, 6, 0, NULL);
565 rt_band_set_pixel(band, 1, 8, 0, NULL);
566 rt_band_set_pixel(band, 4, 8, 0, NULL);
567 rt_band_set_pixel(band, 7, 8, 0, NULL);
568
569 poly = rt_raster_pixel_as_polygon(rast, 1, 1);
570 CU_ASSERT(poly != NULL);
571 lwpoly_free(poly);
572
573 cu_free_raster(rast);
574}
575
576
577
578
579static void test_raster_get_pixel_bilinear() {
580 uint32_t width = 2;
581 uint32_t height = 2;
582 double ul_x = 0.0;
583 double ul_y = 0.0;
584 double scale_x = 1;
585 double scale_y = 1;
586 int err;
587
588 double xr, yr;
589 double igt[6];
590
591 rt_raster rast = rt_raster_new(width, height);
592 rt_raster_set_offsets(rast, ul_x, ul_y);
593 rt_raster_set_scale(rast, scale_x, scale_y);
594
595 double xw = 1.5, yw = 0.5;
596
597 rt_raster_generate_new_band(
598 rast, // rt_raster raster,
599 PT_64BF, // rt_pixtype pixtype,
600 1.0, // double initialvalue,
601 1, // uint32_t hasnodata,
602 -99.0, // double nodatavalue,
603 0 // int index
604 );
605
606 rt_raster_geopoint_to_rasterpoint(
607 rast,
608 xw, yw,
609 &xr, &yr, igt);
610
611 printf("xw = %g, yw = %g, xr = %g, yr = %g\n", xw, yw, xr, yr);
612
613 // err = rt_raster_cell_to_geopoint(
614 // rast,
615 // xr, yr,
616 // &xw, &yw, igt);
617
618 // printf("xw = %g, yw = %g, xr = %g, yr = %g\n", xw, yw, xr, yr);
619
620 rt_band band = rt_raster_get_band(rast, 0);
621 rt_band_set_pixel(band, 0, 0, 10.0, NULL);
622 rt_band_set_pixel(band, 0, 1, 10.0, NULL);
623 rt_band_set_pixel(band, 1, 0, 20.0, NULL);
624 rt_band_set_pixel(band, 1, 1, 40.0, NULL);
625
626
627 double value = 0.0;
628 int nodata = 0;
629 err = rt_band_get_pixel_bilinear(
630 band,
631 xw, yw, // double xw, double yw,
632 &value, &nodata // double *r_value, int *r_nodata)
633 );
634
635 CU_ASSERT_EQUAL(err, ES_NONE);
636 // printf("xw = %g, yw = %g, value = %g, nodata = %d\n", xr, yr, value, nodata);
637
638
639}
640
641/* register tests */
642void raster_geometry_suite_setup(void);
643void raster_geometry_suite_setup(void)
644{
645 CU_pSuite suite = CU_add_suite("raster_geometry", NULL, NULL);
646 PG_ADD_TEST(suite, test_raster_envelope);
647 PG_ADD_TEST(suite, test_raster_envelope_geom);
648 PG_ADD_TEST(suite, test_raster_convex_hull);
649 PG_ADD_TEST(suite, test_raster_surface);
650 PG_ADD_TEST(suite, test_raster_perimeter);
651 PG_ADD_TEST(suite, test_raster_pixel_as_polygon);
652 PG_ADD_TEST(suite, test_raster_get_pixel_bilinear);
653}
654
test_raster_surface
static void test_raster_surface()
Definition cu_raster_geometry.c:255
test_raster_envelope
static void test_raster_envelope()
Definition cu_raster_geometry.c:27
test_raster_convex_hull
static void test_raster_convex_hull()
Definition cu_raster_geometry.c:166
test_raster_envelope_geom
static void test_raster_envelope_geom()
Definition cu_raster_geometry.c:88
lwgeom_to_text
static char * lwgeom_to_text(const LWGEOM *lwgeom)
Definition cu_raster_geometry.c:246
test_raster_pixel_as_polygon
static void test_raster_pixel_as_polygon()
Definition cu_raster_geometry.c:531
test_raster_get_pixel_bilinear
static void test_raster_get_pixel_bilinear()
Definition cu_raster_geometry.c:579
test_raster_perimeter
static void test_raster_perimeter()
Definition cu_raster_geometry.c:377
raster_geometry_suite_setup
void raster_geometry_suite_setup(void)
Definition cu_raster_geometry.c:643
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:1246
LW_PARSER_CHECK_NONE
#define LW_PARSER_CHECK_NONE
Definition liblwgeom.h:2149
lwmpoly_free
void lwmpoly_free(LWMPOLY *mpoly)
Definition lwmpoly.c:53
LINETYPE
#define LINETYPE
Definition liblwgeom.h:103
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:708
POINTTYPE
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition liblwgeom.h:102
lwgeom_as_lwpoly
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition lwgeom.c:243
lwgeom_area
double lwgeom_area(const LWGEOM *geom)
Definition lwgeom.c:1999
WKT_ISO
#define WKT_ISO
Definition liblwgeom.h:2219
getPoint4d_p
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition lwgeom_api.c:125
lwpoly_free
void lwpoly_free(LWPOLY *poly)
Definition lwpoly.c:175
lwgeom_perimeter
double lwgeom_perimeter(const LWGEOM *geom)
Definition lwgeom.c:2022
lwmpoly_as_lwgeom
LWGEOM * lwmpoly_as_lwgeom(const LWMPOLY *obj)
Definition lwgeom.c:322
lwgeom_from_wkt
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition lwin_wkt.c:940
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition rt_raster.c:360
rt_raster_generate_new_band
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
Definition rt_raster.c:489
rt_raster_pixel_as_polygon
LWPOLY * rt_raster_pixel_as_polygon(rt_raster raster, int x, int y)
Get a raster pixel as a polygon.
Definition rt_geometry.c:608
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
PT_32BUI
@ PT_32BUI
Definition librtcore.h:197
PT_64BF
@ PT_64BF
Definition librtcore.h:200
rt_raster_set_skews
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition rt_raster.c:172
rt_raster_surface
rt_errorstate rt_raster_surface(rt_raster raster, int nband, LWMPOLY **surface)
Get a raster as a surface (multipolygon).
Definition rt_geometry.c:356
rt_raster_geopoint_to_rasterpoint
rt_errorstate rt_raster_geopoint_to_rasterpoint(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw,yw map point to a xr,yr raster point.
Definition rt_raster.c:730
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_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_get_convex_hull
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition rt_geometry.c:839
rtdealloc
void rtdealloc(void *mem)
Definition rt_context.c:206
rt_band_get_pixel_bilinear
rt_errorstate rt_band_get_pixel_bilinear(rt_band band, double xr, double yr, double *r_value, int *r_nodata)
Retrieve a point value from the raster using a world coordinate and bilinear interpolation.
Definition rt_band.c:1435
rt_raster_get_envelope_geom
rt_errorstate rt_raster_get_envelope_geom(rt_raster raster, LWGEOM **env)
Get raster's envelope as a geometry.
Definition rt_geometry.c:706
rt_raster_get_perimeter
rt_errorstate rt_raster_get_perimeter(rt_raster raster, int nband, LWGEOM **perimeter)
Get raster perimeter.
Definition rt_geometry.c:183
rt_raster_get_envelope
rt_errorstate rt_raster_get_envelope(rt_raster raster, rt_envelope *env)
Get raster's envelope.
Definition rt_raster.c:782
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_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
LWGEOM::type
uint8_t type
Definition liblwgeom.h:462
LWPOLY::rings
POINTARRAY ** rings
Definition liblwgeom.h:519
LWPOLY::nrings
uint32_t nrings
Definition liblwgeom.h:524
LWPOLY::srid
int32_t srid
Definition liblwgeom.h:520
POINT4D::x
double x
Definition liblwgeom.h:414
POINT4D::y
double y
Definition liblwgeom.h:414
POINTARRAY::npoints
uint32_t npoints
Definition liblwgeom.h:427
rt_band_t
Definition librtcore.h:2515
rt_envelope::MinX
double MinX
Definition librtcore.h:167
rt_envelope::MaxX
double MaxX
Definition librtcore.h:168
rt_envelope::MinY
double MinY
Definition librtcore.h:169
rt_envelope::MaxY
double MaxY
Definition librtcore.h:170
rt_envelope
Definition librtcore.h:166
rt_raster_t
Definition librtcore.h:2486