PostGIS  3.4.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 3464 of file rt_raster.c.

3468  {
3469  int i;
3470 
3471  rt_raster _rast[2] = {rast1, rast2};
3472  double _offset[2][4] = {{0.}};
3473  uint16_t _dim[2][2] = {{0}};
3474 
3475  rt_raster raster = NULL;
3476  int aligned = 0;
3477  int dim[2] = {0};
3478  double gt[6] = {0.};
3479 
3480  assert(NULL != rast1);
3481  assert(NULL != rast2);
3482  assert(NULL != rtnraster);
3483 
3484  /* set rtnraster to NULL */
3485  *rtnraster = NULL;
3486 
3487  /* rasters must be aligned */
3488  if (rt_raster_same_alignment(rast1, rast2, &aligned, NULL) != ES_NONE) {
3489  rterror("rt_raster_from_two_rasters: Could not test for alignment on the two rasters");
3490  return ES_ERROR;
3491  }
3492  if (!aligned) {
3493  rterror("rt_raster_from_two_rasters: The two rasters provided do not have the same alignment");
3494  return ES_ERROR;
3495  }
3496 
3497  /* dimensions */
3498  _dim[0][0] = rast1->width;
3499  _dim[0][1] = rast1->height;
3500  _dim[1][0] = rast2->width;
3501  _dim[1][1] = rast2->height;
3502 
3503  /* get raster offsets */
3504  if (rt_raster_geopoint_to_cell(
3505  _rast[1],
3506  _rast[0]->ipX, _rast[0]->ipY,
3507  &(_offset[1][0]), &(_offset[1][1]),
3508  NULL
3509  ) != ES_NONE) {
3510  rterror("rt_raster_from_two_rasters: Could not compute offsets of the second raster relative to the first raster");
3511  return ES_ERROR;
3512  }
3513  _offset[1][0] = -1 * _offset[1][0];
3514  _offset[1][1] = -1 * _offset[1][1];
3515  _offset[1][2] = _offset[1][0] + _dim[1][0] - 1;
3516  _offset[1][3] = _offset[1][1] + _dim[1][1] - 1;
3517 
3518  i = -1;
3519  switch (extenttype) {
3520  case ET_FIRST:
3521  i = 0;
3522  _offset[0][0] = 0.;
3523  _offset[0][1] = 0.;
3524  /* FALLTHROUGH */
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);
3547  rt_raster_set_geotransform_matrix(raster, 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 */
3579  if (rt_raster_cell_to_geopoint(
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);
3608  rt_raster_set_geotransform_matrix(raster, gt);
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 */
3619  if (rt_raster_geopoint_to_cell(
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 
3632  if (rt_raster_geopoint_to_cell(
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);
3700  if (rt_raster_cell_to_geopoint(
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 
3711  rt_raster_set_geotransform_matrix(raster, gt);
3712 
3713  /* get offsets */
3714  if (rt_raster_geopoint_to_cell(
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 
3727  if (rt_raster_geopoint_to_cell(
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 }
rterror
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:219
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:302
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:306
ES_NONE
@ ES_NONE
Definition: librtcore.h:182
ES_ERROR
@ ES_ERROR
Definition: librtcore.h:183
ET_CUSTOM
@ ET_CUSTOM
Definition: librtcore.h:208
ET_LAST
@ ET_LAST
Definition: librtcore.h:207
ET_INTERSECTION
@ ET_INTERSECTION
Definition: librtcore.h:203
ET_UNION
@ ET_UNION
Definition: librtcore.h:204
ET_SECOND
@ ET_SECOND
Definition: librtcore.h:206
ET_FIRST
@ ET_FIRST
Definition: librtcore.h:205
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's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
window.gt
gt
Definition: window.py:78
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:759
rt_raster_set_geotransform_matrix
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster's geotransform using 6-element array.
Definition: rt_raster.c:731
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:141
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 cell coordinate.
Definition: rt_raster.c:808
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:86
rt_raster_new
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:52
rt_raster_set_srid
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_raster.c:367
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:710
rt_raster_t::width
uint16_t width
Definition: librtcore.h:2454
rt_raster_t::height
uint16_t height
Definition: librtcore.h:2455
rt_raster_t
Definition: librtcore.h:2437

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().

Here is the call graph for this function:
Here is the caller graph for this function: