PostGIS  3.4.0dev-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 991 of file rt_raster.c.

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

References covers(), ES_NONE, FLT_EQ, 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, pixval::x, and pixval::y.

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

Here is the call graph for this function:
Here is the caller graph for this function: