PostGIS  2.5.0dev-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  case ET_LAST:
3526  case ET_SECOND:
3527  if (i < 0) {
3528  i = 1;
3529  _offset[0][0] = -1 * _offset[1][0];
3530  _offset[0][1] = -1 * _offset[1][1];
3531  _offset[1][0] = 0.;
3532  _offset[1][1] = 0.;
3533  }
3534 
3535  dim[0] = _dim[i][0];
3536  dim[1] = _dim[i][1];
3537  raster = rt_raster_new(
3538  dim[0],
3539  dim[1]
3540  );
3541  if (raster == NULL) {
3542  rterror("rt_raster_from_two_rasters: Could not create output raster");
3543  return ES_ERROR;
3544  }
3545  rt_raster_set_srid(raster, _rast[i]->srid);
3546  rt_raster_get_geotransform_matrix(_rast[i], gt);
3548  break;
3549  case ET_UNION: {
3550  double off[4] = {0};
3551 
3552  rt_raster_get_geotransform_matrix(_rast[0], gt);
3553  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3554  gt[0],
3555  gt[1],
3556  gt[2],
3557  gt[3],
3558  gt[4],
3559  gt[5]
3560  );
3561 
3562  /* new raster upper left offset */
3563  off[0] = 0;
3564  if (_offset[1][0] < 0)
3565  off[0] = _offset[1][0];
3566  off[1] = 0;
3567  if (_offset[1][1] < 0)
3568  off[1] = _offset[1][1];
3569 
3570  /* new raster lower right offset */
3571  off[2] = _dim[0][0] - 1;
3572  if ((int) _offset[1][2] >= _dim[0][0])
3573  off[2] = _offset[1][2];
3574  off[3] = _dim[0][1] - 1;
3575  if ((int) _offset[1][3] >= _dim[0][1])
3576  off[3] = _offset[1][3];
3577 
3578  /* upper left corner */
3580  _rast[0],
3581  off[0], off[1],
3582  &(gt[0]), &(gt[3]),
3583  NULL
3584  ) != ES_NONE) {
3585  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3586  return ES_ERROR;
3587  }
3588 
3589  dim[0] = off[2] - off[0] + 1;
3590  dim[1] = off[3] - off[1] + 1;
3591  RASTER_DEBUGF(4, "off = (%f, %f, %f, %f)",
3592  off[0],
3593  off[1],
3594  off[2],
3595  off[3]
3596  );
3597  RASTER_DEBUGF(4, "dim = (%d, %d)", dim[0], dim[1]);
3598 
3599  raster = rt_raster_new(
3600  dim[0],
3601  dim[1]
3602  );
3603  if (raster == NULL) {
3604  rterror("rt_raster_from_two_rasters: Could not create output raster");
3605  return ES_ERROR;
3606  }
3607  rt_raster_set_srid(raster, _rast[0]->srid);
3609  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3610  gt[0],
3611  gt[1],
3612  gt[2],
3613  gt[3],
3614  gt[4],
3615  gt[5]
3616  );
3617 
3618  /* get offsets */
3620  _rast[0],
3621  gt[0], gt[3],
3622  &(_offset[0][0]), &(_offset[0][1]),
3623  NULL
3624  ) != ES_NONE) {
3625  rterror("rt_raster_from_two_rasters: Could not get offsets of the FIRST raster relative to the output raster");
3626  rt_raster_destroy(raster);
3627  return ES_ERROR;
3628  }
3629  _offset[0][0] *= -1;
3630  _offset[0][1] *= -1;
3631 
3633  _rast[1],
3634  gt[0], gt[3],
3635  &(_offset[1][0]), &(_offset[1][1]),
3636  NULL
3637  ) != ES_NONE) {
3638  rterror("rt_raster_from_two_rasters: Could not get offsets of the SECOND raster relative to the output raster");
3639  rt_raster_destroy(raster);
3640  return ES_ERROR;
3641  }
3642  _offset[1][0] *= -1;
3643  _offset[1][1] *= -1;
3644  break;
3645  }
3646  case ET_INTERSECTION: {
3647  double off[4] = {0};
3648 
3649  /* no intersection */
3650  if (
3651  (_offset[1][2] < 0 || _offset[1][0] > (_dim[0][0] - 1)) ||
3652  (_offset[1][3] < 0 || _offset[1][1] > (_dim[0][1] - 1))
3653  ) {
3654  RASTER_DEBUG(3, "The two rasters provided have no intersection. Returning no band raster");
3655 
3656  raster = rt_raster_new(0, 0);
3657  if (raster == NULL) {
3658  rterror("rt_raster_from_two_rasters: Could not create output raster");
3659  return ES_ERROR;
3660  }
3661  rt_raster_set_srid(raster, _rast[0]->srid);
3662  rt_raster_set_scale(raster, 0, 0);
3663 
3664  /* set offsets if provided */
3665  if (NULL != offset) {
3666  for (i = 0; i < 4; i++)
3667  offset[i] = _offset[i / 2][i % 2];
3668  }
3669 
3670  *rtnraster = raster;
3671  return ES_NONE;
3672  }
3673 
3674  if (_offset[1][0] > 0)
3675  off[0] = _offset[1][0];
3676  if (_offset[1][1] > 0)
3677  off[1] = _offset[1][1];
3678 
3679  off[2] = _dim[0][0] - 1;
3680  if (_offset[1][2] < _dim[0][0])
3681  off[2] = _offset[1][2];
3682  off[3] = _dim[0][1] - 1;
3683  if (_offset[1][3] < _dim[0][1])
3684  off[3] = _offset[1][3];
3685 
3686  dim[0] = off[2] - off[0] + 1;
3687  dim[1] = off[3] - off[1] + 1;
3688  raster = rt_raster_new(
3689  dim[0],
3690  dim[1]
3691  );
3692  if (raster == NULL) {
3693  rterror("rt_raster_from_two_rasters: Could not create output raster");
3694  return ES_ERROR;
3695  }
3696  rt_raster_set_srid(raster, _rast[0]->srid);
3697 
3698  /* get upper-left corner */
3699  rt_raster_get_geotransform_matrix(_rast[0], gt);
3701  _rast[0],
3702  off[0], off[1],
3703  &(gt[0]), &(gt[3]),
3704  gt
3705  ) != ES_NONE) {
3706  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3707  rt_raster_destroy(raster);
3708  return ES_ERROR;
3709  }
3710 
3712 
3713  /* get offsets */
3715  _rast[0],
3716  gt[0], gt[3],
3717  &(_offset[0][0]), &(_offset[0][1]),
3718  NULL
3719  ) != ES_NONE) {
3720  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the FIRST raster relative to the output raster");
3721  rt_raster_destroy(raster);
3722  return ES_ERROR;
3723  }
3724  _offset[0][0] *= -1;
3725  _offset[0][1] *= -1;
3726 
3728  _rast[1],
3729  gt[0], gt[3],
3730  &(_offset[1][0]), &(_offset[1][1]),
3731  NULL
3732  ) != ES_NONE) {
3733  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the SECOND raster relative to the output raster");
3734  rt_raster_destroy(raster);
3735  return ES_ERROR;
3736  }
3737  _offset[1][0] *= -1;
3738  _offset[1][1] *= -1;
3739  break;
3740  }
3741  case ET_CUSTOM:
3742  rterror("rt_raster_from_two_rasters: Extent type ET_CUSTOM is not supported by this function");
3743  break;
3744  }
3745 
3746  /* set offsets if provided */
3747  if (NULL != offset) {
3748  for (i = 0; i < 4; i++)
3749  offset[i] = _offset[i / 2][i % 2];
3750  }
3751 
3752  *rtnraster = raster;
3753  return ES_NONE;
3754 }
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:2282
uint16_t width
Definition: librtcore.h:2281
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: