PostGIS 3.7.0dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ 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:304
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:308
ES_NONE
@ ES_NONE
Definition librtcore.h:182
ES_ERROR
@ ES_ERROR
Definition librtcore.h:183
ET_CUSTOM
@ ET_CUSTOM
Definition librtcore.h:210
ET_LAST
@ ET_LAST
Definition librtcore.h:209
ET_INTERSECTION
@ ET_INTERSECTION
Definition librtcore.h:205
ET_UNION
@ ET_UNION
Definition librtcore.h:206
ET_SECOND
@ ET_SECOND
Definition librtcore.h:208
ET_FIRST
@ ET_FIRST
Definition librtcore.h:207
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:125
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:2503
rt_raster_t::height
uint16_t height
Definition librtcore.h:2504
rt_raster_t
Definition librtcore.h:2486

References ES_ERROR, ES_NONE, ET_CUSTOM, ET_FIRST, ET_INTERSECTION, ET_LAST, ET_SECOND, ET_UNION, rt_raster_t::height, 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: