1 /*
2 *
3 * WKTRaster - Raster Types for PostGIS
4 * http://trac.osgeo.org/postgis/wiki/WKTRaster
5 *
6 * Copyright (C) 2013 Bborie Park <dustymugs@gmail.com>
7 * Copyright (C) 2011-2013 Regents of the University of California
8 * <bkpark@ucdavis.edu>
9 * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
10 * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
11 * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
12 * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
13 * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
14 *
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * This program is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with this program; if not, write to the Free Software Foundation,
27 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28 *
29 */
30
31 #include "../../postgis_config.h"
32 /* #define POSTGIS_DEBUG_LEVEL 4 */
33
34 #include "librtcore.h"
35 #include "librtcore_internal.h"
36
37 /******************************************************************************
38 * rt_raster_gdal_warp()
39 ******************************************************************************/
40
41 typedef struct _rti_warp_arg_t* _rti_warp_arg;
42 struct _rti_warp_arg_t {
43
44 struct {
45 GDALDriverH drv;
46 GDALDatasetH ds;
47 char *srs;
48 int destroy_drv;
49 } src, dst;
50
51 GDALWarpOptions *wopts;
52
53 struct {
54 struct {
55 char **item;
56 int len;
57 } option;
58
59 struct {
60 void *transform;
61 void *imgproj;
62 void *approx;
63 } arg;
64
65 GDALTransformerFunc func;
66 } transform;
67
68 };
69
70 static _rti_warp_arg
_rti_warp_arg_init()71 _rti_warp_arg_init() {
72 _rti_warp_arg arg = NULL;
73
74 arg = rtalloc(sizeof(struct _rti_warp_arg_t));
75 if (arg == NULL) {
76 rterror("_rti_warp_arg_init: Could not allocate memory for _rti_warp_arg");
77 return NULL;
78 }
79
80 arg->src.drv = NULL;
81 arg->src.destroy_drv = 0;
82 arg->src.ds = NULL;
83 arg->src.srs = NULL;
84
85 arg->dst.drv = NULL;
86 arg->dst.destroy_drv = 0;
87 arg->dst.ds = NULL;
88 arg->dst.srs = NULL;
89
90 arg->wopts = NULL;
91
92 arg->transform.option.item = NULL;
93 arg->transform.option.len = 0;
94
95 arg->transform.arg.transform = NULL;
96 arg->transform.arg.imgproj = NULL;
97 arg->transform.arg.approx = NULL;
98
99 arg->transform.func = NULL;
100
101 return arg;
102 }
103
104 static void
_rti_warp_arg_destroy(_rti_warp_arg arg)105 _rti_warp_arg_destroy(_rti_warp_arg arg) {
106 int i = 0;
107
108 if (arg->dst.ds != NULL)
109 GDALClose(arg->dst.ds);
110 if (arg->dst.srs != NULL)
111 CPLFree(arg->dst.srs);
112
113 if (arg->dst.drv != NULL && arg->dst.destroy_drv) {
114 GDALDeregisterDriver(arg->dst.drv);
115 GDALDestroyDriver(arg->dst.drv);
116 }
117
118 if (arg->src.ds != NULL)
119 GDALClose(arg->src.ds);
120 if (arg->src.srs != NULL)
121 CPLFree(arg->src.srs);
122
123 if (arg->src.drv != NULL && arg->src.destroy_drv) {
124 GDALDeregisterDriver(arg->src.drv);
125 GDALDestroyDriver(arg->src.drv);
126 }
127
128 if (arg->transform.func == GDALApproxTransform) {
129 if (arg->transform.arg.imgproj != NULL)
130 GDALDestroyGenImgProjTransformer(arg->transform.arg.imgproj);
131 }
132
133 if (arg->wopts != NULL)
134 GDALDestroyWarpOptions(arg->wopts);
135
136 if (arg->transform.option.len > 0 && arg->transform.option.item != NULL) {
137 for (i = 0; i < arg->transform.option.len; i++) {
138 if (arg->transform.option.item[i] != NULL)
139 rtdealloc(arg->transform.option.item[i]);
140 }
141 rtdealloc(arg->transform.option.item);
142 }
143
144 rtdealloc(arg);
145 arg = NULL;
146 }
147
148 /**
149 * Return a warped raster using GDAL Warp API
150 *
151 * @param raster : raster to transform
152 * @param src_srs : the raster's coordinate system in OGC WKT
153 * @param dst_srs : the warped raster's coordinate system in OGC WKT
154 * @param scale_x : the x size of pixels of the warped raster's pixels in
155 * units of dst_srs
156 * @param scale_y : the y size of pixels of the warped raster's pixels in
157 * units of dst_srs
158 * @param width : the number of columns of the warped raster. note that
159 * width/height CANNOT be used with scale_x/scale_y
160 * @param height : the number of rows of the warped raster. note that
161 * width/height CANNOT be used with scale_x/scale_y
162 * @param ul_xw : the X value of upper-left corner of the warped raster in
163 * units of dst_srs
164 * @param ul_yw : the Y value of upper-left corner of the warped raster in
165 * units of dst_srs
166 * @param grid_xw : the X value of point on a grid to align warped raster
167 * to in units of dst_srs
168 * @param grid_yw : the Y value of point on a grid to align warped raster
169 * to in units of dst_srs
170 * @param skew_x : the X skew of the warped raster in units of dst_srs
171 * @param skew_y : the Y skew of the warped raster in units of dst_srs
172 * @param resample_alg : the resampling algorithm
173 * @param max_err : maximum error measured in input pixels permitted
174 * (0.0 for exact calculations)
175 *
176 * @return the warped raster or NULL
177 */
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)178 rt_raster rt_raster_gdal_warp(
179 rt_raster raster,
180 const char *src_srs, const char *dst_srs,
181 double *scale_x, double *scale_y,
182 int *width, int *height,
183 double *ul_xw, double *ul_yw,
184 double *grid_xw, double *grid_yw,
185 double *skew_x, double *skew_y,
186 GDALResampleAlg resample_alg, double max_err
187 ) {
188 CPLErr cplerr;
189 char *dst_options[] = {"SUBCLASS=VRTWarpedDataset", NULL};
190 _rti_warp_arg arg = NULL;
191
192 int hasnodata = 0;
193
194 GDALRasterBandH band;
195 rt_band rtband = NULL;
196 rt_pixtype pt = PT_END;
197 GDALDataType gdal_pt = GDT_Unknown;
198 double nodata = 0.0;
199
200 double _gt[6] = {0};
201 double dst_extent[4];
202 rt_envelope extent;
203
204 int _dim[2] = {0};
205 double _skew[2] = {0};
206 double _scale[2] = {0};
207 int ul_user = 0;
208
209 rt_raster rast = NULL;
210 int i = 0;
211 int numBands = 0;
212
213 /* flag indicating that the spatial info is being substituted */
214 int subspatial = 0;
215
216 RASTER_DEBUG(3, "starting");
217
218 assert(NULL != raster);
219
220 /* internal variables */
221 arg = _rti_warp_arg_init();
222 if (arg == NULL) {
223 rterror("rt_raster_gdal_warp: Could not initialize internal variables");
224 return NULL;
225 }
226
227 /*
228 max_err must be gte zero
229
230 the value 0.125 is the default used in gdalwarp.cpp on line 283
231 */
232 if (max_err < 0.) max_err = 0.125;
233 RASTER_DEBUGF(4, "max_err = %f", max_err);
234
235 /* handle srs */
236 if (src_srs != NULL) {
237 /* reprojection taking place */
238 if (dst_srs != NULL && strcmp(src_srs, dst_srs) != 0) {
239 RASTER_DEBUG(4, "Warp operation does include a reprojection");
240 arg->src.srs = rt_util_gdal_convert_sr(src_srs, 0);
241 arg->dst.srs = rt_util_gdal_convert_sr(dst_srs, 0);
242
243 if (arg->src.srs == NULL || arg->dst.srs == NULL) {
244 rterror("rt_raster_gdal_warp: Could not convert srs values to GDAL accepted format");
245 _rti_warp_arg_destroy(arg);
246 return NULL;
247 }
248 }
249 /* no reprojection, a stub just for clarity */
250 else {
251 RASTER_DEBUG(4, "Warp operation does NOT include reprojection");
252 }
253 }
254 else if (dst_srs != NULL) {
255 /* dst_srs provided but not src_srs */
256 rterror("rt_raster_gdal_warp: SRS required for input raster if SRS provided for warped raster");
257 _rti_warp_arg_destroy(arg);
258 return NULL;
259 }
260
261 /* load raster into a GDAL MEM dataset */
262 arg->src.ds = rt_raster_to_gdal_mem(raster, arg->src.srs, NULL, NULL, 0, &(arg->src.drv), &(arg->src.destroy_drv));
263 if (NULL == arg->src.ds) {
264 rterror("rt_raster_gdal_warp: Could not convert raster to GDAL MEM format");
265 _rti_warp_arg_destroy(arg);
266 return NULL;
267 }
268 RASTER_DEBUG(3, "raster loaded into GDAL MEM dataset");
269
270 /* special case when src_srs and dst_srs is NULL and raster's geotransform matrix is default */
271 if (
272 src_srs == NULL && dst_srs == NULL &&
273 rt_raster_get_srid(raster) == SRID_UNKNOWN
274 ) {
275 double gt[6];
276
277 #if POSTGIS_DEBUG_LEVEL > 3
278 GDALGetGeoTransform(arg->src.ds, gt);
279 RASTER_DEBUGF(3, "GDAL MEM geotransform: %f, %f, %f, %f, %f, %f",
280 gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
281 #endif
282
283 /* default geotransform */
284 rt_raster_get_geotransform_matrix(raster, gt);
285 RASTER_DEBUGF(3, "raster geotransform: %f, %f, %f, %f, %f, %f",
286 gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
287
288 /* substitute spatial info (lack of) with a real one EPSG:32731 (WGS84/UTM zone 31s) */
289 if (FLT_EQ(gt[0], 0.0) && FLT_EQ(gt[3], 0.0) && FLT_EQ(gt[1], 1.0) && FLT_EQ(gt[5], -1.0) &&
290 FLT_EQ(gt[2], 0.0) && FLT_EQ(gt[4], 0.0))
291 {
292 double ngt[6] = {166021.4431, 0.1, 0, 10000000.0000, 0, -0.1};
293
294 rtwarn("Raster has default geotransform. Adjusting metadata for use of GDAL Warp API");
295
296 subspatial = 1;
297
298 GDALSetGeoTransform(arg->src.ds, ngt);
299 GDALFlushCache(arg->src.ds);
300
301 /* EPSG:32731 */
302 arg->src.srs = rt_util_gdal_convert_sr("EPSG:32731", 0);
303 arg->dst.srs = rt_util_gdal_convert_sr("EPSG:32731", 0);
304
305 #if POSTGIS_DEBUG_LEVEL > 3
306 GDALGetGeoTransform(arg->src.ds, gt);
307 RASTER_DEBUGF(3, "GDAL MEM geotransform: %f, %f, %f, %f, %f, %f",
308 gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
309 #endif
310 }
311 }
312
313 /* set transform options */
314 if (arg->src.srs != NULL || arg->dst.srs != NULL) {
315 arg->transform.option.len = 2;
316 arg->transform.option.item = rtalloc(sizeof(char *) * (arg->transform.option.len + 1));
317 if (NULL == arg->transform.option.item) {
318 rterror("rt_raster_gdal_warp: Could not allocation memory for transform options");
319 _rti_warp_arg_destroy(arg);
320 return NULL;
321 }
322 memset(arg->transform.option.item, 0, sizeof(char *) * (arg->transform.option.len + 1));
323
324 for (i = 0; i < arg->transform.option.len; i++) {
325 const char *srs = i ? arg->dst.srs : arg->src.srs;
326 const char *lbl = i ? "DST_SRS=" : "SRC_SRS=";
327 size_t sz = sizeof(char) * (strlen(lbl) + 1);
328 if ( srs ) sz += strlen(srs);
329 arg->transform.option.item[i] = (char *) rtalloc(sz);
330 if (NULL == arg->transform.option.item[i]) {
331 rterror("rt_raster_gdal_warp: Could not allocation memory for transform options");
332 _rti_warp_arg_destroy(arg);
333 return NULL;
334 }
335 sprintf(arg->transform.option.item[i], "%s%s", lbl, srs ? srs : "");
336 RASTER_DEBUGF(4, "arg->transform.option.item[%d] = %s", i, arg->transform.option.item[i]);
337 }
338 }
339 else
340 arg->transform.option.len = 0;
341
342 /* transformation object for building dst dataset */
343 arg->transform.arg.transform = GDALCreateGenImgProjTransformer2(arg->src.ds, NULL, arg->transform.option.item);
344 if (NULL == arg->transform.arg.transform) {
345 rterror("rt_raster_gdal_warp: Could not create GDAL transformation object for output dataset creation");
346 _rti_warp_arg_destroy(arg);
347 return NULL;
348 }
349
350 /* get approximate output georeferenced bounds and resolution */
351 cplerr = GDALSuggestedWarpOutput2(
352 arg->src.ds, GDALGenImgProjTransform,
353 arg->transform.arg.transform, _gt, &(_dim[0]), &(_dim[1]), dst_extent, 0);
354 if (cplerr != CE_None) {
355 rterror("rt_raster_gdal_warp: Could not get GDAL suggested warp output for output dataset creation");
356 _rti_warp_arg_destroy(arg);
357 return NULL;
358 }
359 GDALDestroyGenImgProjTransformer(arg->transform.arg.transform);
360 arg->transform.arg.transform = NULL;
361
362 /*
363 don't use suggested dimensions as use of suggested scales
364 on suggested extent will result in suggested dimensions
365 */
366 _dim[0] = 0;
367 _dim[1] = 0;
368
369 RASTER_DEBUGF(3, "Suggested geotransform: %f, %f, %f, %f, %f, %f",
370 _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
371
372 /* store extent in easier-to-use object */
373 extent.MinX = dst_extent[0];
374 extent.MinY = dst_extent[1];
375 extent.MaxX = dst_extent[2];
376 extent.MaxY = dst_extent[3];
377
378 extent.UpperLeftX = dst_extent[0];
379 extent.UpperLeftY = dst_extent[3];
380
381 RASTER_DEBUGF(3, "Suggested extent: %f, %f, %f, %f",
382 extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
383
384 /* scale and width/height are mutually exclusive */
385 if (
386 ((NULL != scale_x) || (NULL != scale_y)) &&
387 ((NULL != width) || (NULL != height))
388 ) {
389 rterror("rt_raster_gdal_warp: Scale X/Y and width/height are mutually exclusive. Only provide one");
390 _rti_warp_arg_destroy(arg);
391 return NULL;
392 }
393
394 /* user-defined width */
395 if (NULL != width) {
396 _dim[0] = abs(*width);
397 _scale[0] = fabs((extent.MaxX - extent.MinX) / ((double) _dim[0]));
398 }
399 /* user-defined height */
400 if (NULL != height) {
401 _dim[1] = abs(*height);
402 _scale[1] = fabs((extent.MaxY - extent.MinY) / ((double) _dim[1]));
403 }
404
405 /* user-defined scale */
406 if (
407 ((NULL != scale_x) && (FLT_NEQ(*scale_x, 0.0))) &&
408 ((NULL != scale_y) && (FLT_NEQ(*scale_y, 0.0)))
409 ) {
410 _scale[0] = fabs(*scale_x);
411 _scale[1] = fabs(*scale_y);
412
413 /* special override since we changed the original GT scales */
414 if (subspatial) {
415 /*
416 _scale[0] *= 10;
417 _scale[1] *= 10;
418 */
419 _scale[0] /= 10;
420 _scale[1] /= 10;
421 }
422 }
423 else if (
424 ((NULL != scale_x) && (NULL == scale_y)) ||
425 ((NULL == scale_x) && (NULL != scale_y))
426 ) {
427 rterror("rt_raster_gdal_warp: Both X and Y scale values must be provided for scale");
428 _rti_warp_arg_destroy(arg);
429 return NULL;
430 }
431
432 /* scale not defined, use suggested */
433 if (FLT_EQ(_scale[0], 0.0) && FLT_EQ(_scale[1], 0.0))
434 {
435 _scale[0] = fabs(_gt[1]);
436 _scale[1] = fabs(_gt[5]);
437 }
438
439 RASTER_DEBUGF(4, "Using scale: %f x %f", _scale[0], -1 * _scale[1]);
440
441 /* user-defined skew */
442 if (NULL != skew_x) {
443 _skew[0] = *skew_x;
444
445 /*
446 negative scale-x affects skew
447 for now, force skew to be in left-right, top-down orientation
448 */
449 if (
450 NULL != scale_x &&
451 *scale_x < 0.
452 ) {
453 _skew[0] *= -1;
454 }
455 }
456 if (NULL != skew_y) {
457 _skew[1] = *skew_y;
458
459 /*
460 positive scale-y affects skew
461 for now, force skew to be in left-right, top-down orientation
462 */
463 if (
464 NULL != scale_y &&
465 *scale_y > 0.
466 ) {
467 _skew[1] *= -1;
468 }
469 }
470
471 RASTER_DEBUGF(4, "Using skew: %f x %f", _skew[0], _skew[1]);
472
473 /* reprocess extent if skewed */
474 if (FLT_NEQ(_skew[0], 0.0) || FLT_NEQ(_skew[1], 0.0))
475 {
476 rt_raster skewedrast;
477
478 RASTER_DEBUG(3, "Computing skewed extent's envelope");
479
480 skewedrast = rt_raster_compute_skewed_raster(
481 extent,
482 _skew,
483 _scale,
484 0.01
485 );
486 if (skewedrast == NULL) {
487 rterror("rt_raster_gdal_warp: Could not compute skewed raster");
488 _rti_warp_arg_destroy(arg);
489 return NULL;
490 }
491
492 if (_dim[0] == 0)
493 _dim[0] = skewedrast->width;
494 if (_dim[1] == 0)
495 _dim[1] = skewedrast->height;
496
497 extent.UpperLeftX = skewedrast->ipX;
498 extent.UpperLeftY = skewedrast->ipY;
499
500 rt_raster_destroy(skewedrast);
501 }
502
503 /* dimensions not defined, compute */
504 if (!_dim[0])
505 _dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
506 if (!_dim[1])
507 _dim[1] = (int) fmax((fabs(extent.MaxY - extent.MinY) + (_scale[1] / 2.)) / _scale[1], 1);
508
509 /* temporary raster */
510 rast = rt_raster_new(_dim[0], _dim[1]);
511 if (rast == NULL) {
512 rterror("rt_raster_gdal_warp: Out of memory allocating temporary raster");
513 _rti_warp_arg_destroy(arg);
514 return NULL;
515 }
516
517 /* set raster's spatial attributes */
518 rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
519 rt_raster_set_scale(rast, _scale[0], -1 * _scale[1]);
520 rt_raster_set_skews(rast, _skew[0], _skew[1]);
521
522 rt_raster_get_geotransform_matrix(rast, _gt);
523 RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
524 _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
525 RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
526 _dim[0], _dim[1]);
527
528 /* user-defined upper-left corner */
529 if (
530 NULL != ul_xw &&
531 NULL != ul_yw
532 ) {
533 ul_user = 1;
534
535 RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
536
537 /* set upper-left corner */
538 rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
539 extent.UpperLeftX = *ul_xw;
540 extent.UpperLeftY = *ul_yw;
541 }
542 else if (
543 ((NULL != ul_xw) && (NULL == ul_yw)) ||
544 ((NULL == ul_xw) && (NULL != ul_yw))
545 ) {
546 rterror("rt_raster_gdal_warp: Both X and Y upper-left corner values must be provided");
547 rt_raster_destroy(rast);
548 _rti_warp_arg_destroy(arg);
549 return NULL;
550 }
551
552 /* alignment only considered if upper-left corner not provided */
553 if (
554 !ul_user && (
555 (NULL != grid_xw) || (NULL != grid_yw)
556 )
557 ) {
558
559 if (
560 ((NULL != grid_xw) && (NULL == grid_yw)) ||
561 ((NULL == grid_xw) && (NULL != grid_yw))
562 ) {
563 rterror("rt_raster_gdal_warp: Both X and Y alignment values must be provided");
564 rt_raster_destroy(rast);
565 _rti_warp_arg_destroy(arg);
566 return NULL;
567 }
568
569 RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw);
570
571 do {
572 double _r[2] = {0};
573 double _w[2] = {0};
574
575 /* raster is already aligned */
576 if (FLT_EQ(*grid_xw, extent.UpperLeftX) && FLT_EQ(*grid_yw, extent.UpperLeftY)) {
577 RASTER_DEBUG(3, "Skipping raster alignment as it is already aligned to grid");
578 break;
579 }
580
581 extent.UpperLeftX = rast->ipX;
582 extent.UpperLeftY = rast->ipY;
583 rt_raster_set_offsets(rast, *grid_xw, *grid_yw);
584
585 /* process upper-left corner */
586 if (rt_raster_geopoint_to_cell(
587 rast,
588 extent.UpperLeftX, extent.UpperLeftY,
589 &(_r[0]), &(_r[1]),
590 NULL
591 ) != ES_NONE) {
592 rterror("rt_raster_gdal_warp: Could not compute raster pixel for spatial coordinates");
593 rt_raster_destroy(rast);
594 _rti_warp_arg_destroy(arg);
595 return NULL;
596 }
597
598 if (rt_raster_cell_to_geopoint(
599 rast,
600 _r[0], _r[1],
601 &(_w[0]), &(_w[1]),
602 NULL
603 ) != ES_NONE) {
604 rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
605
606 rt_raster_destroy(rast);
607 _rti_warp_arg_destroy(arg);
608 return NULL;
609 }
610
611 /* shift occurred */
612 if (FLT_NEQ(_w[0], extent.UpperLeftX)) {
613 if (NULL == width)
614 rast->width++;
615 else if (NULL == scale_x) {
616 double _c[2] = {0};
617
618 rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
619
620 /* get upper-right corner */
621 if (rt_raster_cell_to_geopoint(
622 rast,
623 rast->width, 0,
624 &(_c[0]), &(_c[1]),
625 NULL
626 ) != ES_NONE) {
627 rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
628 rt_raster_destroy(rast);
629 _rti_warp_arg_destroy(arg);
630 return NULL;
631 }
632
633 rast->scaleX = fabs((_c[0] - _w[0]) / ((double) rast->width));
634 }
635 }
636 if (FLT_NEQ(_w[1], extent.UpperLeftY)) {
637 if (NULL == height)
638 rast->height++;
639 else if (NULL == scale_y) {
640 double _c[2] = {0};
641
642 rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
643
644 /* get upper-right corner */
645 if (rt_raster_cell_to_geopoint(
646 rast,
647 0, rast->height,
648 &(_c[0]), &(_c[1]),
649 NULL
650 ) != ES_NONE) {
651 rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
652
653 rt_raster_destroy(rast);
654 _rti_warp_arg_destroy(arg);
655 return NULL;
656 }
657
658 rast->scaleY = -1 * fabs((_c[1] - _w[1]) / ((double) rast->height));
659 }
660 }
661
662 rt_raster_set_offsets(rast, _w[0], _w[1]);
663 RASTER_DEBUGF(4, "aligned offsets: %f x %f", _w[0], _w[1]);
664 }
665 while (0);
666 }
667
668 /*
669 after this point, rt_envelope extent is no longer used
670 */
671
672 /* get key attributes from rast */
673 _dim[0] = rast->width;
674 _dim[1] = rast->height;
675 rt_raster_get_geotransform_matrix(rast, _gt);
676
677 /* scale-x is negative or scale-y is positive */
678 if ((
679 (NULL != scale_x) && (*scale_x < 0.)
680 ) || (
681 (NULL != scale_y) && (*scale_y > 0)
682 )) {
683 double _w[2] = {0};
684
685 /* negative scale-x */
686 if (
687 (NULL != scale_x) &&
688 (*scale_x < 0.)
689 ) {
690 if (rt_raster_cell_to_geopoint(
691 rast,
692 rast->width, 0,
693 &(_w[0]), &(_w[1]),
694 NULL
695 ) != ES_NONE) {
696 rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
697 rt_raster_destroy(rast);
698 _rti_warp_arg_destroy(arg);
699 return NULL;
700 }
701
702 _gt[0] = _w[0];
703 _gt[1] = *scale_x;
704
705 /* check for skew */
706 if (NULL != skew_x && FLT_NEQ(*skew_x, 0.0))
707 _gt[2] = *skew_x;
708 }
709 /* positive scale-y */
710 if (
711 (NULL != scale_y) &&
712 (*scale_y > 0)
713 ) {
714 if (rt_raster_cell_to_geopoint(
715 rast,
716 0, rast->height,
717 &(_w[0]), &(_w[1]),
718 NULL
719 ) != ES_NONE) {
720 rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
721 rt_raster_destroy(rast);
722 _rti_warp_arg_destroy(arg);
723 return NULL;
724 }
725
726 _gt[3] = _w[1];
727 _gt[5] = *scale_y;
728
729 /* check for skew */
730 if (NULL != skew_y && FLT_NEQ(*skew_y, 0.0))
731 _gt[4] = *skew_y;
732 }
733 }
734
735 rt_raster_destroy(rast);
736 rast = NULL;
737
738 RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
739 _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
740 RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
741 _dim[0], _dim[1]);
742
743 if ( _dim[0] == 0 || _dim[1] == 0 ) {
744 rterror("rt_raster_gdal_warp: The width (%d) or height (%d) of the warped raster is zero", _dim[0], _dim[1]);
745 _rti_warp_arg_destroy(arg);
746 return NULL;
747 }
748
749 /* load VRT driver */
750 if (!rt_util_gdal_driver_registered("VRT")) {
751 RASTER_DEBUG(3, "Registering VRT driver");
752 GDALRegister_VRT();
753 arg->dst.destroy_drv = 1;
754 }
755 arg->dst.drv = GDALGetDriverByName("VRT");
756 if (NULL == arg->dst.drv) {
757 rterror("rt_raster_gdal_warp: Could not load the output GDAL VRT driver");
758 _rti_warp_arg_destroy(arg);
759 return NULL;
760 }
761
762 /* create dst dataset */
763 arg->dst.ds = GDALCreate(arg->dst.drv, "", _dim[0], _dim[1], 0, GDT_Byte, dst_options);
764 if (NULL == arg->dst.ds) {
765 rterror("rt_raster_gdal_warp: Could not create GDAL VRT dataset");
766 _rti_warp_arg_destroy(arg);
767 return NULL;
768 }
769
770 /* set dst srs */
771 if (arg->dst.srs != NULL) {
772 cplerr = GDALSetProjection(arg->dst.ds, arg->dst.srs);
773 if (cplerr != CE_None) {
774 rterror("rt_raster_gdal_warp: Could not set projection");
775 _rti_warp_arg_destroy(arg);
776 return NULL;
777 }
778 RASTER_DEBUGF(3, "Applied SRS: %s", GDALGetProjectionRef(arg->dst.ds));
779 }
780
781 /* set dst geotransform */
782 cplerr = GDALSetGeoTransform(arg->dst.ds, _gt);
783 if (cplerr != CE_None) {
784 rterror("rt_raster_gdal_warp: Could not set geotransform");
785 _rti_warp_arg_destroy(arg);
786 return NULL;
787 }
788
789 /* add bands to dst dataset */
790 numBands = rt_raster_get_num_bands(raster);
791 for (i = 0; i < numBands; i++) {
792 rtband = rt_raster_get_band(raster, i);
793 if (NULL == rtband) {
794 rterror("rt_raster_gdal_warp: Could not get band %d for adding to VRT dataset", i);
795 _rti_warp_arg_destroy(arg);
796 return NULL;
797 }
798
799 pt = rt_band_get_pixtype(rtband);
800 gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
801 if (gdal_pt == GDT_Unknown)
802 rtwarn("rt_raster_gdal_warp: Unknown pixel type for band %d", i);
803
804 cplerr = GDALAddBand(arg->dst.ds, gdal_pt, NULL);
805 if (cplerr != CE_None) {
806 rterror("rt_raster_gdal_warp: Could not add band to VRT dataset");
807 _rti_warp_arg_destroy(arg);
808 return NULL;
809 }
810
811 /* get band to write data to */
812 band = NULL;
813 band = GDALGetRasterBand(arg->dst.ds, i + 1);
814 if (NULL == band) {
815 rterror("rt_raster_gdal_warp: Could not get GDAL band for additional processing");
816 _rti_warp_arg_destroy(arg);
817 return NULL;
818 }
819
820 /* set nodata */
821 if (rt_band_get_hasnodata_flag(rtband) != FALSE) {
822 hasnodata = 1;
823 rt_band_get_nodata(rtband, &nodata);
824 if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
825 rtwarn("rt_raster_gdal_warp: Could not set nodata value for band %d", i);
826 RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
827 }
828 }
829
830 /* create transformation object */
831 arg->transform.arg.transform = arg->transform.arg.imgproj = GDALCreateGenImgProjTransformer2(
832 arg->src.ds, arg->dst.ds,
833 arg->transform.option.item
834 );
835 if (NULL == arg->transform.arg.transform) {
836 rterror("rt_raster_gdal_warp: Could not create GDAL transformation object");
837 _rti_warp_arg_destroy(arg);
838 return NULL;
839 }
840 arg->transform.func = GDALGenImgProjTransform;
841
842 /* use approximate transformation object */
843 if (max_err > 0.0) {
844 arg->transform.arg.transform = arg->transform.arg.approx = GDALCreateApproxTransformer(
845 GDALGenImgProjTransform,
846 arg->transform.arg.imgproj, max_err
847 );
848 if (NULL == arg->transform.arg.transform) {
849 rterror("rt_raster_gdal_warp: Could not create GDAL approximate transformation object");
850 _rti_warp_arg_destroy(arg);
851 return NULL;
852 }
853
854 arg->transform.func = GDALApproxTransform;
855 }
856
857 /* warp options */
858 arg->wopts = GDALCreateWarpOptions();
859 if (NULL == arg->wopts) {
860 rterror("rt_raster_gdal_warp: Could not create GDAL warp options object");
861 _rti_warp_arg_destroy(arg);
862 return NULL;
863 }
864
865 /* set options */
866 arg->wopts->eResampleAlg = resample_alg;
867 arg->wopts->hSrcDS = arg->src.ds;
868 arg->wopts->hDstDS = arg->dst.ds;
869 arg->wopts->pfnTransformer = arg->transform.func;
870 arg->wopts->pTransformerArg = arg->transform.arg.transform;
871 arg->wopts->papszWarpOptions = (char **) CPLMalloc(sizeof(char *) * 2);
872 arg->wopts->papszWarpOptions[0] = (char *) CPLMalloc(sizeof(char) * (strlen("INIT_DEST=NO_DATA") + 1));
873 strcpy(arg->wopts->papszWarpOptions[0], "INIT_DEST=NO_DATA");
874 arg->wopts->papszWarpOptions[1] = NULL;
875
876 /* set band mapping */
877 arg->wopts->nBandCount = numBands;
878 arg->wopts->panSrcBands = (int *) CPLMalloc(sizeof(int) * arg->wopts->nBandCount);
879 arg->wopts->panDstBands = (int *) CPLMalloc(sizeof(int) * arg->wopts->nBandCount);
880 for (i = 0; i < arg->wopts->nBandCount; i++)
881 arg->wopts->panDstBands[i] = arg->wopts->panSrcBands[i] = i + 1;
882
883 /* set nodata mapping */
884 if (hasnodata) {
885 RASTER_DEBUG(3, "Setting nodata mapping");
886 arg->wopts->padfSrcNoDataReal = (double *) CPLMalloc(numBands * sizeof(double));
887 arg->wopts->padfDstNoDataReal = (double *) CPLMalloc(numBands * sizeof(double));
888 arg->wopts->padfSrcNoDataImag = (double *) CPLMalloc(numBands * sizeof(double));
889 arg->wopts->padfDstNoDataImag = (double *) CPLMalloc(numBands * sizeof(double));
890 if (
891 NULL == arg->wopts->padfSrcNoDataReal ||
892 NULL == arg->wopts->padfDstNoDataReal ||
893 NULL == arg->wopts->padfSrcNoDataImag ||
894 NULL == arg->wopts->padfDstNoDataImag
895 ) {
896 rterror("rt_raster_gdal_warp: Out of memory allocating nodata mapping");
897 _rti_warp_arg_destroy(arg);
898 return NULL;
899 }
900 for (i = 0; i < numBands; i++) {
901 band = rt_raster_get_band(raster, i);
902 if (!band) {
903 rterror("rt_raster_gdal_warp: Could not process bands for nodata values");
904 _rti_warp_arg_destroy(arg);
905 return NULL;
906 }
907
908 if (!rt_band_get_hasnodata_flag(band)) {
909 /*
910 based on line 1004 of gdalwarp.cpp
911 the problem is that there is a chance that this number is a legitimate value
912 */
913 arg->wopts->padfSrcNoDataReal[i] = -123456.789;
914 }
915 else {
916 rt_band_get_nodata(band, &(arg->wopts->padfSrcNoDataReal[i]));
917 }
918
919 arg->wopts->padfDstNoDataReal[i] = arg->wopts->padfSrcNoDataReal[i];
920 arg->wopts->padfDstNoDataImag[i] = arg->wopts->padfSrcNoDataImag[i] = 0.0;
921 RASTER_DEBUGF(4, "Mapped nodata value for band %d: %f (%f) => %f (%f)",
922 i,
923 arg->wopts->padfSrcNoDataReal[i], arg->wopts->padfSrcNoDataImag[i],
924 arg->wopts->padfDstNoDataReal[i], arg->wopts->padfDstNoDataImag[i]
925 );
926 }
927 }
928
929 /* warp raster */
930 RASTER_DEBUG(3, "Warping raster");
931 cplerr = GDALInitializeWarpedVRT(arg->dst.ds, arg->wopts);
932 if (cplerr != CE_None) {
933 rterror("rt_raster_gdal_warp: Could not warp raster");
934 _rti_warp_arg_destroy(arg);
935 return NULL;
936 }
937 /*
938 GDALSetDescription(arg->dst.ds, "/tmp/warped.vrt");
939 */
940 GDALFlushCache(arg->dst.ds);
941 RASTER_DEBUG(3, "Raster warped");
942
943 /* convert gdal dataset to raster */
944 RASTER_DEBUG(3, "Converting GDAL dataset to raster");
945 rast = rt_raster_from_gdal_dataset(arg->dst.ds);
946
947 _rti_warp_arg_destroy(arg);
948
949 if (NULL == rast) {
950 rterror("rt_raster_gdal_warp: Could not warp raster");
951 return NULL;
952 }
953
954 /* substitute spatial, reset back to default */
955 if (subspatial) {
956 double gt[6] = {0, 1, 0, 0, 0, -1};
957 /* See http://trac.osgeo.org/postgis/ticket/2911 */
958 /* We should proably also tweak rotation here */
959 /* NOTE: the times 10 is because it was divided by 10 in a section above,
960 * I'm not sure the above division was needed */
961 gt[1] = _scale[0] * 10;
962 gt[5] = -1 * _scale[1] * 10;
963
964 rt_raster_set_geotransform_matrix(rast, gt);
965 rt_raster_set_srid(rast, SRID_UNKNOWN);
966 }
967
968 RASTER_DEBUG(3, "done");
969
970 return rast;
971 }
972
973