PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ rt_raster_compute_skewed_raster()

rt_raster rt_raster_compute_skewed_raster ( rt_envelope  extent,
double *  skew,
double *  scale,
double  tolerance 
)

Definition at line 960 of file rt_raster.c.

References covers(), ES_NONE, FLT_EQ, rt_raster_t::height, LWGEOM2GEOS(), lwgeom_free(), lwgeom_geos_error(), 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(), rtinfo(), 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().

965  {
966  uint32_t run = 0;
967  uint32_t max_run = 1;
968  double dbl_run = 0;
969 
970  int rtn;
971  int covers = 0;
972  rt_raster raster;
973  double _gt[6] = {0};
974  double _igt[6] = {0};
975  int _d[2] = {1, -1};
976  int _dlast = 0;
977  int _dlastpos = 0;
978  double _w[2] = {0};
979  double _r[2] = {0};
980  double _xy[2] = {0};
981  int i;
982  int j;
983  int x;
984  int y;
985 
986  LWGEOM *geom = NULL;
987  GEOSGeometry *sgeom = NULL;
988  GEOSGeometry *ngeom = NULL;
989 
990  if (
991  (tolerance < 0.) ||
992  FLT_EQ(tolerance, 0.)
993  ) {
994  tolerance = 0.1;
995  }
996  else if (tolerance > 1.)
997  tolerance = 1;
998 
999  dbl_run = tolerance;
1000  while (dbl_run < 10) {
1001  dbl_run *= 10.;
1002  max_run *= 10;
1003  }
1004 
1005  /* scale must be provided */
1006  if (scale == NULL)
1007  return NULL;
1008  for (i = 0; i < 2; i++) {
1009  if (FLT_EQ(scale[i], 0)) {
1010  rterror("rt_raster_compute_skewed_raster: Scale cannot be zero");
1011  return 0;
1012  }
1013 
1014  if (i < 1)
1015  _gt[1] = fabs(scale[i] * tolerance);
1016  else
1017  _gt[5] = fabs(scale[i] * tolerance);
1018  }
1019  /* conform scale-y to be negative */
1020  _gt[5] *= -1;
1021 
1022  /* skew not provided or skew is zero, return raster of correct dim and spatial attributes */
1023  if (
1024  (skew == NULL) || (
1025  FLT_EQ(skew[0], 0) &&
1026  FLT_EQ(skew[1], 0)
1027  )
1028  ) {
1029  int _dim[2] = {
1030  (int) fmax((fabs(extent.MaxX - extent.MinX) + (fabs(scale[0]) / 2.)) / fabs(scale[0]), 1),
1031  (int) fmax((fabs(extent.MaxY - extent.MinY) + (fabs(scale[1]) / 2.)) / fabs(scale[1]), 1)
1032  };
1033 
1034  raster = rt_raster_new(_dim[0], _dim[1]);
1035  if (raster == NULL) {
1036  rterror("rt_raster_compute_skewed_raster: Could not create output raster");
1037  return NULL;
1038  }
1039 
1040  rt_raster_set_offsets(raster, extent.MinX, extent.MaxY);
1041  rt_raster_set_scale(raster, fabs(scale[0]), -1 * fabs(scale[1]));
1042  rt_raster_set_skews(raster, skew[0], skew[1]);
1043 
1044  return raster;
1045  }
1046 
1047  /* direction to shift upper-left corner */
1048  if (skew[0] > 0.)
1049  _d[0] = -1;
1050  if (skew[1] < 0.)
1051  _d[1] = 1;
1052 
1053  /* geotransform */
1054  _gt[0] = extent.UpperLeftX;
1055  _gt[2] = skew[0] * tolerance;
1056  _gt[3] = extent.UpperLeftY;
1057  _gt[4] = skew[1] * tolerance;
1058 
1059  RASTER_DEBUGF(4, "Initial geotransform: %f, %f, %f, %f, %f, %f",
1060  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
1061  );
1062  RASTER_DEBUGF(4, "Delta: %d, %d", _d[0], _d[1]);
1063 
1064  /* simple raster */
1065  if ((raster = rt_raster_new(1, 1)) == NULL) {
1066  rterror("rt_raster_compute_skewed_raster: Out of memory allocating extent raster");
1067  return NULL;
1068  }
1069  rt_raster_set_geotransform_matrix(raster, _gt);
1070 
1071  /* get inverse geotransform matrix */
1072  if (!GDALInvGeoTransform(_gt, _igt)) {
1073  rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
1074  rt_raster_destroy(raster);
1075  return NULL;
1076  }
1077  RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
1078  _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
1079  );
1080 
1081  /* shift along axis */
1082  for (i = 0; i < 2; i++) {
1083  covers = 0;
1084  run = 0;
1085 
1086  RASTER_DEBUGF(3, "Shifting along %s axis", i < 1 ? "X" : "Y");
1087 
1088  do {
1089 
1090  /* prevent possible infinite loop */
1091  if (run > max_run) {
1092  rterror("rt_raster_compute_skewed_raster: Could not compute skewed extent due to check preventing infinite loop");
1093  rt_raster_destroy(raster);
1094  return NULL;
1095  }
1096 
1097  /*
1098  check the four corners that they are covered along the specific axis
1099  pixel column should be >= 0
1100  */
1101  for (j = 0; j < 4; j++) {
1102  switch (j) {
1103  /* upper-left */
1104  case 0:
1105  _xy[0] = extent.MinX;
1106  _xy[1] = extent.MaxY;
1107  break;
1108  /* lower-left */
1109  case 1:
1110  _xy[0] = extent.MinX;
1111  _xy[1] = extent.MinY;
1112  break;
1113  /* lower-right */
1114  case 2:
1115  _xy[0] = extent.MaxX;
1116  _xy[1] = extent.MinY;
1117  break;
1118  /* upper-right */
1119  case 3:
1120  _xy[0] = extent.MaxX;
1121  _xy[1] = extent.MaxY;
1122  break;
1123  }
1124 
1125  rtn = rt_raster_geopoint_to_cell(
1126  raster,
1127  _xy[0], _xy[1],
1128  &(_r[0]), &(_r[1]),
1129  _igt
1130  );
1131  if (rtn != ES_NONE) {
1132  rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
1133  rt_raster_destroy(raster);
1134  return NULL;
1135  }
1136 
1137  RASTER_DEBUGF(4, "Point %d at cell %d x %d", j, (int) _r[0], (int) _r[1]);
1138 
1139  /* raster doesn't cover point */
1140  if ((int) _r[i] < 0) {
1141  RASTER_DEBUGF(4, "Point outside of skewed extent: %d", j);
1142  covers = 0;
1143 
1144  if (_dlastpos != j) {
1145  _dlast = (int) _r[i];
1146  _dlastpos = j;
1147  }
1148  else if ((int) _r[i] < _dlast) {
1149  RASTER_DEBUG(4, "Point going in wrong direction. Reversing direction");
1150  _d[i] *= -1;
1151  _dlastpos = -1;
1152  run = 0;
1153  }
1154 
1155  break;
1156  }
1157 
1158  covers++;
1159  }
1160 
1161  if (!covers) {
1162  x = 0;
1163  y = 0;
1164  if (i < 1)
1165  x = _d[i] * fabs(_r[i]);
1166  else
1167  y = _d[i] * fabs(_r[i]);
1168 
1169  rtn = rt_raster_cell_to_geopoint(
1170  raster,
1171  x, y,
1172  &(_w[0]), &(_w[1]),
1173  _gt
1174  );
1175  if (rtn != ES_NONE) {
1176  rterror("rt_raster_compute_skewed_raster: Could not compute spatial coordinates for raster pixel");
1177  rt_raster_destroy(raster);
1178  return NULL;
1179  }
1180 
1181  /* adjust ul */
1182  if (i < 1)
1183  _gt[0] = _w[i];
1184  else
1185  _gt[3] = _w[i];
1186  rt_raster_set_geotransform_matrix(raster, _gt);
1187  RASTER_DEBUGF(4, "Shifted geotransform: %f, %f, %f, %f, %f, %f",
1188  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
1189  );
1190 
1191  /* get inverse geotransform matrix */
1192  if (!GDALInvGeoTransform(_gt, _igt)) {
1193  rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
1194  rt_raster_destroy(raster);
1195  return NULL;
1196  }
1197  RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
1198  _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
1199  );
1200  }
1201 
1202  run++;
1203  }
1204  while (!covers);
1205  }
1206 
1207  /* covers test */
1208  rtn = rt_raster_geopoint_to_cell(
1209  raster,
1210  extent.MaxX, extent.MinY,
1211  &(_r[0]), &(_r[1]),
1212  _igt
1213  );
1214  if (rtn != ES_NONE) {
1215  rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
1216  rt_raster_destroy(raster);
1217  return NULL;
1218  }
1219 
1220  RASTER_DEBUGF(4, "geopoint %f x %f at cell %d x %d", extent.MaxX, extent.MinY, (int) _r[0], (int) _r[1]);
1221 
1222  raster->width = _r[0];
1223  raster->height = _r[1];
1224 
1225  /* initialize GEOS */
1226  initGEOS(rtinfo, lwgeom_geos_error);
1227 
1228  /* create reference LWPOLY */
1229  {
1230  LWPOLY *npoly = rt_util_envelope_to_lwpoly(extent);
1231  if (npoly == NULL) {
1232  rterror("rt_raster_compute_skewed_raster: Could not build extent's geometry for covers test");
1233  rt_raster_destroy(raster);
1234  return NULL;
1235  }
1236 
1237  ngeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(npoly), 0);
1238  lwpoly_free(npoly);
1239  }
1240 
1241  do {
1242  covers = 0;
1243 
1244  /* construct sgeom from raster */
1245  if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
1246  rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for covers test");
1247  GEOSGeom_destroy(ngeom);
1248  rt_raster_destroy(raster);
1249  return NULL;
1250  }
1251 
1252  sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
1253  lwgeom_free(geom);
1254 
1255  covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
1256  GEOSGeom_destroy(sgeom);
1257 
1258  if (covers == 2) {
1259  rterror("rt_raster_compute_skewed_raster: Could not run covers test");
1260  GEOSGeom_destroy(ngeom);
1261  rt_raster_destroy(raster);
1262  return NULL;
1263  }
1264 
1265  if (covers)
1266  break;
1267 
1268  raster->width++;
1269  raster->height++;
1270  }
1271  while (!covers);
1272 
1273  RASTER_DEBUGF(4, "Skewed extent does cover normal extent with dimensions %d x %d", raster->width, raster->height);
1274 
1275  raster->width = (int) ((((double) raster->width) * fabs(_gt[1]) + fabs(scale[0] / 2.)) / fabs(scale[0]));
1276  raster->height = (int) ((((double) raster->height) * fabs(_gt[5]) + fabs(scale[1] / 2.)) / fabs(scale[1]));
1277  _gt[1] = fabs(scale[0]);
1278  _gt[5] = -1 * fabs(scale[1]);
1279  _gt[2] = skew[0];
1280  _gt[4] = skew[1];
1281  rt_raster_set_geotransform_matrix(raster, _gt);
1282 
1283  /* minimize width/height */
1284  for (i = 0; i < 2; i++) {
1285  covers = 1;
1286  do {
1287  if (i < 1)
1288  raster->width--;
1289  else
1290  raster->height--;
1291 
1292  /* construct sgeom from raster */
1293  if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
1294  rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for minimizing dimensions");
1295  GEOSGeom_destroy(ngeom);
1296  rt_raster_destroy(raster);
1297  return NULL;
1298  }
1299 
1300  sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
1301  lwgeom_free(geom);
1302 
1303  covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
1304  GEOSGeom_destroy(sgeom);
1305 
1306  if (covers == 2) {
1307  rterror("rt_raster_compute_skewed_raster: Could not run covers test for minimizing dimensions");
1308  GEOSGeom_destroy(ngeom);
1309  rt_raster_destroy(raster);
1310  return NULL;
1311  }
1312 
1313  if (!covers) {
1314  if (i < 1)
1315  raster->width++;
1316  else
1317  raster->height++;
1318 
1319  break;
1320  }
1321  }
1322  while (covers);
1323  }
1324 
1325  GEOSGeom_destroy(ngeom);
1326 
1327  return raster;
1328 }
rt_envelope::MinY
double MinY
Definition: librtcore.h:167
rt_envelope::UpperLeftX
double UpperLeftX
Definition: librtcore.h:170
covers
Datum covers(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2023
rtrowdump.raster
raster
Be careful!! Zeros function&#39;s input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
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_envelope::MaxY
double MaxY
Definition: librtcore.h:168
rt_envelope::MaxX
double MaxX
Definition: librtcore.h:166
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1099
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
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2185
rterror
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
rt_util_envelope_to_lwpoly
LWPOLY * rt_util_envelope_to_lwpoly(rt_envelope ext)
Definition: rt_util.c:435
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:806
rt_raster_set_geotransform_matrix
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster&#39;s geotransform using 6-element array.
Definition: rt_raster.c:727
lwpoly_as_lwgeom
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:288
rt_raster_t::height
uint16_t height
Definition: librtcore.h:2253
rt_raster_get_convex_hull
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster&#39;s convex hull.
Definition: rt_geometry.c:803
uint32_t
unsigned int uint32_t
Definition: uthash.h:78
lwgeom_geos_error
void lwgeom_geos_error(const char *fmt,...)
Definition: liblwgeom/lwgeom_geos.c:42
lwpoly_free
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:174
rt_raster_t::width
uint16_t width
Definition: librtcore.h:2252
rtinfo
void rtinfo(const char *fmt,...)
Definition: rt_context.c:211
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
LWGEOM2GEOS
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
Definition: liblwgeom/lwgeom_geos.c:357
rt_raster_set_offsets
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
Definition: rt_raster.c:199
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_skews
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
Definition: rt_raster.c:168
rt_envelope::MinX
double MinX
Definition: librtcore.h:165
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
rt_raster_t
Definition: librtcore.h:2235
rt_envelope::UpperLeftY
double UpperLeftY
Definition: librtcore.h:171
rt_raster_destroy
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
Here is the call graph for this function:
Here is the caller graph for this function: