PostGIS  2.5.7dev-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 3361 of file rt_raster.c.

3365  {
3366  int i;
3367 
3368  rt_raster _rast[2] = {rast1, rast2};
3369  double _offset[2][4] = {{0.}};
3370  uint16_t _dim[2][2] = {{0}};
3371 
3372  rt_raster raster = NULL;
3373  int aligned = 0;
3374  int dim[2] = {0};
3375  double gt[6] = {0.};
3376 
3377  assert(NULL != rast1);
3378  assert(NULL != rast2);
3379  assert(NULL != rtnraster);
3380 
3381  /* set rtnraster to NULL */
3382  *rtnraster = NULL;
3383 
3384  /* rasters must be aligned */
3385  if (rt_raster_same_alignment(rast1, rast2, &aligned, NULL) != ES_NONE) {
3386  rterror("rt_raster_from_two_rasters: Could not test for alignment on the two rasters");
3387  return ES_ERROR;
3388  }
3389  if (!aligned) {
3390  rterror("rt_raster_from_two_rasters: The two rasters provided do not have the same alignment");
3391  return ES_ERROR;
3392  }
3393 
3394  /* dimensions */
3395  _dim[0][0] = rast1->width;
3396  _dim[0][1] = rast1->height;
3397  _dim[1][0] = rast2->width;
3398  _dim[1][1] = rast2->height;
3399 
3400  /* get raster offsets */
3401  if (rt_raster_geopoint_to_cell(
3402  _rast[1],
3403  _rast[0]->ipX, _rast[0]->ipY,
3404  &(_offset[1][0]), &(_offset[1][1]),
3405  NULL
3406  ) != ES_NONE) {
3407  rterror("rt_raster_from_two_rasters: Could not compute offsets of the second raster relative to the first raster");
3408  return ES_ERROR;
3409  }
3410  _offset[1][0] = -1 * _offset[1][0];
3411  _offset[1][1] = -1 * _offset[1][1];
3412  _offset[1][2] = _offset[1][0] + _dim[1][0] - 1;
3413  _offset[1][3] = _offset[1][1] + _dim[1][1] - 1;
3414 
3415  i = -1;
3416  switch (extenttype) {
3417  case ET_FIRST:
3418  i = 0;
3419  _offset[0][0] = 0.;
3420  _offset[0][1] = 0.;
3421  /* FALLTHROUGH */
3422  case ET_LAST:
3423  case ET_SECOND:
3424  if (i < 0) {
3425  i = 1;
3426  _offset[0][0] = -1 * _offset[1][0];
3427  _offset[0][1] = -1 * _offset[1][1];
3428  _offset[1][0] = 0.;
3429  _offset[1][1] = 0.;
3430  }
3431 
3432  dim[0] = _dim[i][0];
3433  dim[1] = _dim[i][1];
3434  raster = rt_raster_new(
3435  dim[0],
3436  dim[1]
3437  );
3438  if (raster == NULL) {
3439  rterror("rt_raster_from_two_rasters: Could not create output raster");
3440  return ES_ERROR;
3441  }
3442  rt_raster_set_srid(raster, _rast[i]->srid);
3443  rt_raster_get_geotransform_matrix(_rast[i], gt);
3444  rt_raster_set_geotransform_matrix(raster, gt);
3445  break;
3446  case ET_UNION: {
3447  double off[4] = {0};
3448 
3449  rt_raster_get_geotransform_matrix(_rast[0], gt);
3450  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3451  gt[0],
3452  gt[1],
3453  gt[2],
3454  gt[3],
3455  gt[4],
3456  gt[5]
3457  );
3458 
3459  /* new raster upper left offset */
3460  off[0] = 0;
3461  if (_offset[1][0] < 0)
3462  off[0] = _offset[1][0];
3463  off[1] = 0;
3464  if (_offset[1][1] < 0)
3465  off[1] = _offset[1][1];
3466 
3467  /* new raster lower right offset */
3468  off[2] = _dim[0][0] - 1;
3469  if ((int) _offset[1][2] >= _dim[0][0])
3470  off[2] = _offset[1][2];
3471  off[3] = _dim[0][1] - 1;
3472  if ((int) _offset[1][3] >= _dim[0][1])
3473  off[3] = _offset[1][3];
3474 
3475  /* upper left corner */
3476  if (rt_raster_cell_to_geopoint(
3477  _rast[0],
3478  off[0], off[1],
3479  &(gt[0]), &(gt[3]),
3480  NULL
3481  ) != ES_NONE) {
3482  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3483  return ES_ERROR;
3484  }
3485 
3486  dim[0] = off[2] - off[0] + 1;
3487  dim[1] = off[3] - off[1] + 1;
3488  RASTER_DEBUGF(4, "off = (%f, %f, %f, %f)",
3489  off[0],
3490  off[1],
3491  off[2],
3492  off[3]
3493  );
3494  RASTER_DEBUGF(4, "dim = (%d, %d)", dim[0], dim[1]);
3495 
3496  raster = rt_raster_new(
3497  dim[0],
3498  dim[1]
3499  );
3500  if (raster == NULL) {
3501  rterror("rt_raster_from_two_rasters: Could not create output raster");
3502  return ES_ERROR;
3503  }
3504  rt_raster_set_srid(raster, _rast[0]->srid);
3505  rt_raster_set_geotransform_matrix(raster, gt);
3506  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
3507  gt[0],
3508  gt[1],
3509  gt[2],
3510  gt[3],
3511  gt[4],
3512  gt[5]
3513  );
3514 
3515  /* get offsets */
3516  if (rt_raster_geopoint_to_cell(
3517  _rast[0],
3518  gt[0], gt[3],
3519  &(_offset[0][0]), &(_offset[0][1]),
3520  NULL
3521  ) != ES_NONE) {
3522  rterror("rt_raster_from_two_rasters: Could not get offsets of the FIRST raster relative to the output raster");
3523  rt_raster_destroy(raster);
3524  return ES_ERROR;
3525  }
3526  _offset[0][0] *= -1;
3527  _offset[0][1] *= -1;
3528 
3529  if (rt_raster_geopoint_to_cell(
3530  _rast[1],
3531  gt[0], gt[3],
3532  &(_offset[1][0]), &(_offset[1][1]),
3533  NULL
3534  ) != ES_NONE) {
3535  rterror("rt_raster_from_two_rasters: Could not get offsets of the SECOND raster relative to the output raster");
3536  rt_raster_destroy(raster);
3537  return ES_ERROR;
3538  }
3539  _offset[1][0] *= -1;
3540  _offset[1][1] *= -1;
3541  break;
3542  }
3543  case ET_INTERSECTION: {
3544  double off[4] = {0};
3545 
3546  /* no intersection */
3547  if (
3548  (_offset[1][2] < 0 || _offset[1][0] > (_dim[0][0] - 1)) ||
3549  (_offset[1][3] < 0 || _offset[1][1] > (_dim[0][1] - 1))
3550  ) {
3551  RASTER_DEBUG(3, "The two rasters provided have no intersection. Returning no band raster");
3552 
3553  raster = rt_raster_new(0, 0);
3554  if (raster == NULL) {
3555  rterror("rt_raster_from_two_rasters: Could not create output raster");
3556  return ES_ERROR;
3557  }
3558  rt_raster_set_srid(raster, _rast[0]->srid);
3559  rt_raster_set_scale(raster, 0, 0);
3560 
3561  /* set offsets if provided */
3562  if (NULL != offset) {
3563  for (i = 0; i < 4; i++)
3564  offset[i] = _offset[i / 2][i % 2];
3565  }
3566 
3567  *rtnraster = raster;
3568  return ES_NONE;
3569  }
3570 
3571  if (_offset[1][0] > 0)
3572  off[0] = _offset[1][0];
3573  if (_offset[1][1] > 0)
3574  off[1] = _offset[1][1];
3575 
3576  off[2] = _dim[0][0] - 1;
3577  if (_offset[1][2] < _dim[0][0])
3578  off[2] = _offset[1][2];
3579  off[3] = _dim[0][1] - 1;
3580  if (_offset[1][3] < _dim[0][1])
3581  off[3] = _offset[1][3];
3582 
3583  dim[0] = off[2] - off[0] + 1;
3584  dim[1] = off[3] - off[1] + 1;
3585  raster = rt_raster_new(
3586  dim[0],
3587  dim[1]
3588  );
3589  if (raster == NULL) {
3590  rterror("rt_raster_from_two_rasters: Could not create output raster");
3591  return ES_ERROR;
3592  }
3593  rt_raster_set_srid(raster, _rast[0]->srid);
3594 
3595  /* get upper-left corner */
3596  rt_raster_get_geotransform_matrix(_rast[0], gt);
3597  if (rt_raster_cell_to_geopoint(
3598  _rast[0],
3599  off[0], off[1],
3600  &(gt[0]), &(gt[3]),
3601  gt
3602  ) != ES_NONE) {
3603  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
3604  rt_raster_destroy(raster);
3605  return ES_ERROR;
3606  }
3607 
3608  rt_raster_set_geotransform_matrix(raster, gt);
3609 
3610  /* get offsets */
3611  if (rt_raster_geopoint_to_cell(
3612  _rast[0],
3613  gt[0], gt[3],
3614  &(_offset[0][0]), &(_offset[0][1]),
3615  NULL
3616  ) != ES_NONE) {
3617  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the FIRST raster relative to the output raster");
3618  rt_raster_destroy(raster);
3619  return ES_ERROR;
3620  }
3621  _offset[0][0] *= -1;
3622  _offset[0][1] *= -1;
3623 
3624  if (rt_raster_geopoint_to_cell(
3625  _rast[1],
3626  gt[0], gt[3],
3627  &(_offset[1][0]), &(_offset[1][1]),
3628  NULL
3629  ) != ES_NONE) {
3630  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the SECOND raster relative to the output raster");
3631  rt_raster_destroy(raster);
3632  return ES_ERROR;
3633  }
3634  _offset[1][0] *= -1;
3635  _offset[1][1] *= -1;
3636  break;
3637  }
3638  case ET_CUSTOM:
3639  rterror("rt_raster_from_two_rasters: Extent type ET_CUSTOM is not supported by this function");
3640  break;
3641  }
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 }
rterror
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
ES_NONE
@ ES_NONE
Definition: librtcore.h:180
ES_ERROR
@ ES_ERROR
Definition: librtcore.h:181
ET_CUSTOM
@ ET_CUSTOM
Definition: librtcore.h:206
ET_LAST
@ ET_LAST
Definition: librtcore.h:205
ET_INTERSECTION
@ ET_INTERSECTION
Definition: librtcore.h:201
ET_UNION
@ ET_UNION
Definition: librtcore.h:202
ET_SECOND
@ ET_SECOND
Definition: librtcore.h:204
ET_FIRST
@ ET_FIRST
Definition: librtcore.h:203
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:77
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
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:727
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_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_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
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
rt_raster_set_srid
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_raster.c:363
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
rt_raster_t::width
uint16_t width
Definition: librtcore.h:2301
rt_raster_t::height
uint16_t height
Definition: librtcore.h:2302
rt_raster_t
Definition: librtcore.h:2284

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: