PostGIS 3.0.6dev-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 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
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, 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: