PostGIS  3.1.6dev-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 3337 of file rt_raster.c.

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

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: