PostGIS  3.0.6dev-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 958 of file rt_raster.c.

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

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: