PostGIS  2.1.10dev-r@@SVN_REVISION@@
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 9710 of file rt_api.c.

References _rti_warp_arg_destroy(), _rti_warp_arg_init(), 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_to_gdal_mem(), rt_util_gdal_convert_sr(), rt_util_gdal_driver_registered(), rt_util_pixtype_to_gdal_datatype(), rtalloc(), rterror(), rtinfo(), 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().

9719  {
9720  CPLErr cplerr;
9721  char *dst_options[] = {"SUBCLASS=VRTWarpedDataset", NULL};
9722  _rti_warp_arg arg = NULL;
9723 
9724  int hasnodata = 0;
9725 
9726  GDALRasterBandH band;
9727  rt_band rtband = NULL;
9728  rt_pixtype pt = PT_END;
9729  GDALDataType gdal_pt = GDT_Unknown;
9730  double nodata = 0.0;
9731 
9732  double _gt[6] = {0};
9733  double dst_extent[4];
9734  rt_envelope extent;
9735 
9736  int _dim[2] = {0};
9737  double _skew[2] = {0};
9738  double _scale[2] = {0};
9739  int ul_user = 0;
9740 
9741  rt_raster rast = NULL;
9742  int i = 0;
9743  int numBands = 0;
9744 
9745  int subgt = 0;
9746 
9747  RASTER_DEBUG(3, "starting");
9748 
9749  assert(NULL != raster);
9750 
9751  /* internal variables */
9752  arg = _rti_warp_arg_init();
9753  if (arg == NULL) {
9754  rterror("rt_raster_gdal_warp: Could not initialize internal variables");
9755  return NULL;
9756  }
9757 
9758  /*
9759  max_err must be gte zero
9760 
9761  the value 0.125 is the default used in gdalwarp.cpp on line 283
9762  */
9763  if (max_err < 0.) max_err = 0.125;
9764  RASTER_DEBUGF(4, "max_err = %f", max_err);
9765 
9766  /* handle srs */
9767  if (src_srs != NULL) {
9768  /* reprojection taking place */
9769  if (dst_srs != NULL && strcmp(src_srs, dst_srs) != 0) {
9770  RASTER_DEBUG(4, "Warp operation does include a reprojection");
9771  arg->src.srs = rt_util_gdal_convert_sr(src_srs, 0);
9772  arg->dst.srs = rt_util_gdal_convert_sr(dst_srs, 0);
9773 
9774  if (arg->src.srs == NULL || arg->dst.srs == NULL) {
9775  rterror("rt_raster_gdal_warp: Could not convert srs values to GDAL accepted format");
9776  _rti_warp_arg_destroy(arg);
9777  return NULL;
9778  }
9779  }
9780  /* no reprojection, a stub just for clarity */
9781  else {
9782  RASTER_DEBUG(4, "Warp operation does NOT include reprojection");
9783  }
9784  }
9785  else if (dst_srs != NULL) {
9786  /* dst_srs provided but not src_srs */
9787  rterror("rt_raster_gdal_warp: SRS required for input raster if SRS provided for warped raster");
9788  _rti_warp_arg_destroy(arg);
9789  return NULL;
9790  }
9791 
9792  /* load raster into a GDAL MEM dataset */
9793  arg->src.ds = rt_raster_to_gdal_mem(raster, arg->src.srs, NULL, NULL, 0, &(arg->src.drv), &(arg->src.destroy_drv));
9794  if (NULL == arg->src.ds) {
9795  rterror("rt_raster_gdal_warp: Could not convert raster to GDAL MEM format");
9796  _rti_warp_arg_destroy(arg);
9797  return NULL;
9798  }
9799  RASTER_DEBUG(3, "raster loaded into GDAL MEM dataset");
9800 
9801  /* special case when src_srs and dst_srs is NULL and raster's geotransform matrix is default */
9802  if (src_srs == NULL && dst_srs == NULL && rt_raster_get_srid(raster) == SRID_UNKNOWN) {
9803  double gt[6];
9804 
9805 #if POSTGIS_DEBUG_LEVEL > 3
9806  GDALGetGeoTransform(arg->src.ds, gt);
9807  RASTER_DEBUGF(3, "GDAL MEM geotransform: %f, %f, %f, %f, %f, %f",
9808  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
9809 #endif
9810 
9811  /* default geotransform */
9813  RASTER_DEBUGF(3, "raster geotransform: %f, %f, %f, %f, %f, %f",
9814  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
9815 
9816  if (
9817  FLT_EQ(gt[0], 0) &&
9818  FLT_EQ(gt[1], 1) &&
9819  FLT_EQ(gt[2], 0) &&
9820  FLT_EQ(gt[3], 0) &&
9821  FLT_EQ(gt[4], 0) &&
9822  FLT_EQ(gt[5], -1)
9823  ) {
9824  double ngt[6] = {0, 10, 0, 0, 0, -10};
9825 
9826  rtinfo("Raster has default geotransform. Adjusting metadata for use of GDAL Warp API");
9827 
9828  GDALSetGeoTransform(arg->src.ds, ngt);
9829  GDALFlushCache(arg->src.ds);
9830 
9831  subgt = 1;
9832 
9833 #if POSTGIS_DEBUG_LEVEL > 3
9834  GDALGetGeoTransform(arg->src.ds, gt);
9835  RASTER_DEBUGF(3, "GDAL MEM geotransform: %f, %f, %f, %f, %f, %f",
9836  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
9837 #endif
9838  }
9839  }
9840 
9841  /* set transform options */
9842  if (arg->src.srs != NULL || arg->dst.srs != NULL) {
9843  arg->transform.option.len = 2;
9844  arg->transform.option.item = rtalloc(sizeof(char *) * (arg->transform.option.len + 1));
9845  if (NULL == arg->transform.option.item) {
9846  rterror("rt_raster_gdal_warp: Could not allocation memory for transform options");
9847  _rti_warp_arg_destroy(arg);
9848  return NULL;
9849  }
9850  memset(arg->transform.option.item, 0, sizeof(char *) * (arg->transform.option.len + 1));
9851 
9852  for (i = 0; i < arg->transform.option.len; i++) {
9853  switch (i) {
9854  case 1:
9855  if (arg->dst.srs != NULL)
9856  arg->transform.option.item[i] = (char *) rtalloc(sizeof(char) * (strlen("DST_SRS=") + strlen(arg->dst.srs) + 1));
9857  else
9858  arg->transform.option.item[i] = (char *) rtalloc(sizeof(char) * (strlen("DST_SRS=") + 1));
9859  break;
9860  case 0:
9861  if (arg->src.srs != NULL)
9862  arg->transform.option.item[i] = (char *) rtalloc(sizeof(char) * (strlen("SRC_SRS=") + strlen(arg->src.srs) + 1));
9863  else
9864  arg->transform.option.item[i] = (char *) rtalloc(sizeof(char) * (strlen("SRC_SRS=") + 1));
9865  break;
9866  }
9867  if (NULL == arg->transform.option.item[i]) {
9868  rterror("rt_raster_gdal_warp: Could not allocation memory for transform options");
9869  _rti_warp_arg_destroy(arg);
9870  return NULL;
9871  }
9872 
9873  switch (i) {
9874  case 1:
9875  if (arg->dst.srs != NULL) {
9876  snprintf(
9877  arg->transform.option.item[i],
9878  sizeof(char) * (strlen("DST_SRS=") + strlen(arg->dst.srs) + 1),
9879  "DST_SRS=%s",
9880  arg->dst.srs
9881  );
9882  }
9883  else
9884  sprintf(arg->transform.option.item[i], "%s", "DST_SRS=");
9885  break;
9886  case 0:
9887  if (arg->src.srs != NULL) {
9888  snprintf(
9889  arg->transform.option.item[i],
9890  sizeof(char) * (strlen("SRC_SRS=") + strlen(arg->src.srs) + 1),
9891  "SRC_SRS=%s",
9892  arg->src.srs
9893  );
9894  }
9895  else
9896  sprintf(arg->transform.option.item[i], "%s", "SRC_SRS=");
9897  break;
9898  }
9899  RASTER_DEBUGF(4, "arg->transform.option.item[%d] = %s", i, arg->transform.option.item[i]);
9900  }
9901  }
9902  else
9903  arg->transform.option.len = 0;
9904 
9905  /* transformation object for building dst dataset */
9906  arg->transform.arg.transform = GDALCreateGenImgProjTransformer2(arg->src.ds, NULL, arg->transform.option.item);
9907  if (NULL == arg->transform.arg.transform) {
9908  rterror("rt_raster_gdal_warp: Could not create GDAL transformation object for output dataset creation");
9909  _rti_warp_arg_destroy(arg);
9910  return NULL;
9911  }
9912 
9913  /* get approximate output georeferenced bounds and resolution */
9914  cplerr = GDALSuggestedWarpOutput2(
9915  arg->src.ds, GDALGenImgProjTransform,
9916  arg->transform.arg.transform, _gt, &(_dim[0]), &(_dim[1]), dst_extent, 0);
9917  if (cplerr != CE_None) {
9918  rterror("rt_raster_gdal_warp: Could not get GDAL suggested warp output for output dataset creation");
9919  _rti_warp_arg_destroy(arg);
9920  return NULL;
9921  }
9922  GDALDestroyGenImgProjTransformer(arg->transform.arg.transform);
9923  arg->transform.arg.transform = NULL;
9924 
9925  /*
9926  don't use suggested dimensions as use of suggested scales
9927  on suggested extent will result in suggested dimensions
9928  */
9929  _dim[0] = 0;
9930  _dim[1] = 0;
9931 
9932  RASTER_DEBUGF(3, "Suggested geotransform: %f, %f, %f, %f, %f, %f",
9933  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
9934 
9935  /* store extent in easier-to-use object */
9936  extent.MinX = dst_extent[0];
9937  extent.MinY = dst_extent[1];
9938  extent.MaxX = dst_extent[2];
9939  extent.MaxY = dst_extent[3];
9940 
9941  extent.UpperLeftX = dst_extent[0];
9942  extent.UpperLeftY = dst_extent[3];
9943 
9944  RASTER_DEBUGF(3, "Suggested extent: %f, %f, %f, %f",
9945  extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
9946 
9947  /* scale and width/height are mutually exclusive */
9948  if (
9949  ((NULL != scale_x) || (NULL != scale_y)) &&
9950  ((NULL != width) || (NULL != height))
9951  ) {
9952  rterror("rt_raster_gdal_warp: Scale X/Y and width/height are mutually exclusive. Only provide one");
9953  _rti_warp_arg_destroy(arg);
9954  return NULL;
9955  }
9956 
9957  /* user-defined width */
9958  if (NULL != width) {
9959  _dim[0] = abs(*width);
9960  _scale[0] = fabs((extent.MaxX - extent.MinX) / ((double) _dim[0]));
9961  }
9962  /* user-defined height */
9963  if (NULL != height) {
9964  _dim[1] = abs(*height);
9965  _scale[1] = fabs((extent.MaxY - extent.MinY) / ((double) _dim[1]));
9966  }
9967 
9968  /* user-defined scale */
9969  if (
9970  ((NULL != scale_x) && (FLT_NEQ(*scale_x, 0.0))) &&
9971  ((NULL != scale_y) && (FLT_NEQ(*scale_y, 0.0)))
9972  ) {
9973  _scale[0] = fabs(*scale_x);
9974  _scale[1] = fabs(*scale_y);
9975 
9976  /* special override */
9977  if (subgt) {
9978  _scale[0] *= 10;
9979  _scale[1] *= 10;
9980  }
9981  }
9982  else if (
9983  ((NULL != scale_x) && (NULL == scale_y)) ||
9984  ((NULL == scale_x) && (NULL != scale_y))
9985  ) {
9986  rterror("rt_raster_gdal_warp: Both X and Y scale values must be provided for scale");
9987  _rti_warp_arg_destroy(arg);
9988  return NULL;
9989  }
9990 
9991  /* scale not defined, use suggested */
9992  if (FLT_EQ(_scale[0], 0) && FLT_EQ(_scale[1], 0)) {
9993  _scale[0] = fabs(_gt[1]);
9994  _scale[1] = fabs(_gt[5]);
9995  }
9996 
9997  RASTER_DEBUGF(4, "Using scale: %f x %f", _scale[0], -1 * _scale[1]);
9998 
9999  /* user-defined skew */
10000  if (NULL != skew_x) {
10001  _skew[0] = *skew_x;
10002 
10003  /*
10004  negative scale-x affects skew
10005  for now, force skew to be in left-right, top-down orientation
10006  */
10007  if (
10008  NULL != scale_x &&
10009  *scale_x < 0.
10010  ) {
10011  _skew[0] *= -1;
10012  }
10013  }
10014  if (NULL != skew_y) {
10015  _skew[1] = *skew_y;
10016 
10017  /*
10018  positive scale-y affects skew
10019  for now, force skew to be in left-right, top-down orientation
10020  */
10021  if (
10022  NULL != scale_y &&
10023  *scale_y > 0.
10024  ) {
10025  _skew[1] *= -1;
10026  }
10027  }
10028 
10029  RASTER_DEBUGF(4, "Using skew: %f x %f", _skew[0], _skew[1]);
10030 
10031  /* reprocess extent if skewed */
10032  if (
10033  FLT_NEQ(_skew[0], 0) ||
10034  FLT_NEQ(_skew[1], 0)
10035  ) {
10036  rt_raster skewedrast;
10037 
10038  RASTER_DEBUG(3, "Computing skewed extent's envelope");
10039 
10040  skewedrast = rt_raster_compute_skewed_raster(
10041  extent,
10042  _skew,
10043  _scale,
10044  0.01
10045  );
10046  if (skewedrast == NULL) {
10047  rterror("rt_raster_gdal_warp: Could not compute skewed raster");
10048  _rti_warp_arg_destroy(arg);
10049  return NULL;
10050  }
10051 
10052  if (_dim[0] == 0)
10053  _dim[0] = skewedrast->width;
10054  if (_dim[1] == 0)
10055  _dim[1] = skewedrast->height;
10056 
10057  extent.UpperLeftX = skewedrast->ipX;
10058  extent.UpperLeftY = skewedrast->ipY;
10059 
10060  rt_raster_destroy(skewedrast);
10061  }
10062 
10063  /* dimensions not defined, compute */
10064  if (!_dim[0])
10065  _dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
10066  if (!_dim[1])
10067  _dim[1] = (int) fmax((fabs(extent.MaxY - extent.MinY) + (_scale[1] / 2.)) / _scale[1], 1);
10068 
10069  /* temporary raster */
10070  rast = rt_raster_new(_dim[0], _dim[1]);
10071  if (rast == NULL) {
10072  rterror("rt_raster_gdal_warp: Out of memory allocating temporary raster");
10073  _rti_warp_arg_destroy(arg);
10074  return NULL;
10075  }
10076 
10077  /* set raster's spatial attributes */
10078  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
10079  rt_raster_set_scale(rast, _scale[0], -1 * _scale[1]);
10080  rt_raster_set_skews(rast, _skew[0], _skew[1]);
10081 
10083  RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
10084  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
10085  RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
10086  _dim[0], _dim[1]);
10087 
10088  /* user-defined upper-left corner */
10089  if (
10090  NULL != ul_xw &&
10091  NULL != ul_yw
10092  ) {
10093  ul_user = 1;
10094 
10095  RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
10096 
10097  /* set upper-left corner */
10098  rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
10099  extent.UpperLeftX = *ul_xw;
10100  extent.UpperLeftY = *ul_yw;
10101  }
10102  else if (
10103  ((NULL != ul_xw) && (NULL == ul_yw)) ||
10104  ((NULL == ul_xw) && (NULL != ul_yw))
10105  ) {
10106  rterror("rt_raster_gdal_warp: Both X and Y upper-left corner values must be provided");
10107  rt_raster_destroy(rast);
10108  _rti_warp_arg_destroy(arg);
10109  return NULL;
10110  }
10111 
10112  /* alignment only considered if upper-left corner not provided */
10113  if (
10114  !ul_user && (
10115  (NULL != grid_xw) || (NULL != grid_yw)
10116  )
10117  ) {
10118 
10119  if (
10120  ((NULL != grid_xw) && (NULL == grid_yw)) ||
10121  ((NULL == grid_xw) && (NULL != grid_yw))
10122  ) {
10123  rterror("rt_raster_gdal_warp: Both X and Y alignment values must be provided");
10124  rt_raster_destroy(rast);
10125  _rti_warp_arg_destroy(arg);
10126  return NULL;
10127  }
10128 
10129  RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw);
10130 
10131  do {
10132  double _r[2] = {0};
10133  double _w[2] = {0};
10134 
10135  /* raster is already aligned */
10136  if (FLT_EQ(*grid_xw, extent.UpperLeftX) && FLT_EQ(*grid_yw, extent.UpperLeftY)) {
10137  RASTER_DEBUG(3, "Skipping raster alignment as it is already aligned to grid");
10138  break;
10139  }
10140 
10141  extent.UpperLeftX = rast->ipX;
10142  extent.UpperLeftY = rast->ipY;
10143  rt_raster_set_offsets(rast, *grid_xw, *grid_yw);
10144 
10145  /* process upper-left corner */
10147  rast,
10148  extent.UpperLeftX, extent.UpperLeftY,
10149  &(_r[0]), &(_r[1]),
10150  NULL
10151  ) != ES_NONE) {
10152  rterror("rt_raster_gdal_warp: Could not compute raster pixel for spatial coordinates");
10153  rt_raster_destroy(rast);
10154  _rti_warp_arg_destroy(arg);
10155  return NULL;
10156  }
10157 
10159  rast,
10160  _r[0], _r[1],
10161  &(_w[0]), &(_w[1]),
10162  NULL
10163  ) != ES_NONE) {
10164  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
10165 
10166  rt_raster_destroy(rast);
10167  _rti_warp_arg_destroy(arg);
10168  return NULL;
10169  }
10170 
10171  /* shift occurred */
10172  if (FLT_NEQ(_w[0], extent.UpperLeftX)) {
10173  if (NULL == width)
10174  rast->width++;
10175  else if (NULL == scale_x) {
10176  double _c[2] = {0};
10177 
10178  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
10179 
10180  /* get upper-right corner */
10182  rast,
10183  rast->width, 0,
10184  &(_c[0]), &(_c[1]),
10185  NULL
10186  ) != ES_NONE) {
10187  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
10188  rt_raster_destroy(rast);
10189  _rti_warp_arg_destroy(arg);
10190  return NULL;
10191  }
10192 
10193  rast->scaleX = fabs((_c[0] - _w[0]) / ((double) rast->width));
10194  }
10195  }
10196  if (FLT_NEQ(_w[1], extent.UpperLeftY)) {
10197  if (NULL == height)
10198  rast->height++;
10199  else if (NULL == scale_y) {
10200  double _c[2] = {0};
10201 
10202  rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
10203 
10204  /* get upper-right corner */
10206  rast,
10207  0, rast->height,
10208  &(_c[0]), &(_c[1]),
10209  NULL
10210  ) != ES_NONE) {
10211  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
10212 
10213  rt_raster_destroy(rast);
10214  _rti_warp_arg_destroy(arg);
10215  return NULL;
10216  }
10217 
10218  rast->scaleY = -1 * fabs((_c[1] - _w[1]) / ((double) rast->height));
10219  }
10220  }
10221 
10222  rt_raster_set_offsets(rast, _w[0], _w[1]);
10223  RASTER_DEBUGF(4, "aligned offsets: %f x %f", _w[0], _w[1]);
10224  }
10225  while (0);
10226  }
10227 
10228  /*
10229  after this point, rt_envelope extent is no longer used
10230  */
10231 
10232  /* get key attributes from rast */
10233  _dim[0] = rast->width;
10234  _dim[1] = rast->height;
10236 
10237  /* scale-x is negative or scale-y is positive */
10238  if ((
10239  (NULL != scale_x) && (*scale_x < 0.)
10240  ) || (
10241  (NULL != scale_y) && (*scale_y > 0)
10242  )) {
10243  double _w[2] = {0};
10244 
10245  /* negative scale-x */
10246  if (
10247  (NULL != scale_x) &&
10248  (*scale_x < 0.)
10249  ) {
10251  rast,
10252  rast->width, 0,
10253  &(_w[0]), &(_w[1]),
10254  NULL
10255  ) != ES_NONE) {
10256  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
10257  rt_raster_destroy(rast);
10258  _rti_warp_arg_destroy(arg);
10259  return NULL;
10260  }
10261 
10262  _gt[0] = _w[0];
10263  _gt[1] = *scale_x;
10264 
10265  /* check for skew */
10266  if (NULL != skew_x && FLT_NEQ(*skew_x, 0))
10267  _gt[2] = *skew_x;
10268  }
10269  /* positive scale-y */
10270  if (
10271  (NULL != scale_y) &&
10272  (*scale_y > 0)
10273  ) {
10275  rast,
10276  0, rast->height,
10277  &(_w[0]), &(_w[1]),
10278  NULL
10279  ) != ES_NONE) {
10280  rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
10281  rt_raster_destroy(rast);
10282  _rti_warp_arg_destroy(arg);
10283  return NULL;
10284  }
10285 
10286  _gt[3] = _w[1];
10287  _gt[5] = *scale_y;
10288 
10289  /* check for skew */
10290  if (NULL != skew_y && FLT_NEQ(*skew_y, 0))
10291  _gt[4] = *skew_y;
10292  }
10293  }
10294 
10295  rt_raster_destroy(rast);
10296  rast = NULL;
10297 
10298  RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
10299  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
10300  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
10301  _dim[0], _dim[1]);
10302 
10303  if (_dim[0] == 0 || _dim[1] == 0) {
10304  rterror("rt_raster_gdal_warp: The width (%f) or height (%f) of the warped raster is zero", _dim[0], _dim[1]);
10305  _rti_warp_arg_destroy(arg);
10306  return NULL;
10307  }
10308 
10309  /* load VRT driver */
10310  if (!rt_util_gdal_driver_registered("VRT")) {
10311  GDALRegister_VRT();
10312  arg->dst.destroy_drv = 1;
10313  }
10314  arg->dst.drv = GDALGetDriverByName("VRT");
10315  if (NULL == arg->dst.drv) {
10316  rterror("rt_raster_gdal_warp: Could not load the output GDAL VRT driver");
10317  _rti_warp_arg_destroy(arg);
10318  return NULL;
10319  }
10320 
10321  /* create dst dataset */
10322  arg->dst.ds = GDALCreate(arg->dst.drv, "", _dim[0], _dim[1], 0, GDT_Byte, dst_options);
10323  if (NULL == arg->dst.ds) {
10324  rterror("rt_raster_gdal_warp: Could not create GDAL VRT dataset");
10325  _rti_warp_arg_destroy(arg);
10326  return NULL;
10327  }
10328 
10329  /* set dst srs */
10330  if (arg->dst.srs != NULL) {
10331  cplerr = GDALSetProjection(arg->dst.ds, arg->dst.srs);
10332  if (cplerr != CE_None) {
10333  rterror("rt_raster_gdal_warp: Could not set projection");
10334  _rti_warp_arg_destroy(arg);
10335  return NULL;
10336  }
10337  RASTER_DEBUGF(3, "Applied SRS: %s", GDALGetProjectionRef(arg->dst.ds));
10338  }
10339 
10340  /* set dst geotransform */
10341  cplerr = GDALSetGeoTransform(arg->dst.ds, _gt);
10342  if (cplerr != CE_None) {
10343  rterror("rt_raster_gdal_warp: Could not set geotransform");
10344  _rti_warp_arg_destroy(arg);
10345  return NULL;
10346  }
10347 
10348  /* add bands to dst dataset */
10349  numBands = rt_raster_get_num_bands(raster);
10350  for (i = 0; i < numBands; i++) {
10351  rtband = rt_raster_get_band(raster, i);
10352  if (NULL == rtband) {
10353  rterror("rt_raster_gdal_warp: Could not get band %d for adding to VRT dataset", i);
10354  _rti_warp_arg_destroy(arg);
10355  return NULL;
10356  }
10357 
10358  pt = rt_band_get_pixtype(rtband);
10359  gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
10360  if (gdal_pt == GDT_Unknown)
10361  rtwarn("rt_raster_gdal_warp: Unknown pixel type for band %d", i);
10362 
10363  cplerr = GDALAddBand(arg->dst.ds, gdal_pt, NULL);
10364  if (cplerr != CE_None) {
10365  rterror("rt_raster_gdal_warp: Could not add band to VRT dataset");
10366  _rti_warp_arg_destroy(arg);
10367  return NULL;
10368  }
10369 
10370  /* get band to write data to */
10371  band = NULL;
10372  band = GDALGetRasterBand(arg->dst.ds, i + 1);
10373  if (NULL == band) {
10374  rterror("rt_raster_gdal_warp: Could not get GDAL band for additional processing");
10375  _rti_warp_arg_destroy(arg);
10376  return NULL;
10377  }
10378 
10379  /* set nodata */
10380  if (rt_band_get_hasnodata_flag(rtband) != FALSE) {
10381  hasnodata = 1;
10382  rt_band_get_nodata(rtband, &nodata);
10383  if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
10384  rtwarn("rt_raster_gdal_warp: Could not set nodata value for band %d", i);
10385  RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
10386  }
10387  }
10388 
10389  /* create transformation object */
10390  arg->transform.arg.transform = arg->transform.arg.imgproj = GDALCreateGenImgProjTransformer2(
10391  arg->src.ds, arg->dst.ds,
10392  arg->transform.option.item
10393  );
10394  if (NULL == arg->transform.arg.transform) {
10395  rterror("rt_raster_gdal_warp: Could not create GDAL transformation object");
10396  _rti_warp_arg_destroy(arg);
10397  return NULL;
10398  }
10399  arg->transform.func = GDALGenImgProjTransform;
10400 
10401  /* use approximate transformation object */
10402  if (max_err > 0.0) {
10403  arg->transform.arg.transform = arg->transform.arg.approx = GDALCreateApproxTransformer(
10404  GDALGenImgProjTransform,
10405  arg->transform.arg.imgproj, max_err
10406  );
10407  if (NULL == arg->transform.arg.transform) {
10408  rterror("rt_raster_gdal_warp: Could not create GDAL approximate transformation object");
10409  _rti_warp_arg_destroy(arg);
10410  return NULL;
10411  }
10412 
10413  arg->transform.func = GDALApproxTransform;
10414  }
10415 
10416  /* warp options */
10417  arg->wopts = GDALCreateWarpOptions();
10418  if (NULL == arg->wopts) {
10419  rterror("rt_raster_gdal_warp: Could not create GDAL warp options object");
10420  _rti_warp_arg_destroy(arg);
10421  return NULL;
10422  }
10423 
10424  /* set options */
10425  arg->wopts->eResampleAlg = resample_alg;
10426  arg->wopts->hSrcDS = arg->src.ds;
10427  arg->wopts->hDstDS = arg->dst.ds;
10428  arg->wopts->pfnTransformer = arg->transform.func;
10429  arg->wopts->pTransformerArg = arg->transform.arg.transform;
10430  arg->wopts->papszWarpOptions = (char **) CPLMalloc(sizeof(char *) * 2);
10431  arg->wopts->papszWarpOptions[0] = (char *) CPLMalloc(sizeof(char) * (strlen("INIT_DEST=NO_DATA") + 1));
10432  strcpy(arg->wopts->papszWarpOptions[0], "INIT_DEST=NO_DATA");
10433  arg->wopts->papszWarpOptions[1] = NULL;
10434 
10435  /* set band mapping */
10436  arg->wopts->nBandCount = numBands;
10437  arg->wopts->panSrcBands = (int *) CPLMalloc(sizeof(int) * arg->wopts->nBandCount);
10438  arg->wopts->panDstBands = (int *) CPLMalloc(sizeof(int) * arg->wopts->nBandCount);
10439  for (i = 0; i < arg->wopts->nBandCount; i++)
10440  arg->wopts->panDstBands[i] = arg->wopts->panSrcBands[i] = i + 1;
10441 
10442  /* set nodata mapping */
10443  if (hasnodata) {
10444  RASTER_DEBUG(3, "Setting nodata mapping");
10445  arg->wopts->padfSrcNoDataReal = (double *) CPLMalloc(numBands * sizeof(double));
10446  arg->wopts->padfDstNoDataReal = (double *) CPLMalloc(numBands * sizeof(double));
10447  arg->wopts->padfSrcNoDataImag = (double *) CPLMalloc(numBands * sizeof(double));
10448  arg->wopts->padfDstNoDataImag = (double *) CPLMalloc(numBands * sizeof(double));
10449  if (
10450  NULL == arg->wopts->padfSrcNoDataReal ||
10451  NULL == arg->wopts->padfDstNoDataReal ||
10452  NULL == arg->wopts->padfSrcNoDataImag ||
10453  NULL == arg->wopts->padfDstNoDataImag
10454  ) {
10455  rterror("rt_raster_gdal_warp: Out of memory allocating nodata mapping");
10456  _rti_warp_arg_destroy(arg);
10457  return NULL;
10458  }
10459  for (i = 0; i < numBands; i++) {
10460  band = rt_raster_get_band(raster, i);
10461  if (!band) {
10462  rterror("rt_raster_gdal_warp: Could not process bands for nodata values");
10463  _rti_warp_arg_destroy(arg);
10464  return NULL;
10465  }
10466 
10467  if (!rt_band_get_hasnodata_flag(band)) {
10468  /*
10469  based on line 1004 of gdalwarp.cpp
10470  the problem is that there is a chance that this number is a legitimate value
10471  */
10472  arg->wopts->padfSrcNoDataReal[i] = -123456.789;
10473  }
10474  else {
10475  rt_band_get_nodata(band, &(arg->wopts->padfSrcNoDataReal[i]));
10476  }
10477 
10478  arg->wopts->padfDstNoDataReal[i] = arg->wopts->padfSrcNoDataReal[i];
10479  arg->wopts->padfDstNoDataImag[i] = arg->wopts->padfSrcNoDataImag[i] = 0.0;
10480  RASTER_DEBUGF(4, "Mapped nodata value for band %d: %f (%f) => %f (%f)",
10481  i,
10482  arg->wopts->padfSrcNoDataReal[i], arg->wopts->padfSrcNoDataImag[i],
10483  arg->wopts->padfDstNoDataReal[i], arg->wopts->padfDstNoDataImag[i]
10484  );
10485  }
10486  }
10487 
10488  /* warp raster */
10489  RASTER_DEBUG(3, "Warping raster");
10490  cplerr = GDALInitializeWarpedVRT(arg->dst.ds, arg->wopts);
10491  if (cplerr != CE_None) {
10492  rterror("rt_raster_gdal_warp: Could not warp raster");
10493  _rti_warp_arg_destroy(arg);
10494  return NULL;
10495  }
10496  /*
10497  GDALSetDescription(arg->dst.ds, "/tmp/warped.vrt");
10498  */
10499  GDALFlushCache(arg->dst.ds);
10500  RASTER_DEBUG(3, "Raster warped");
10501 
10502  /* convert gdal dataset to raster */
10503  RASTER_DEBUG(3, "Converting GDAL dataset to raster");
10504  rast = rt_raster_from_gdal_dataset(arg->dst.ds);
10505 
10506  _rti_warp_arg_destroy(arg);
10507 
10508  if (NULL == rast) {
10509  rterror("rt_raster_gdal_warp: Could not warp raster");
10510  return NULL;
10511  }
10512 
10513  /* substitute geotransform matrix, reset back to default */
10514  if (subgt) {
10515  double gt[6] = {0, 1, 0, 0, 0, -1};
10516  /* See http://trac.osgeo.org/postgis/ticket/2911 */
10517  /* We should proably also tweak rotation here */
10518  /* NOTE: the division by 10 is because it was multiplied by 10 in
10519  * a section above. I'm not sure the above division was needed */
10520  gt[1] = _scale[0] / 10;
10521  gt[5] = -1 * _scale[1] / 10;
10522 
10524  }
10525 
10526  RASTER_DEBUG(3, "done");
10527 
10528  return rast;
10529 }
rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds)
Return a raster from a GDAL dataset.
Definition: rt_api.c:9339
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_api.c:5677
#define FLT_NEQ(x, y)
Definition: rt_api.h:2158
double MinY
Definition: rt_api.h:154
double UpperLeftX
Definition: rt_api.h:157
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
GDALWarpOptions * wopts
Definition: rt_api.c:9583
rt_raster rt_raster_compute_skewed_raster(rt_envelope extent, double *skew, double *scale, double tolerance)
Definition: rt_api.c:6764
double MaxY
Definition: rt_api.h:155
double MaxX
Definition: rt_api.h:153
tuple gt
Definition: window.py:79
Definition: rt_api.h:184
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_api.c:1900
char * srs
Definition: rt_api.c:9579
tuple band
Definition: ovdump.py:57
tuple rast
Definition: rtpixdump.py:62
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_api.c:8989
GDALDatasetH ds
Definition: rt_api.c:9578
double ipY
Definition: rt_api.h:2221
static _rti_warp_arg _rti_warp_arg_init()
Definition: rt_api.c:9603
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_api.c:5661
uint16_t height
Definition: rt_api.h:2227
static void _rti_warp_arg_destroy(_rti_warp_arg arg)
Definition: rt_api.c:9637
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_api.c:5473
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_api.c:6054
void rtinfo(const char *fmt,...)
Definition: rt_api.c:907
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_api.c:6005
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
int rt_util_gdal_driver_registered(const char *drv)
Definition: rt_api.c:461
struct _rti_warp_arg_t::@6 src
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_api.c:5442
void rtwarn(const char *fmt,...)
Definition: rt_api.c:920
rt_pixtype
Definition: rt_api.h:172
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
uint16_t width
Definition: rt_api.h:2226
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:154
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_api.c:3058
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster's geotransform using 6-element array.
Definition: rt_api.c:6026
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
struct _rti_warp_arg_t::@6 dst
void * rtalloc(size_t size)
Raster core memory management functions.
Definition: rt_api.c:867
void rterror(const char *fmt,...)
Raster core error and info handlers.
Definition: rt_api.c:895
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
#define FALSE
Definition: dbfopen.c:169
double MinX
Definition: rt_api.h:152
GDALDriverH drv
Definition: rt_api.c:9577
char * rt_util_gdal_convert_sr(const char *srs, int proj4)
Definition: rt_api.c:323
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_api.c:6105
void * transform
Definition: rt_api.c:9592
double ipX
Definition: rt_api.h:2220
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_api.c:5353
double scaleX
Definition: rt_api.h:2218
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_api.c:229
double UpperLeftY
Definition: rt_api.h:158
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_api.c:2002
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_api.c:5504
if(!(yy_init))
Definition: lwin_wkt_lex.c:860
double scaleY
Definition: rt_api.h:2219

Here is the call graph for this function:

Here is the caller graph for this function: