PostGIS  2.4.9dev-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 3444 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().

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