PostGIS  2.1.10dev-r@@SVN_REVISION@@
rt_errorstate rt_raster_from_two_rasters ( rt_raster  rast1,
rt_raster  rast2,
rt_extenttype  extenttype,
rt_raster rtnraster,
double *  offset 
)

Definition at line 12876 of file rt_api.c.

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().

12880  {
12881  int i;
12882 
12883  rt_raster _rast[2] = {rast1, rast2};
12884  double _offset[2][4] = {{0.}};
12885  uint16_t _dim[2][2] = {{0}};
12886 
12887  rt_raster raster = NULL;
12888  int aligned = 0;
12889  int dim[2] = {0};
12890  double gt[6] = {0.};
12891 
12892  assert(NULL != rast1);
12893  assert(NULL != rast2);
12894  assert(NULL != rtnraster);
12895 
12896  /* set rtnraster to NULL */
12897  *rtnraster = NULL;
12898 
12899  /* rasters must be aligned */
12900  if (rt_raster_same_alignment(rast1, rast2, &aligned, NULL) != ES_NONE) {
12901  rterror("rt_raster_from_two_rasters: Could not test for alignment on the two rasters");
12902  return ES_ERROR;
12903  }
12904  if (!aligned) {
12905  rterror("rt_raster_from_two_rasters: The two rasters provided do not have the same alignment");
12906  return ES_ERROR;
12907  }
12908 
12909  /* dimensions */
12910  _dim[0][0] = rast1->width;
12911  _dim[0][1] = rast1->height;
12912  _dim[1][0] = rast2->width;
12913  _dim[1][1] = rast2->height;
12914 
12915  /* get raster offsets */
12917  _rast[1],
12918  _rast[0]->ipX, _rast[0]->ipY,
12919  &(_offset[1][0]), &(_offset[1][1]),
12920  NULL
12921  ) != ES_NONE) {
12922  rterror("rt_raster_from_two_rasters: Could not compute offsets of the second raster relative to the first raster");
12923  return ES_ERROR;
12924  }
12925  _offset[1][0] = -1 * _offset[1][0];
12926  _offset[1][1] = -1 * _offset[1][1];
12927  _offset[1][2] = _offset[1][0] + _dim[1][0] - 1;
12928  _offset[1][3] = _offset[1][1] + _dim[1][1] - 1;
12929 
12930  i = -1;
12931  switch (extenttype) {
12932  case ET_FIRST:
12933  i = 0;
12934  _offset[0][0] = 0.;
12935  _offset[0][1] = 0.;
12936  case ET_LAST:
12937  case ET_SECOND:
12938  if (i < 0) {
12939  i = 1;
12940  _offset[0][0] = -1 * _offset[1][0];
12941  _offset[0][1] = -1 * _offset[1][1];
12942  _offset[1][0] = 0.;
12943  _offset[1][1] = 0.;
12944  }
12945 
12946  dim[0] = _dim[i][0];
12947  dim[1] = _dim[i][1];
12948  raster = rt_raster_new(
12949  dim[0],
12950  dim[1]
12951  );
12952  if (raster == NULL) {
12953  rterror("rt_raster_from_two_rasters: Could not create output raster");
12954  return ES_ERROR;
12955  }
12956  rt_raster_set_srid(raster, _rast[i]->srid);
12957  rt_raster_get_geotransform_matrix(_rast[i], gt);
12959  break;
12960  case ET_UNION: {
12961  double off[4] = {0};
12962 
12963  rt_raster_get_geotransform_matrix(_rast[0], gt);
12964  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
12965  gt[0],
12966  gt[1],
12967  gt[2],
12968  gt[3],
12969  gt[4],
12970  gt[5]
12971  );
12972 
12973  /* new raster upper left offset */
12974  off[0] = 0;
12975  if (_offset[1][0] < 0)
12976  off[0] = _offset[1][0];
12977  off[1] = 0;
12978  if (_offset[1][1] < 0)
12979  off[1] = _offset[1][1];
12980 
12981  /* new raster lower right offset */
12982  off[2] = _dim[0][0] - 1;
12983  if ((int) _offset[1][2] >= _dim[0][0])
12984  off[2] = _offset[1][2];
12985  off[3] = _dim[0][1] - 1;
12986  if ((int) _offset[1][3] >= _dim[0][1])
12987  off[3] = _offset[1][3];
12988 
12989  /* upper left corner */
12991  _rast[0],
12992  off[0], off[1],
12993  &(gt[0]), &(gt[3]),
12994  NULL
12995  ) != ES_NONE) {
12996  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
12997  return ES_ERROR;
12998  }
12999 
13000  dim[0] = off[2] - off[0] + 1;
13001  dim[1] = off[3] - off[1] + 1;
13002  RASTER_DEBUGF(4, "off = (%f, %f, %f, %f)",
13003  off[0],
13004  off[1],
13005  off[2],
13006  off[3]
13007  );
13008  RASTER_DEBUGF(4, "dim = (%d, %d)", dim[0], dim[1]);
13009 
13010  raster = rt_raster_new(
13011  dim[0],
13012  dim[1]
13013  );
13014  if (raster == NULL) {
13015  rterror("rt_raster_from_two_rasters: Could not create output raster");
13016  return ES_ERROR;
13017  }
13018  rt_raster_set_srid(raster, _rast[0]->srid);
13020  RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
13021  gt[0],
13022  gt[1],
13023  gt[2],
13024  gt[3],
13025  gt[4],
13026  gt[5]
13027  );
13028 
13029  /* get offsets */
13031  _rast[0],
13032  gt[0], gt[3],
13033  &(_offset[0][0]), &(_offset[0][1]),
13034  NULL
13035  ) != ES_NONE) {
13036  rterror("rt_raster_from_two_rasters: Could not get offsets of the FIRST raster relative to the output raster");
13037  rt_raster_destroy(raster);
13038  return ES_ERROR;
13039  }
13040  _offset[0][0] *= -1;
13041  _offset[0][1] *= -1;
13042 
13044  _rast[1],
13045  gt[0], gt[3],
13046  &(_offset[1][0]), &(_offset[1][1]),
13047  NULL
13048  ) != ES_NONE) {
13049  rterror("rt_raster_from_two_rasters: Could not get offsets of the SECOND raster relative to the output raster");
13050  rt_raster_destroy(raster);
13051  return ES_ERROR;
13052  }
13053  _offset[1][0] *= -1;
13054  _offset[1][1] *= -1;
13055  break;
13056  }
13057  case ET_INTERSECTION: {
13058  double off[4] = {0};
13059 
13060  /* no intersection */
13061  if (
13062  (_offset[1][2] < 0 || _offset[1][0] > (_dim[0][0] - 1)) ||
13063  (_offset[1][3] < 0 || _offset[1][1] > (_dim[0][1] - 1))
13064  ) {
13065  RASTER_DEBUG(3, "The two rasters provided have no intersection. Returning no band raster");
13066 
13067  raster = rt_raster_new(0, 0);
13068  if (raster == NULL) {
13069  rterror("rt_raster_from_two_rasters: Could not create output raster");
13070  return ES_ERROR;
13071  }
13072  rt_raster_set_srid(raster, _rast[0]->srid);
13073  rt_raster_set_scale(raster, 0, 0);
13074 
13075  /* set offsets if provided */
13076  if (NULL != offset) {
13077  for (i = 0; i < 4; i++)
13078  offset[i] = _offset[i / 2][i % 2];
13079  }
13080 
13081  *rtnraster = raster;
13082  return ES_NONE;
13083  }
13084 
13085  if (_offset[1][0] > 0)
13086  off[0] = _offset[1][0];
13087  if (_offset[1][1] > 0)
13088  off[1] = _offset[1][1];
13089 
13090  off[2] = _dim[0][0] - 1;
13091  if (_offset[1][2] < _dim[0][0])
13092  off[2] = _offset[1][2];
13093  off[3] = _dim[0][1] - 1;
13094  if (_offset[1][3] < _dim[0][1])
13095  off[3] = _offset[1][3];
13096 
13097  dim[0] = off[2] - off[0] + 1;
13098  dim[1] = off[3] - off[1] + 1;
13099  raster = rt_raster_new(
13100  dim[0],
13101  dim[1]
13102  );
13103  if (raster == NULL) {
13104  rterror("rt_raster_from_two_rasters: Could not create output raster");
13105  return ES_ERROR;
13106  }
13107  rt_raster_set_srid(raster, _rast[0]->srid);
13108 
13109  /* get upper-left corner */
13110  rt_raster_get_geotransform_matrix(_rast[0], gt);
13112  _rast[0],
13113  off[0], off[1],
13114  &(gt[0]), &(gt[3]),
13115  gt
13116  ) != ES_NONE) {
13117  rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
13118  rt_raster_destroy(raster);
13119  return ES_ERROR;
13120  }
13121 
13123 
13124  /* get offsets */
13126  _rast[0],
13127  gt[0], gt[3],
13128  &(_offset[0][0]), &(_offset[0][1]),
13129  NULL
13130  ) != ES_NONE) {
13131  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the FIRST raster relative to the output raster");
13132  rt_raster_destroy(raster);
13133  return ES_ERROR;
13134  }
13135  _offset[0][0] *= -1;
13136  _offset[0][1] *= -1;
13137 
13139  _rast[1],
13140  gt[0], gt[3],
13141  &(_offset[1][0]), &(_offset[1][1]),
13142  NULL
13143  ) != ES_NONE) {
13144  rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the SECOND raster relative to the output raster");
13145  rt_raster_destroy(raster);
13146  return ES_ERROR;
13147  }
13148  _offset[1][0] *= -1;
13149  _offset[1][1] *= -1;
13150  break;
13151  }
13152  case ET_CUSTOM:
13153  rterror("rt_raster_from_two_rasters: Extent type ET_CUSTOM is not supported by this function");
13154  break;
13155  }
13156 
13157  /* set offsets if provided */
13158  if (NULL != offset) {
13159  for (i = 0; i < 4; i++)
13160  offset[i] = _offset[i / 2][i % 2];
13161  }
13162 
13163  *rtnraster = raster;
13164  return ES_NONE;
13165 }
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
tuple gt
Definition: window.py:79
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:123
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_api.c:5668
uint16_t height
Definition: rt_api.h:2227
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_api.c:6054
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_api.c:6005
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
Definition: rt_api.c:5442
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
uint16_t width
Definition: rt_api.h:2226
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster's geotransform using 6-element array.
Definition: rt_api.c:6026
void rterror(const char *fmt,...)
Raster core error and info handlers.
Definition: rt_api.c:895
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_api.c:6105
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_api.c:5353
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
Definition: rt_api.c:12774

Here is the call graph for this function:

Here is the caller graph for this function: