PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ 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.

Parameters
raster: raster to transform
src_srs: the raster's coordinate system in OGC WKT
dst_srs: the warped raster's coordinate system in OGC WKT
scale_x: the x size of pixels of the warped raster's pixels in units of dst_srs
scale_y: the y size of pixels of the warped raster's pixels in units of dst_srs
width: the number of columns of the warped raster. note that width/height CANNOT be used with scale_x/scale_y
height: the number of rows of the warped raster. note that width/height CANNOT be used with scale_x/scale_y
ul_xw: the X value of upper-left corner of the warped raster in units of dst_srs
ul_yw: the Y value of upper-left corner of the warped raster in units of dst_srs
grid_xw: the X value of point on a grid to align warped raster to in units of dst_srs
grid_yw: the Y value of point on a grid to align warped raster to in units of dst_srs
skew_x: the X skew of the warped raster in units of dst_srs
skew_y: the Y skew of the warped raster in units of dst_srs
resample_alg: the resampling algorithm
max_err: maximum error measured in input pixels permitted (0.0 for exact calculations)
Returns
the warped raster or NULL

Definition at line 178 of file rt_warp.c.

References _rti_warp_arg_destroy(), _rti_warp_arg_init(), _rti_warp_arg_t::arg, ovdump::band, _rti_warp_arg_t::destroy_drv, _rti_warp_arg_t::drv, _rti_warp_arg_t::ds, _rti_warp_arg_t::dst, ES_NONE, FALSE, FLT_EQ, FLT_NEQ, window::gt, rt_raster_t::height, if(), rt_raster_t::ipX, rt_raster_t::ipY, rt_envelope::MaxX, rt_envelope::MaxY, rt_envelope::MinX, rt_envelope::MinY, PT_END, rtpixdump::rast, RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_hasnodata_flag(), rt_band_get_nodata(), rt_band_get_pixtype(), rt_raster_cell_to_geopoint(), rt_raster_compute_skewed_raster(), rt_raster_destroy(), rt_raster_from_gdal_dataset(), rt_raster_geopoint_to_cell(), rt_raster_get_band(), rt_raster_get_geotransform_matrix(), rt_raster_get_num_bands(), rt_raster_get_srid(), rt_raster_new(), rt_raster_set_geotransform_matrix(), rt_raster_set_offsets(), rt_raster_set_scale(), rt_raster_set_skews(), rt_raster_set_srid(), rt_raster_to_gdal_mem(), rt_util_gdal_convert_sr(), rt_util_gdal_driver_registered(), rt_util_pixtype_to_gdal_datatype(), rtalloc(), rterror(), rtwarn(), rt_raster_t::scaleX, rt_raster_t::scaleY, _rti_warp_arg_t::src, SRID_UNKNOWN, _rti_warp_arg_t::srs, _rti_warp_arg_t::transform, rt_envelope::UpperLeftX, rt_envelope::UpperLeftY, rt_raster_t::width, and _rti_warp_arg_t::wopts.

Referenced by RASTER_GDALWarp(), and test_gdal_warp().

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 (
290  FLT_EQ(gt[0], 0) && FLT_EQ(gt[3], 0) &&
291  FLT_EQ(gt[1], 1) && FLT_EQ(gt[5], -1) &&
292  FLT_EQ(gt[2], 0) && FLT_EQ(gt[4], 0)
293  ) {
294  double ngt[6] = {166021.4431, 0.1, 0, 10000000.0000, 0, -0.1};
295 
296  rtwarn("Raster has default geotransform. Adjusting metadata for use of GDAL Warp API");
297 
298  subspatial = 1;
299 
300  GDALSetGeoTransform(arg->src.ds, ngt);
301  GDALFlushCache(arg->src.ds);
302 
303  /* EPSG:32731 */
304  arg->src.srs = rt_util_gdal_convert_sr("EPSG:32731", 0);
305  arg->dst.srs = rt_util_gdal_convert_sr("EPSG:32731", 0);
306 
307 #if POSTGIS_DEBUG_LEVEL > 3
308  GDALGetGeoTransform(arg->src.ds, gt);
309  RASTER_DEBUGF(3, "GDAL MEM geotransform: %f, %f, %f, %f, %f, %f",
310  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
311 #endif
312  }
313  }
314 
315  /* set transform options */
316  if (arg->src.srs != NULL || arg->dst.srs != NULL) {
317  arg->transform.option.len = 2;
318  arg->transform.option.item = rtalloc(sizeof(char *) * (arg->transform.option.len + 1));
319  if (NULL == arg->transform.option.item) {
320  rterror("rt_raster_gdal_warp: Could not allocation memory for transform options");
321  _rti_warp_arg_destroy(arg);
322  return NULL;
323  }
324  memset(arg->transform.option.item, 0, sizeof(char *) * (arg->transform.option.len + 1));
325 
326  for (i = 0; i < arg->transform.option.len; i++) {
327  const char *srs = i ? arg->dst.srs : arg->src.srs;
328  const char *lbl = i ? "DST_SRS=" : "SRC_SRS=";
329  size_t sz = sizeof(char) * (strlen(lbl) + 1);
330  if ( srs ) sz += strlen(srs);
331  arg->transform.option.item[i] = (char *) rtalloc(sz);
332  if (NULL == arg->transform.option.item[i]) {
333  rterror("rt_raster_gdal_warp: Could not allocation memory for transform options");
334  _rti_warp_arg_destroy(arg);
335  return NULL;
336  }
337  sprintf(arg->transform.option.item[i], "%s%s", lbl, srs ? srs : "");
338  RASTER_DEBUGF(4, "arg->transform.option.item[%d] = %s", i, arg->transform.option.item[i]);
339  }
340  }
341  else
342  arg->transform.option.len = 0;
343 
344  /* transformation object for building dst dataset */
345  arg->transform.arg.transform = GDALCreateGenImgProjTransformer2(arg->src.ds, NULL, arg->transform.option.item);
346  if (NULL == arg->transform.arg.transform) {
347  rterror("rt_raster_gdal_warp: Could not create GDAL transformation object for output dataset creation");
348  _rti_warp_arg_destroy(arg);
349  return NULL;
350  }
351 
352  /* get approximate output georeferenced bounds and resolution */
353  cplerr = GDALSuggestedWarpOutput2(
354  arg->src.ds, GDALGenImgProjTransform,
355  arg->transform.arg.transform, _gt, &(_dim[0]), &(_dim[1]), dst_extent, 0);
356  if (cplerr != CE_None) {
357  rterror("rt_raster_gdal_warp: Could not get GDAL suggested warp output for output dataset creation");
358  _rti_warp_arg_destroy(arg);
359  return NULL;
360  }
361  GDALDestroyGenImgProjTransformer(arg->transform.arg.transform);
362  arg->transform.arg.transform = NULL;
363 
364  /*
365  don't use suggested dimensions as use of suggested scales
366  on suggested extent will result in suggested dimensions
367  */
368  _dim[0] = 0;
369  _dim[1] = 0;
370 
371  RASTER_DEBUGF(3, "Suggested geotransform: %f, %f, %f, %f, %f, %f",
372  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
373 
374  /* store extent in easier-to-use object */
375  extent.MinX = dst_extent[0];
376  extent.MinY = dst_extent[1];
377  extent.MaxX = dst_extent[2];
378  extent.MaxY = dst_extent[3];
379 
380  extent.UpperLeftX = dst_extent[0];
381  extent.UpperLeftY = dst_extent[3];
382 
383  RASTER_DEBUGF(3, "Suggested extent: %f, %f, %f, %f",
384  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
385 
386  /* scale and width/height are mutually exclusive */
387  if (
388  ((NULL != scale_x) || (NULL != scale_y)) &&
389  ((NULL != width) || (NULL != height))
390  ) {
391  rterror("rt_raster_gdal_warp: Scale X/Y and width/height are mutually exclusive. Only provide one");
392  _rti_warp_arg_destroy(arg);
393  return NULL;
394  }
395 
396  /* user-defined width */
397  if (NULL != width) {
398  _dim[0] = abs(*width);
399  _scale[0] = fabs((extent.MaxX - extent.MinX) / ((double) _dim[0]));
400  }
401  /* user-defined height */
402  if (NULL != height) {
403  _dim[1] = abs(*height);
404  _scale[1] = fabs((extent.MaxY - extent.MinY) / ((double) _dim[1]));
405  }
406 
407  /* user-defined scale */
408  if (
409  ((NULL != scale_x) && (FLT_NEQ(*scale_x, 0.0))) &&
410  ((NULL != scale_y) && (FLT_NEQ(*scale_y, 0.0)))
411  ) {
412  _scale[0] = fabs(*scale_x);
413  _scale[1] = fabs(*scale_y);
414 
415  /* special override since we changed the original GT scales */
416  if (subspatial) {
417  /*
418  _scale[0] *= 10;
419  _scale[1] *= 10;
420  */
421  _scale[0] /= 10;
422  _scale[1] /= 10;
423  }
424  }
425  else if (
426  ((NULL != scale_x) && (NULL == scale_y)) ||
427  ((NULL == scale_x) && (NULL != scale_y))
428  ) {
429  rterror("rt_raster_gdal_warp: Both X and Y scale values must be provided for scale");
430  _rti_warp_arg_destroy(arg);
431  return NULL;
432  }
433 
434  /* scale not defined, use suggested */
435  if (FLT_EQ(_scale[0], 0) && FLT_EQ(_scale[1], 0)) {
436  _scale[0] = fabs(_gt[1]);
437  _scale[1] = fabs(_gt[5]);
438  }
439 
440  RASTER_DEBUGF(4, "Using scale: %f x %f", _scale[0], -1 * _scale[1]);
441 
442  /* user-defined skew */
443  if (NULL != skew_x) {
444  _skew[0] = *skew_x;
445 
446  /*
447  negative scale-x affects skew
448  for now, force skew to be in left-right, top-down orientation
449  */
450  if (
451  NULL != scale_x &&
452  *scale_x < 0.
453  ) {
454  _skew[0] *= -1;
455  }
456  }
457  if (NULL != skew_y) {
458  _skew[1] = *skew_y;
459 
460  /*
461  positive scale-y affects skew
462  for now, force skew to be in left-right, top-down orientation
463  */
464  if (
465  NULL != scale_y &&
466  *scale_y > 0.
467  ) {
468  _skew[1] *= -1;
469  }
470  }
471 
472  RASTER_DEBUGF(4, "Using skew: %f x %f", _skew[0], _skew[1]);
473 
474  /* reprocess extent if skewed */
475  if (
476  FLT_NEQ(_skew[0], 0) ||
477  FLT_NEQ(_skew[1], 0)
478  ) {
479  rt_raster skewedrast;
480 
481  RASTER_DEBUG(3, "Computing skewed extent's envelope");
482 
483  skewedrast = rt_raster_compute_skewed_raster(
484  extent,
485  _skew,
486  _scale,
487  0.01
488  );
489  if (skewedrast == NULL) {
490  rterror("rt_raster_gdal_warp: Could not compute skewed raster");
491  _rti_warp_arg_destroy(arg);
492  return NULL;
493  }
494 
495  if (_dim[0] == 0)
496  _dim[0] = skewedrast->width;
497  if (_dim[1] == 0)
498  _dim[1] = skewedrast->height;
499 
500  extent.UpperLeftX = skewedrast->ipX;
501  extent.UpperLeftY = skewedrast->ipY;
502 
503  rt_raster_destroy(skewedrast);
504  }
505 
506  /* dimensions not defined, compute */
507  if (!_dim[0])
508  _dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
509  if (!_dim[1])
510  _dim[1] = (int) fmax((fabs(extent.MaxY - extent.MinY) + (_scale[1] / 2.)) / _scale[1], 1);
511 
512  /* temporary raster */
513  rast = rt_raster_new(_dim[0], _dim[1]);
514  if (rast == NULL) {
515  rterror("rt_raster_gdal_warp: Out of memory allocating temporary raster");
516  _rti_warp_arg_destroy(arg);
517  return NULL;
518  }
519 
520  /* set raster's spatial attributes */
521  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
522  rt_raster_set_scale(rast, _scale[0], -1 * _scale[1]);
523  rt_raster_set_skews(rast, _skew[0], _skew[1]);
524 
525  rt_raster_get_geotransform_matrix(rast, _gt);
526  RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
527  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
528  RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
529  _dim[0], _dim[1]);
530 
531  /* user-defined upper-left corner */
532  if (
533  NULL != ul_xw &&
534  NULL != ul_yw
535  ) {
536  ul_user = 1;
537 
538  RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
539 
540  /* set upper-left corner */
541  rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
542  extent.UpperLeftX = *ul_xw;
543  extent.UpperLeftY = *ul_yw;
544  }
545  else if (
546  ((NULL != ul_xw) && (NULL == ul_yw)) ||
547  ((NULL == ul_xw) && (NULL != ul_yw))
548  ) {
549  rterror("rt_raster_gdal_warp: Both X and Y upper-left corner values must be provided");
550  rt_raster_destroy(rast);
551  _rti_warp_arg_destroy(arg);
552  return NULL;
553  }
554 
555  /* alignment only considered if upper-left corner not provided */
556  if (
557  !ul_user && (
558  (NULL != grid_xw) || (NULL != grid_yw)
559  )
560  ) {
561 
562  if (
563  ((NULL != grid_xw) && (NULL == grid_yw)) ||
564  ((NULL == grid_xw) && (NULL != grid_yw))
565  ) {
566  rterror("rt_raster_gdal_warp: Both X and Y alignment values must be provided");
567  rt_raster_destroy(rast);
568  _rti_warp_arg_destroy(arg);
569  return NULL;
570  }
571 
572  RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw);
573 
574  do {
575  double _r[2] = {0};
576  double _w[2] = {0};
577 
578  /* raster is already aligned */
579  if (FLT_EQ(*grid_xw, extent.UpperLeftX) && FLT_EQ(*grid_yw, extent.UpperLeftY)) {
580  RASTER_DEBUG(3, "Skipping raster alignment as it is already aligned to grid");
581  break;
582  }
583 
584  extent.UpperLeftX = rast->ipX;
585  extent.UpperLeftY = rast->ipY;
586  rt_raster_set_offsets(rast, *grid_xw, *grid_yw);
587 
588  /* process upper-left corner */
589  if (rt_raster_geopoint_to_cell(
590  rast,
591  extent.UpperLeftX, extent.UpperLeftY,
592  &(_r[0]), &(_r[1]),
593  NULL
594  ) != ES_NONE) {
595  rterror("rt_raster_gdal_warp: Could not compute raster pixel for spatial coordinates");
596  rt_raster_destroy(rast);
597  _rti_warp_arg_destroy(arg);
598  return NULL;
599  }
600 
601  if (rt_raster_cell_to_geopoint(
602  rast,
603  _r[0], _r[1],
604  &(_w[0]), &(_w[1]),
605  NULL
606  ) != ES_NONE) {
607  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
608 
609  rt_raster_destroy(rast);
610  _rti_warp_arg_destroy(arg);
611  return NULL;
612  }
613 
614  /* shift occurred */
615  if (FLT_NEQ(_w[0], extent.UpperLeftX)) {
616  if (NULL == width)
617  rast->width++;
618  else if (NULL == scale_x) {
619  double _c[2] = {0};
620 
621  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
622 
623  /* get upper-right corner */
624  if (rt_raster_cell_to_geopoint(
625  rast,
626  rast->width, 0,
627  &(_c[0]), &(_c[1]),
628  NULL
629  ) != ES_NONE) {
630  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
631  rt_raster_destroy(rast);
632  _rti_warp_arg_destroy(arg);
633  return NULL;
634  }
635 
636  rast->scaleX = fabs((_c[0] - _w[0]) / ((double) rast->width));
637  }
638  }
639  if (FLT_NEQ(_w[1], extent.UpperLeftY)) {
640  if (NULL == height)
641  rast->height++;
642  else if (NULL == scale_y) {
643  double _c[2] = {0};
644 
645  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
646 
647  /* get upper-right corner */
648  if (rt_raster_cell_to_geopoint(
649  rast,
650  0, rast->height,
651  &(_c[0]), &(_c[1]),
652  NULL
653  ) != ES_NONE) {
654  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
655 
656  rt_raster_destroy(rast);
657  _rti_warp_arg_destroy(arg);
658  return NULL;
659  }
660 
661  rast->scaleY = -1 * fabs((_c[1] - _w[1]) / ((double) rast->height));
662  }
663  }
664 
665  rt_raster_set_offsets(rast, _w[0], _w[1]);
666  RASTER_DEBUGF(4, "aligned offsets: %f x %f", _w[0], _w[1]);
667  }
668  while (0);
669  }
670 
671  /*
672  after this point, rt_envelope extent is no longer used
673  */
674 
675  /* get key attributes from rast */
676  _dim[0] = rast->width;
677  _dim[1] = rast->height;
678  rt_raster_get_geotransform_matrix(rast, _gt);
679 
680  /* scale-x is negative or scale-y is positive */
681  if ((
682  (NULL != scale_x) && (*scale_x < 0.)
683  ) || (
684  (NULL != scale_y) && (*scale_y > 0)
685  )) {
686  double _w[2] = {0};
687 
688  /* negative scale-x */
689  if (
690  (NULL != scale_x) &&
691  (*scale_x < 0.)
692  ) {
693  if (rt_raster_cell_to_geopoint(
694  rast,
695  rast->width, 0,
696  &(_w[0]), &(_w[1]),
697  NULL
698  ) != ES_NONE) {
699  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
700  rt_raster_destroy(rast);
701  _rti_warp_arg_destroy(arg);
702  return NULL;
703  }
704 
705  _gt[0] = _w[0];
706  _gt[1] = *scale_x;
707 
708  /* check for skew */
709  if (NULL != skew_x && FLT_NEQ(*skew_x, 0))
710  _gt[2] = *skew_x;
711  }
712  /* positive scale-y */
713  if (
714  (NULL != scale_y) &&
715  (*scale_y > 0)
716  ) {
717  if (rt_raster_cell_to_geopoint(
718  rast,
719  0, rast->height,
720  &(_w[0]), &(_w[1]),
721  NULL
722  ) != ES_NONE) {
723  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
724  rt_raster_destroy(rast);
725  _rti_warp_arg_destroy(arg);
726  return NULL;
727  }
728 
729  _gt[3] = _w[1];
730  _gt[5] = *scale_y;
731 
732  /* check for skew */
733  if (NULL != skew_y && FLT_NEQ(*skew_y, 0))
734  _gt[4] = *skew_y;
735  }
736  }
737 
738  rt_raster_destroy(rast);
739  rast = NULL;
740 
741  RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
742  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
743  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
744  _dim[0], _dim[1]);
745 
746  if ( _dim[0] == 0 || _dim[1] == 0 ) {
747  rterror("rt_raster_gdal_warp: The width (%d) or height (%d) of the warped raster is zero", _dim[0], _dim[1]);
748  _rti_warp_arg_destroy(arg);
749  return NULL;
750  }
751 
752  /* load VRT driver */
753  if (!rt_util_gdal_driver_registered("VRT")) {
754  RASTER_DEBUG(3, "Registering VRT driver");
755  GDALRegister_VRT();
756  arg->dst.destroy_drv = 1;
757  }
758  arg->dst.drv = GDALGetDriverByName("VRT");
759  if (NULL == arg->dst.drv) {
760  rterror("rt_raster_gdal_warp: Could not load the output GDAL VRT driver");
761  _rti_warp_arg_destroy(arg);
762  return NULL;
763  }
764 
765  /* create dst dataset */
766  arg->dst.ds = GDALCreate(arg->dst.drv, "", _dim[0], _dim[1], 0, GDT_Byte, dst_options);
767  if (NULL == arg->dst.ds) {
768  rterror("rt_raster_gdal_warp: Could not create GDAL VRT dataset");
769  _rti_warp_arg_destroy(arg);
770  return NULL;
771  }
772 
773  /* set dst srs */
774  if (arg->dst.srs != NULL) {
775  cplerr = GDALSetProjection(arg->dst.ds, arg->dst.srs);
776  if (cplerr != CE_None) {
777  rterror("rt_raster_gdal_warp: Could not set projection");
778  _rti_warp_arg_destroy(arg);
779  return NULL;
780  }
781  RASTER_DEBUGF(3, "Applied SRS: %s", GDALGetProjectionRef(arg->dst.ds));
782  }
783 
784  /* set dst geotransform */
785  cplerr = GDALSetGeoTransform(arg->dst.ds, _gt);
786  if (cplerr != CE_None) {
787  rterror("rt_raster_gdal_warp: Could not set geotransform");
788  _rti_warp_arg_destroy(arg);
789  return NULL;
790  }
791 
792  /* add bands to dst dataset */
793  numBands = rt_raster_get_num_bands(raster);
794  for (i = 0; i < numBands; i++) {
795  rtband = rt_raster_get_band(raster, i);
796  if (NULL == rtband) {
797  rterror("rt_raster_gdal_warp: Could not get band %d for adding to VRT dataset", i);
798  _rti_warp_arg_destroy(arg);
799  return NULL;
800  }
801 
802  pt = rt_band_get_pixtype(rtband);
803  gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
804  if (gdal_pt == GDT_Unknown)
805  rtwarn("rt_raster_gdal_warp: Unknown pixel type for band %d", i);
806 
807  cplerr = GDALAddBand(arg->dst.ds, gdal_pt, NULL);
808  if (cplerr != CE_None) {
809  rterror("rt_raster_gdal_warp: Could not add band to VRT dataset");
810  _rti_warp_arg_destroy(arg);
811  return NULL;
812  }
813 
814  /* get band to write data to */
815  band = NULL;
816  band = GDALGetRasterBand(arg->dst.ds, i + 1);
817  if (NULL == band) {
818  rterror("rt_raster_gdal_warp: Could not get GDAL band for additional processing");
819  _rti_warp_arg_destroy(arg);
820  return NULL;
821  }
822 
823  /* set nodata */
824  if (rt_band_get_hasnodata_flag(rtband) != FALSE) {
825  hasnodata = 1;
826  rt_band_get_nodata(rtband, &nodata);
827  if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
828  rtwarn("rt_raster_gdal_warp: Could not set nodata value for band %d", i);
829  RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
830  }
831  }
832 
833  /* create transformation object */
834  arg->transform.arg.transform = arg->transform.arg.imgproj = GDALCreateGenImgProjTransformer2(
835  arg->src.ds, arg->dst.ds,
836  arg->transform.option.item
837  );
838  if (NULL == arg->transform.arg.transform) {
839  rterror("rt_raster_gdal_warp: Could not create GDAL transformation object");
840  _rti_warp_arg_destroy(arg);
841  return NULL;
842  }
843  arg->transform.func = GDALGenImgProjTransform;
844 
845  /* use approximate transformation object */
846  if (max_err > 0.0) {
847  arg->transform.arg.transform = arg->transform.arg.approx = GDALCreateApproxTransformer(
848  GDALGenImgProjTransform,
849  arg->transform.arg.imgproj, max_err
850  );
851  if (NULL == arg->transform.arg.transform) {
852  rterror("rt_raster_gdal_warp: Could not create GDAL approximate transformation object");
853  _rti_warp_arg_destroy(arg);
854  return NULL;
855  }
856 
857  arg->transform.func = GDALApproxTransform;
858  }
859 
860  /* warp options */
861  arg->wopts = GDALCreateWarpOptions();
862  if (NULL == arg->wopts) {
863  rterror("rt_raster_gdal_warp: Could not create GDAL warp options object");
864  _rti_warp_arg_destroy(arg);
865  return NULL;
866  }
867 
868  /* set options */
869  arg->wopts->eResampleAlg = resample_alg;
870  arg->wopts->hSrcDS = arg->src.ds;
871  arg->wopts->hDstDS = arg->dst.ds;
872  arg->wopts->pfnTransformer = arg->transform.func;
873  arg->wopts->pTransformerArg = arg->transform.arg.transform;
874  arg->wopts->papszWarpOptions = (char **) CPLMalloc(sizeof(char *) * 2);
875  arg->wopts->papszWarpOptions[0] = (char *) CPLMalloc(sizeof(char) * (strlen("INIT_DEST=NO_DATA") + 1));
876  strcpy(arg->wopts->papszWarpOptions[0], "INIT_DEST=NO_DATA");
877  arg->wopts->papszWarpOptions[1] = NULL;
878 
879  /* set band mapping */
880  arg->wopts->nBandCount = numBands;
881  arg->wopts->panSrcBands = (int *) CPLMalloc(sizeof(int) * arg->wopts->nBandCount);
882  arg->wopts->panDstBands = (int *) CPLMalloc(sizeof(int) * arg->wopts->nBandCount);
883  for (i = 0; i < arg->wopts->nBandCount; i++)
884  arg->wopts->panDstBands[i] = arg->wopts->panSrcBands[i] = i + 1;
885 
886  /* set nodata mapping */
887  if (hasnodata) {
888  RASTER_DEBUG(3, "Setting nodata mapping");
889  arg->wopts->padfSrcNoDataReal = (double *) CPLMalloc(numBands * sizeof(double));
890  arg->wopts->padfDstNoDataReal = (double *) CPLMalloc(numBands * sizeof(double));
891  arg->wopts->padfSrcNoDataImag = (double *) CPLMalloc(numBands * sizeof(double));
892  arg->wopts->padfDstNoDataImag = (double *) CPLMalloc(numBands * sizeof(double));
893  if (
894  NULL == arg->wopts->padfSrcNoDataReal ||
895  NULL == arg->wopts->padfDstNoDataReal ||
896  NULL == arg->wopts->padfSrcNoDataImag ||
897  NULL == arg->wopts->padfDstNoDataImag
898  ) {
899  rterror("rt_raster_gdal_warp: Out of memory allocating nodata mapping");
900  _rti_warp_arg_destroy(arg);
901  return NULL;
902  }
903  for (i = 0; i < numBands; i++) {
904  band = rt_raster_get_band(raster, i);
905  if (!band) {
906  rterror("rt_raster_gdal_warp: Could not process bands for nodata values");
907  _rti_warp_arg_destroy(arg);
908  return NULL;
909  }
910 
911  if (!rt_band_get_hasnodata_flag(band)) {
912  /*
913  based on line 1004 of gdalwarp.cpp
914  the problem is that there is a chance that this number is a legitimate value
915  */
916  arg->wopts->padfSrcNoDataReal[i] = -123456.789;
917  }
918  else {
919  rt_band_get_nodata(band, &(arg->wopts->padfSrcNoDataReal[i]));
920  }
921 
922  arg->wopts->padfDstNoDataReal[i] = arg->wopts->padfSrcNoDataReal[i];
923  arg->wopts->padfDstNoDataImag[i] = arg->wopts->padfSrcNoDataImag[i] = 0.0;
924  RASTER_DEBUGF(4, "Mapped nodata value for band %d: %f (%f) => %f (%f)",
925  i,
926  arg->wopts->padfSrcNoDataReal[i], arg->wopts->padfSrcNoDataImag[i],
927  arg->wopts->padfDstNoDataReal[i], arg->wopts->padfDstNoDataImag[i]
928  );
929  }
930  }
931 
932  /* warp raster */
933  RASTER_DEBUG(3, "Warping raster");
934  cplerr = GDALInitializeWarpedVRT(arg->dst.ds, arg->wopts);
935  if (cplerr != CE_None) {
936  rterror("rt_raster_gdal_warp: Could not warp raster");
937  _rti_warp_arg_destroy(arg);
938  return NULL;
939  }
940  /*
941  GDALSetDescription(arg->dst.ds, "/tmp/warped.vrt");
942  */
943  GDALFlushCache(arg->dst.ds);
944  RASTER_DEBUG(3, "Raster warped");
945 
946  /* convert gdal dataset to raster */
947  RASTER_DEBUG(3, "Converting GDAL dataset to raster");
948  rast = rt_raster_from_gdal_dataset(arg->dst.ds);
949 
950  _rti_warp_arg_destroy(arg);
951 
952  if (NULL == rast) {
953  rterror("rt_raster_gdal_warp: Could not warp raster");
954  return NULL;
955  }
956 
957  /* substitute spatial, reset back to default */
958  if (subspatial) {
959  double gt[6] = {0, 1, 0, 0, 0, -1};
960  /* See http://trac.osgeo.org/postgis/ticket/2911 */
961  /* We should proably also tweak rotation here */
962  /* NOTE: the times 10 is because it was divided by 10 in a section above,
963  * I'm not sure the above division was needed */
964  gt[1] = _scale[0] * 10;
965  gt[5] = -1 * _scale[1] * 10;
966 
967  rt_raster_set_geotransform_matrix(rast, gt);
968  rt_raster_set_srid(rast, SRID_UNKNOWN);
969  }
970 
971  RASTER_DEBUG(3, "done");
972 
973  return rast;
974 }
rt_envelope::MinY
double MinY
Definition: librtcore.h:167
rt_envelope::UpperLeftX
double UpperLeftX
Definition: librtcore.h:170
rt_band_t
Definition: librtcore.h:2264
rt_raster_get_num_bands
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
_rti_warp_arg_t::wopts
GDALWarpOptions * wopts
Definition: rt_warp.c:51
rt_raster_set_geotransform_matrix
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster&#39;s geotransform using 6-element array.
Definition: rt_raster.c:727
rt_envelope::MaxY
double MaxY
Definition: librtcore.h:168
rt_raster_get_geotransform_matrix
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:706
rt_envelope::MaxX
double MaxX
Definition: librtcore.h:166
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:168
ovdump.band
band
Definition: ovdump.py:57
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:1809
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2185
_rti_warp_arg_t::srs
char * srs
Definition: rt_warp.c:47
rterror
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
_rti_warp_arg_t::destroy_drv
int destroy_drv
Definition: rt_warp.c:48
_rti_warp_arg_t::ds
GDALDatasetH ds
Definition: rt_warp.c:46
rt_raster_t::ipY
double ipY
Definition: librtcore.h:2247
window.gt
gt
Definition: window.py:77
rt_pixtype
rt_pixtype
Definition: librtcore.h:185
rt_raster_serialized_t::width
uint16_t width
Definition: librtcore.h:2227
rt_raster_t::height
uint16_t height
Definition: librtcore.h:2253
_rti_warp_arg_init
static _rti_warp_arg _rti_warp_arg_init()
Definition: rt_warp.c:71
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
rt_raster_cell_to_geopoint
rt_errorstate rt_raster_cell_to_geopoint(rt_raster raster, double xr, double yr, double *xw, double *yw, double *gt)
Convert an xr, yr raster point to an xw, yw point on map.
Definition: rt_raster.c:755
rt_util_pixtype_to_gdal_datatype
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_util.c:118
rt_util_gdal_driver_registered
int rt_util_gdal_driver_registered(const char *drv)
Definition: rt_util.c:353
rtwarn
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
rt_envelope
Definition: librtcore.h:164
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:137
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:199
rt_raster_t::width
uint16_t width
Definition: librtcore.h:2252
SRID_UNKNOWN
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:188
_rti_warp_arg_destroy
static void _rti_warp_arg_destroy(_rti_warp_arg arg)
Definition: rt_warp.c:105
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:381
FLT_NEQ
#define FLT_NEQ(x, y)
Definition: librtcore.h:2184
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
rt_raster_set_srid
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster&#39;s SRID.
Definition: rt_raster.c:363
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
rt_raster_serialized_t::numBands
uint16_t numBands
Definition: librtcore.h:2205
rt_raster_get_srid
int32_t rt_raster_get_srid(rt_raster raster)
Get raster&#39;s SRID.
Definition: rt_raster.c:356
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
rt_raster_geopoint_to_cell
rt_errorstate rt_raster_geopoint_to_cell(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:806
_rti_warp_arg_t
Definition: rt_warp.c:42
rt_raster_new
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
rt_raster_serialized_t::height
uint16_t height
Definition: librtcore.h:2228
FALSE
#define FALSE
Definition: dbfopen.c:168
rt_envelope::MinX
double MinX
Definition: librtcore.h:165
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
_rti_warp_arg_t::src
struct _rti_warp_arg_t::@12 src
rt_util_gdal_convert_sr
char * rt_util_gdal_convert_sr(const char *srs, int proj4)
Definition: rt_util.c:212
rt_band_get_pixtype
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:498
_rti_warp_arg_t::drv
GDALDriverH drv
Definition: rt_warp.c:45
_rti_warp_arg_t::dst
struct _rti_warp_arg_t::@12 dst
rt_raster_compute_skewed_raster
rt_raster rt_raster_compute_skewed_raster(rt_envelope extent, double *skew, double *scale, double tolerance)
Definition: rt_raster.c:960
_rti_warp_arg_t::transform
void * transform
Definition: rt_warp.c:60
rt_raster_t::ipX
double ipX
Definition: librtcore.h:2246
rt_raster_t
Definition: librtcore.h:2235
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:2165
rt_raster_t::scaleX
double scaleX
Definition: librtcore.h:2244
rt_envelope::UpperLeftY
double UpperLeftY
Definition: librtcore.h:171
if
if(!(yy_init))
Definition: lwin_wkt_lex.c:1080
rt_raster_t::scaleY
double scaleY
Definition: librtcore.h:2245
Here is the call graph for this function:
Here is the caller graph for this function: