PostGIS  2.5.0beta1dev-r@@SVN_REVISION@@

◆ rt_raster_from_two_rasters()

rt_errorstate rt_raster_from_two_rasters ( rt_raster  rast1,
rt_raster  rast2,
rt_extenttype  extenttype,
rt_raster rtnraster,
double *  offset 
)

Definition at line 3465 of file rt_raster.c.

References ES_ERROR, ES_NONE, ET_CUSTOM, ET_FIRST, ET_INTERSECTION, ET_LAST, ET_SECOND, ET_UNION, window::gt, rt_raster_t::height, rtrowdump::raster, RASTER_DEBUG, RASTER_DEBUGF, rt_raster_cell_to_geopoint(), rt_raster_destroy(), rt_raster_geopoint_to_cell(), rt_raster_get_geotransform_matrix(), rt_raster_new(), rt_raster_same_alignment(), rt_raster_set_geotransform_matrix(), rt_raster_set_scale(), rt_raster_set_srid(), rterror(), and rt_raster_t::width.

Referenced by RASTER_clip(), RASTER_mapAlgebra2(), RASTER_union_transfn(), rt_raster_iterator(), and test_raster_from_two_rasters().

3469  {
3470  int i;
3471 
3472  rt_raster _rast[2] = {rast1, rast2};
3473  double _offset[2][4] = {{0.}};
3474  uint16_t _dim[2][2] = {{0}};
3475 
3476  rt_raster raster = NULL;
3477  int aligned = 0;
3478  int dim[2] = {0};
3479  double gt[6] = {0.};
3480 
3481  assert(NULL != rast1);
3482  assert(NULL != rast2);
3483  assert(NULL != rtnraster);
3484 
3485  /* set rtnraster to NULL */
3486  *rtnraster = NULL;
3487 
3488  /* rasters must be aligned */
3489  if (rt_raster_same_alignment(rast1, rast2, &aligned, NULL) != ES_NONE) {
3490  rterror("rt_raster_from_two_rasters: Could not test for alignment on the two rasters");
3491  return ES_ERROR;
3492  }
3493  if (!aligned) {
3494  rterror("rt_raster_from_two_rasters: The two rasters provided do not have the same alignment");
3495  return ES_ERROR;
3496  }
3497 
3498  /* dimensions */
3499  _dim[0][0] = rast1->width;
3500  _dim[0][1] = rast1->height;
3501  _dim[1][0] = rast2->width;
3502  _dim[1][1] = rast2->height;
3503 
3504  /* get raster offsets */
3506  _rast[1],
3507  _rast[0]->ipX, _rast[0]->ipY,
3508  &(_offset[1][0]), &(_offset[1][1]),
3509  NULL
3510  ) != ES_NONE) {
3511  rterror("rt_raster_from_two_rasters: Could not compute offsets of the second raster relative to the first raster");
3512  return ES_ERROR;
3513  }
3514  _offset[1][0] = -1 * _offset[1][0];
3515  _offset[1][1] = -1 * _offset[1][1];
3516  _offset[1][2] = _offset[1][0] + _dim[1][0] - 1;
3517  _offset[1][3] = _offset[1][1] + _dim[1][1] - 1;
3518 
3519  i = -1;
3520  switch (extenttype) {
3521  case ET_FIRST:
3522  i = 0;
3523  _offset[0][0] = 0.;
3524  _offset[0][1] = 0.;
3525  /* FALLTHROUGH */
3526  case ET_LAST:
3527  case ET_SECOND:
3528  if (i < 0) {
3529  i = 1;
3530  _offset[0][0] = -1 * _offset[1][0];
3531  _offset[0][1] = -1 * _offset[1][1];
3532  _offset[1][0] = 0.;
3533  _offset[1][1] = 0.;
3534  }
3535 
3536  dim[0] = _dim[i][0];
3537  dim[1] = _dim[i][1];
3538  raster = rt_raster_new(
3539  dim[0],
3540  dim[1]
3541  );
3542  if (raster == NULL) {
3543  rterror("rt_raster_from_two_rasters: Could not create output raster");
3544  return ES_ERROR;
3545  }
3546  rt_raster_set_srid(raster, _rast[i]->srid);
3547  rt_raster_get_geotransform_matrix(_rast[i], gt);
3549  break;
3550  case ET_UNION: {
3551  double off[4] = {0};
3552 
3553  rt_raster_get_geotransform_matrix(_rast[0], gt);
3554  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3555  gt[0],
3556  gt[1],
3557  gt[2],
3558  gt[3],
3559  gt[4],
3560  gt[5]
3561  );
3562 
3563  /* new raster upper left offset */
3564  off[0] = 0;
3565  if (_offset[1][0] < 0)
3566  off[0] = _offset[1][0];
3567  off[1] = 0;
3568  if (_offset[1][1] < 0)
3569  off[1] = _offset[1][1];
3570 
3571  /* new raster lower right offset */
3572  off[2] = _dim[0][0] - 1;
3573  if ((int) _offset[1][2] >= _dim[0][0])
3574  off[2] = _offset[1][2];
3575  off[3] = _dim[0][1] - 1;
3576  if ((int) _offset[1][3] >= _dim[0][1])
3577  off[3] = _offset[1][3];
3578 
3579  /* upper left corner */
3581  _rast[0],
3582  off[0], off[1],
3583  &(gt[0]), &(gt[3]),
3584  NULL
3585  ) != ES_NONE) {
3586  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3587  return ES_ERROR;
3588  }
3589 
3590  dim[0] = off[2] - off[0] + 1;
3591  dim[1] = off[3] - off[1] + 1;
3592  RASTER_DEBUGF(4, "off = (%f, %f, %f, %f)",
3593  off[0],
3594  off[1],
3595  off[2],
3596  off[3]
3597  );
3598  RASTER_DEBUGF(4, "dim = (%d, %d)", dim[0], dim[1]);
3599 
3600  raster = rt_raster_new(
3601  dim[0],
3602  dim[1]
3603  );
3604  if (raster == NULL) {
3605  rterror("rt_raster_from_two_rasters: Could not create output raster");
3606  return ES_ERROR;
3607  }
3608  rt_raster_set_srid(raster, _rast[0]->srid);
3610  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3611  gt[0],
3612  gt[1],
3613  gt[2],
3614  gt[3],
3615  gt[4],
3616  gt[5]
3617  );
3618 
3619  /* get offsets */
3621  _rast[0],
3622  gt[0], gt[3],
3623  &(_offset[0][0]), &(_offset[0][1]),
3624  NULL
3625  ) != ES_NONE) {
3626  rterror("rt_raster_from_two_rasters: Could not get offsets of the FIRST raster relative to the output raster");
3627  rt_raster_destroy(raster);
3628  return ES_ERROR;
3629  }
3630  _offset[0][0] *= -1;
3631  _offset[0][1] *= -1;
3632 
3634  _rast[1],
3635  gt[0], gt[3],
3636  &(_offset[1][0]), &(_offset[1][1]),
3637  NULL
3638  ) != ES_NONE) {
3639  rterror("rt_raster_from_two_rasters: Could not get offsets of the SECOND raster relative to the output raster");
3640  rt_raster_destroy(raster);
3641  return ES_ERROR;
3642  }
3643  _offset[1][0] *= -1;
3644  _offset[1][1] *= -1;
3645  break;
3646  }
3647  case ET_INTERSECTION: {
3648  double off[4] = {0};
3649 
3650  /* no intersection */
3651  if (
3652  (_offset[1][2] < 0 || _offset[1][0] > (_dim[0][0] - 1)) ||
3653  (_offset[1][3] < 0 || _offset[1][1] > (_dim[0][1] - 1))
3654  ) {
3655  RASTER_DEBUG(3, "The two rasters provided have no intersection. Returning no band raster");
3656 
3657  raster = rt_raster_new(0, 0);
3658  if (raster == NULL) {
3659  rterror("rt_raster_from_two_rasters: Could not create output raster");
3660  return ES_ERROR;
3661  }
3662  rt_raster_set_srid(raster, _rast[0]->srid);
3663  rt_raster_set_scale(raster, 0, 0);
3664 
3665  /* set offsets if provided */
3666  if (NULL != offset) {
3667  for (i = 0; i < 4; i++)
3668  offset[i] = _offset[i / 2][i % 2];
3669  }
3670 
3671  *rtnraster = raster;
3672  return ES_NONE;
3673  }
3674 
3675  if (_offset[1][0] > 0)
3676  off[0] = _offset[1][0];
3677  if (_offset[1][1] > 0)
3678  off[1] = _offset[1][1];
3679 
3680  off[2] = _dim[0][0] - 1;
3681  if (_offset[1][2] < _dim[0][0])
3682  off[2] = _offset[1][2];
3683  off[3] = _dim[0][1] - 1;
3684  if (_offset[1][3] < _dim[0][1])
3685  off[3] = _offset[1][3];
3686 
3687  dim[0] = off[2] - off[0] + 1;
3688  dim[1] = off[3] - off[1] + 1;
3689  raster = rt_raster_new(
3690  dim[0],
3691  dim[1]
3692  );
3693  if (raster == NULL) {
3694  rterror("rt_raster_from_two_rasters: Could not create output raster");
3695  return ES_ERROR;
3696  }
3697  rt_raster_set_srid(raster, _rast[0]->srid);
3698 
3699  /* get upper-left corner */
3700  rt_raster_get_geotransform_matrix(_rast[0], gt);
3702  _rast[0],
3703  off[0], off[1],
3704  &(gt[0]), &(gt[3]),
3705  gt
3706  ) != ES_NONE) {
3707  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3708  rt_raster_destroy(raster);
3709  return ES_ERROR;
3710  }
3711 
3713 
3714  /* get offsets */
3716  _rast[0],
3717  gt[0], gt[3],
3718  &(_offset[0][0]), &(_offset[0][1]),
3719  NULL
3720  ) != ES_NONE) {
3721  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the FIRST raster relative to the output raster");
3722  rt_raster_destroy(raster);
3723  return ES_ERROR;
3724  }
3725  _offset[0][0] *= -1;
3726  _offset[0][1] *= -1;
3727 
3729  _rast[1],
3730  gt[0], gt[3],
3731  &(_offset[1][0]), &(_offset[1][1]),
3732  NULL
3733  ) != ES_NONE) {
3734  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the SECOND raster relative to the output raster");
3735  rt_raster_destroy(raster);
3736  return ES_ERROR;
3737  }
3738  _offset[1][0] *= -1;
3739  _offset[1][1] *= -1;
3740  break;
3741  }
3742  case ET_CUSTOM:
3743  rterror("rt_raster_from_two_rasters: Extent type ET_CUSTOM is not supported by this function");
3744  break;
3745  }
3746 
3747  /* set offsets if provided */
3748  if (NULL != offset) {
3749  for (i = 0; i < 4; i++)
3750  offset[i] = _offset[i / 2][i % 2];
3751  }
3752 
3753  *rtnraster = raster;
3754  return ES_NONE;
3755 }
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_raster.c:137
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
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster&#39;s SRID.
Definition: rt_raster.c:363
gt
Definition: window.py:77
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
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
uint16_t height
Definition: librtcore.h:2302
uint16_t width
Definition: librtcore.h:2301
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:706
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
Here is the call graph for this function:
Here is the caller graph for this function: