PostGIS  3.0.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 3350 of file rt_raster.c.

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