PostGIS  2.5.0beta2dev-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 3461 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().

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