PostGIS  3.7.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 3348 of file rt_raster.c.

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

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: