PostGIS  2.1.10dev-r@@SVN_REVISION@@
rt_raster rt_raster_compute_skewed_raster ( rt_envelope  extent,
double *  skew,
double *  scale,
double  tolerance 
)

Definition at line 6764 of file rt_api.c.

References covers(), ES_NONE, FLT_EQ, rt_raster_t::height, LWGEOM2GEOS(), lwgeom_free(), lwgeom_geos_error(), lwnotice(), lwpoly_as_lwgeom(), lwpoly_free(), rt_envelope::MaxX, rt_envelope::MaxY, rt_envelope::MinX, rt_envelope::MinY, rtrowdump::raster, RASTER_DEBUG, RASTER_DEBUGF, rt_raster_cell_to_geopoint(), rt_raster_destroy(), rt_raster_geopoint_to_cell(), rt_raster_get_convex_hull(), rt_raster_new(), rt_raster_set_geotransform_matrix(), rt_raster_set_offsets(), rt_raster_set_scale(), rt_raster_set_skews(), rt_util_envelope_to_lwpoly(), rterror(), rt_envelope::UpperLeftX, rt_envelope::UpperLeftY, rt_raster_t::width, pixval::x, and pixval::y.

Referenced by rt_raster_gdal_rasterize(), rt_raster_gdal_warp(), and test_raster_compute_skewed_raster().

6769  {
6770  uint32_t run = 0;
6771  uint32_t max_run = 1;
6772  double dbl_run = 0;
6773 
6774  int rtn;
6775  int covers = 0;
6776  rt_raster raster;
6777  double _gt[6] = {0};
6778  double _igt[6] = {0};
6779  int _d[2] = {1, -1};
6780  int _dlast = 0;
6781  int _dlastpos = 0;
6782  double _w[2] = {0};
6783  double _r[2] = {0};
6784  double _xy[2] = {0};
6785  int i;
6786  int j;
6787  int x;
6788  int y;
6789 
6790  LWGEOM *geom = NULL;
6791  GEOSGeometry *sgeom = NULL;
6792  GEOSGeometry *ngeom = NULL;
6793 
6794  if (
6795  (tolerance < 0.) ||
6796  FLT_EQ(tolerance, 0.)
6797  ) {
6798  tolerance = 0.1;
6799  }
6800  else if (tolerance > 1.)
6801  tolerance = 1;
6802 
6803  dbl_run = tolerance;
6804  while (dbl_run < 10) {
6805  dbl_run *= 10.;
6806  max_run *= 10;
6807  }
6808 
6809  /* scale must be provided */
6810  if (scale == NULL)
6811  return NULL;
6812  for (i = 0; i < 2; i++) {
6813  if (FLT_EQ(scale[i], 0)) {
6814  rterror("rt_raster_compute_skewed_raster: Scale cannot be zero");
6815  return 0;
6816  }
6817 
6818  if (i < 1)
6819  _gt[1] = fabs(scale[i] * tolerance);
6820  else
6821  _gt[5] = fabs(scale[i] * tolerance);
6822  }
6823  /* conform scale-y to be negative */
6824  _gt[5] *= -1;
6825 
6826  /* skew not provided or skew is zero, return raster of correct dim and spatial attributes */
6827  if (
6828  (skew == NULL) || (
6829  FLT_EQ(skew[0], 0) &&
6830  FLT_EQ(skew[1], 0)
6831  )
6832  ) {
6833  int _dim[2] = {
6834  (int) fmax((fabs(extent.MaxX - extent.MinX) + (fabs(scale[0]) / 2.)) / fabs(scale[0]), 1),
6835  (int) fmax((fabs(extent.MaxY - extent.MinY) + (fabs(scale[1]) / 2.)) / fabs(scale[1]), 1)
6836  };
6837 
6838  raster = rt_raster_new(_dim[0], _dim[1]);
6839  if (raster == NULL) {
6840  rterror("rt_raster_compute_skewed_raster: Could not create output raster");
6841  return NULL;
6842  }
6843 
6844  rt_raster_set_offsets(raster, extent.MinX, extent.MaxY);
6845  rt_raster_set_scale(raster, fabs(scale[0]), -1 * fabs(scale[1]));
6846  rt_raster_set_skews(raster, skew[0], skew[1]);
6847 
6848  return raster;
6849  }
6850 
6851  /* direction to shift upper-left corner */
6852  if (skew[0] > 0.)
6853  _d[0] = -1;
6854  if (skew[1] < 0.)
6855  _d[1] = 1;
6856 
6857  /* geotransform */
6858  _gt[0] = extent.UpperLeftX;
6859  _gt[2] = skew[0] * tolerance;
6860  _gt[3] = extent.UpperLeftY;
6861  _gt[4] = skew[1] * tolerance;
6862 
6863  RASTER_DEBUGF(4, "Initial geotransform: %f, %f, %f, %f, %f, %f",
6864  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
6865  );
6866  RASTER_DEBUGF(4, "Delta: %d, %d", _d[0], _d[1]);
6867 
6868  /* simple raster */
6869  if ((raster = rt_raster_new(1, 1)) == NULL) {
6870  rterror("rt_raster_compute_skewed_raster: Out of memory allocating extent raster");
6871  return NULL;
6872  }
6873  rt_raster_set_geotransform_matrix(raster, _gt);
6874 
6875  /* get inverse geotransform matrix */
6876  if (!GDALInvGeoTransform(_gt, _igt)) {
6877  rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
6878  rt_raster_destroy(raster);
6879  return NULL;
6880  }
6881  RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
6882  _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
6883  );
6884 
6885  /* shift along axis */
6886  for (i = 0; i < 2; i++) {
6887  covers = 0;
6888  run = 0;
6889 
6890  RASTER_DEBUGF(3, "Shifting along %s axis", i < 1 ? "X" : "Y");
6891 
6892  do {
6893 
6894  /* prevent possible infinite loop */
6895  if (run > max_run) {
6896  rterror("rt_raster_compute_skewed_raster: Could not compute skewed extent due to check preventing infinite loop");
6897  rt_raster_destroy(raster);
6898  return NULL;
6899  }
6900 
6901  /*
6902  check the four corners that they are covered along the specific axis
6903  pixel column should be >= 0
6904  */
6905  for (j = 0; j < 4; j++) {
6906  switch (j) {
6907  /* upper-left */
6908  case 0:
6909  _xy[0] = extent.MinX;
6910  _xy[1] = extent.MaxY;
6911  break;
6912  /* lower-left */
6913  case 1:
6914  _xy[0] = extent.MinX;
6915  _xy[1] = extent.MinY;
6916  break;
6917  /* lower-right */
6918  case 2:
6919  _xy[0] = extent.MaxX;
6920  _xy[1] = extent.MinY;
6921  break;
6922  /* upper-right */
6923  case 3:
6924  _xy[0] = extent.MaxX;
6925  _xy[1] = extent.MaxY;
6926  break;
6927  }
6928 
6930  raster,
6931  _xy[0], _xy[1],
6932  &(_r[0]), &(_r[1]),
6933  _igt
6934  );
6935  if (rtn != ES_NONE) {
6936  rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
6937  rt_raster_destroy(raster);
6938  return NULL;
6939  }
6940 
6941  RASTER_DEBUGF(4, "Point %d at cell %d x %d", j, (int) _r[0], (int) _r[1]);
6942 
6943  /* raster doesn't cover point */
6944  if ((int) _r[i] < 0) {
6945  RASTER_DEBUGF(4, "Point outside of skewed extent: %d", j);
6946  covers = 0;
6947 
6948  if (_dlastpos != j) {
6949  _dlast = (int) _r[i];
6950  _dlastpos = j;
6951  }
6952  else if ((int) _r[i] < _dlast) {
6953  RASTER_DEBUG(4, "Point going in wrong direction. Reversing direction");
6954  _d[i] *= -1;
6955  _dlastpos = -1;
6956  run = 0;
6957  }
6958 
6959  break;
6960  }
6961 
6962  covers++;
6963  }
6964 
6965  if (!covers) {
6966  x = 0;
6967  y = 0;
6968  if (i < 1)
6969  x = _d[i] * fabs(_r[i]);
6970  else
6971  y = _d[i] * fabs(_r[i]);
6972 
6974  raster,
6975  x, y,
6976  &(_w[0]), &(_w[1]),
6977  _gt
6978  );
6979  if (rtn != ES_NONE) {
6980  rterror("rt_raster_compute_skewed_raster: Could not compute spatial coordinates for raster pixel");
6981  rt_raster_destroy(raster);
6982  return NULL;
6983  }
6984 
6985  /* adjust ul */
6986  if (i < 1)
6987  _gt[0] = _w[i];
6988  else
6989  _gt[3] = _w[i];
6990  rt_raster_set_geotransform_matrix(raster, _gt);
6991  RASTER_DEBUGF(4, "Shifted geotransform: %f, %f, %f, %f, %f, %f",
6992  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
6993  );
6994 
6995  /* get inverse geotransform matrix */
6996  if (!GDALInvGeoTransform(_gt, _igt)) {
6997  rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
6998  rt_raster_destroy(raster);
6999  return NULL;
7000  }
7001  RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
7002  _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
7003  );
7004  }
7005 
7006  run++;
7007  }
7008  while (!covers);
7009  }
7010 
7011  /* covers test */
7013  raster,
7014  extent.MaxX, extent.MinY,
7015  &(_r[0]), &(_r[1]),
7016  _igt
7017  );
7018  if (rtn != ES_NONE) {
7019  rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
7020  rt_raster_destroy(raster);
7021  return NULL;
7022  }
7023 
7024  RASTER_DEBUGF(4, "geopoint %f x %f at cell %d x %d", extent.MaxX, extent.MinY, (int) _r[0], (int) _r[1]);
7025 
7026  raster->width = _r[0];
7027  raster->height = _r[1];
7028 
7029  /* initialize GEOS */
7030  initGEOS(lwnotice, lwgeom_geos_error);
7031 
7032  /* create reference LWPOLY */
7033  {
7034  LWPOLY *npoly = rt_util_envelope_to_lwpoly(extent);
7035  if (npoly == NULL) {
7036  rterror("rt_raster_compute_skewed_raster: Could not build extent's geometry for covers test");
7037  rt_raster_destroy(raster);
7038  return NULL;
7039  }
7040 
7041  ngeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(npoly));
7042  lwpoly_free(npoly);
7043  }
7044 
7045  do {
7046  covers = 0;
7047 
7048  /* construct sgeom from raster */
7049  if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
7050  rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for covers test");
7051  GEOSGeom_destroy(ngeom);
7052  rt_raster_destroy(raster);
7053  return NULL;
7054  }
7055 
7056  sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom);
7057  lwgeom_free(geom);
7058 
7059  covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
7060  GEOSGeom_destroy(sgeom);
7061 
7062  if (covers == 2) {
7063  rterror("rt_raster_compute_skewed_raster: Could not run covers test");
7064  GEOSGeom_destroy(ngeom);
7065  rt_raster_destroy(raster);
7066  return NULL;
7067  }
7068 
7069  if (covers)
7070  break;
7071 
7072  raster->width++;
7073  raster->height++;
7074  }
7075  while (!covers);
7076 
7077  RASTER_DEBUGF(4, "Skewed extent does cover normal extent with dimensions %d x %d", raster->width, raster->height);
7078 
7079  raster->width = (int) ((((double) raster->width) * fabs(_gt[1]) + fabs(scale[0] / 2.)) / fabs(scale[0]));
7080  raster->height = (int) ((((double) raster->height) * fabs(_gt[5]) + fabs(scale[1] / 2.)) / fabs(scale[1]));
7081  _gt[1] = fabs(scale[0]);
7082  _gt[5] = -1 * fabs(scale[1]);
7083  _gt[2] = skew[0];
7084  _gt[4] = skew[1];
7085  rt_raster_set_geotransform_matrix(raster, _gt);
7086 
7087  /* minimize width/height */
7088  for (i = 0; i < 2; i++) {
7089  covers = 1;
7090  do {
7091  if (i < 1)
7092  raster->width--;
7093  else
7094  raster->height--;
7095 
7096  /* construct sgeom from raster */
7097  if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
7098  rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for minimizing dimensions");
7099  GEOSGeom_destroy(ngeom);
7100  rt_raster_destroy(raster);
7101  return NULL;
7102  }
7103 
7104  sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom);
7105  lwgeom_free(geom);
7106 
7107  covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
7108  GEOSGeom_destroy(sgeom);
7109 
7110  if (covers == 2) {
7111  rterror("rt_raster_compute_skewed_raster: Could not run covers test for minimizing dimensions");
7112  GEOSGeom_destroy(ngeom);
7113  rt_raster_destroy(raster);
7114  return NULL;
7115  }
7116 
7117  if (!covers) {
7118  if (i < 1)
7119  raster->width++;
7120  else
7121  raster->height++;
7122 
7123  break;
7124  }
7125  }
7126  while (covers);
7127  }
7128 
7129  GEOSGeom_destroy(ngeom);
7130 
7131  return raster;
7132 }
double MinY
Definition: rt_api.h:154
double UpperLeftX
Definition: rt_api.h:157
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
Datum covers(PG_FUNCTION_ARGS)
double MaxY
Definition: rt_api.h:155
double MaxX
Definition: rt_api.h:153
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1006
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:123
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom)
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:239
uint16_t height
Definition: rt_api.h:2227
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_api.c:5473
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
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:54
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
void lwgeom_geos_error(const char *fmt,...)
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:79
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
#define FLT_EQ(x, y)
Definition: rt_api.h:2159
tuple x
Definition: pixval.py:53
void rterror(const char *fmt,...)
Raster core error and info handlers.
Definition: rt_api.c:895
double MinX
Definition: rt_api.h:152
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_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
Definition: rt_api.c:6556
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_api.c:5353
LWPOLY * rt_util_envelope_to_lwpoly(rt_envelope env)
Definition: rt_api.c:543
double UpperLeftY
Definition: rt_api.h:158
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_api.c:5504
tuple y
Definition: pixval.py:54

Here is the call graph for this function:

Here is the caller graph for this function: