PostGIS  3.7.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 869 of file rt_raster.c.

874  {
875  uint32_t run = 0;
876  uint32_t max_run = 1;
877  double dbl_run = 0;
878 
879  int rtn;
880  int covers = 0;
881  rt_raster raster;
882  double _gt[6] = {0};
883  double _igt[6] = {0};
884  int _d[2] = {1, -1};
885  int _dlast = 0;
886  int _dlastpos = 0;
887  double _w[2] = {0};
888  double _r[2] = {0};
889  double _xy[2] = {0};
890  int i;
891  int j;
892  int x;
893  int y;
894 
895  LWGEOM *geom = NULL;
896  GEOSGeometry *sgeom = NULL;
897  GEOSGeometry *ngeom = NULL;
898 
899  if (
900  (tolerance < 0.) ||
901  FLT_EQ(tolerance, 0.)
902  ) {
903  tolerance = 0.1;
904  }
905  else if (tolerance > 1.)
906  tolerance = 1;
907 
908  dbl_run = tolerance;
909  while (dbl_run < 10) {
910  dbl_run *= 10.;
911  max_run *= 10;
912  }
913 
914  /* scale must be provided */
915  if (scale == NULL)
916  return NULL;
917  for (i = 0; i < 2; i++) {
918  if (FLT_EQ(scale[i], 0.0))
919  {
920  rterror("rt_raster_compute_skewed_raster: Scale cannot be zero");
921  return 0;
922  }
923 
924  if (i < 1)
925  _gt[1] = fabs(scale[i] * tolerance);
926  else
927  _gt[5] = fabs(scale[i] * tolerance);
928  }
929  /* conform scale-y to be negative */
930  _gt[5] *= -1;
931 
932  /* skew not provided or skew is zero, return raster of correct dim and spatial attributes */
933  if ((skew == NULL) || (FLT_EQ(skew[0], 0.0) && FLT_EQ(skew[1], 0.0)))
934  {
935  int _dim[2] = {
936  (int) fmax((fabs(extent.MaxX - extent.MinX) + (fabs(scale[0]) / 2.)) / fabs(scale[0]), 1),
937  (int) fmax((fabs(extent.MaxY - extent.MinY) + (fabs(scale[1]) / 2.)) / fabs(scale[1]), 1)
938  };
939 
940  raster = rt_raster_new(_dim[0], _dim[1]);
941  if (raster == NULL) {
942  rterror("rt_raster_compute_skewed_raster: Could not create output raster");
943  return NULL;
944  }
945 
946  rt_raster_set_offsets(raster, extent.MinX, extent.MaxY);
947  rt_raster_set_scale(raster, fabs(scale[0]), -1 * fabs(scale[1]));
948  rt_raster_set_skews(raster, skew[0], skew[1]);
949 
950  return raster;
951  }
952 
953  /* direction to shift upper-left corner */
954  if (skew[0] > 0.)
955  _d[0] = -1;
956  if (skew[1] < 0.)
957  _d[1] = 1;
958 
959  /* geotransform */
960  _gt[0] = extent.UpperLeftX;
961  _gt[2] = skew[0] * tolerance;
962  _gt[3] = extent.UpperLeftY;
963  _gt[4] = skew[1] * tolerance;
964 
965  RASTER_DEBUGF(4, "Initial geotransform: %f, %f, %f, %f, %f, %f",
966  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
967  );
968  RASTER_DEBUGF(4, "Delta: %d, %d", _d[0], _d[1]);
969 
970  /* simple raster */
971  if ((raster = rt_raster_new(1, 1)) == NULL) {
972  rterror("rt_raster_compute_skewed_raster: Out of memory allocating extent raster");
973  return NULL;
974  }
975  rt_raster_set_geotransform_matrix(raster, _gt);
976 
977  /* get inverse geotransform matrix */
978  if (!GDALInvGeoTransform(_gt, _igt)) {
979  rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
980  rt_raster_destroy(raster);
981  return NULL;
982  }
983  RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
984  _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
985  );
986 
987  /* shift along axis */
988  for (i = 0; i < 2; i++) {
989  covers = 0;
990  run = 0;
991 
992  RASTER_DEBUGF(3, "Shifting along %s axis", i < 1 ? "X" : "Y");
993 
994  do {
995 
996  /* prevent possible infinite loop */
997  if (run > max_run) {
998  rterror("rt_raster_compute_skewed_raster: Could not compute skewed extent due to check preventing infinite loop");
999  rt_raster_destroy(raster);
1000  return NULL;
1001  }
1002 
1003  /*
1004  check the four corners that they are covered along the specific axis
1005  pixel column should be >= 0
1006  */
1007  for (j = 0; j < 4; j++) {
1008  switch (j) {
1009  /* upper-left */
1010  case 0:
1011  _xy[0] = extent.MinX;
1012  _xy[1] = extent.MaxY;
1013  break;
1014  /* lower-left */
1015  case 1:
1016  _xy[0] = extent.MinX;
1017  _xy[1] = extent.MinY;
1018  break;
1019  /* lower-right */
1020  case 2:
1021  _xy[0] = extent.MaxX;
1022  _xy[1] = extent.MinY;
1023  break;
1024  /* upper-right */
1025  case 3:
1026  _xy[0] = extent.MaxX;
1027  _xy[1] = extent.MaxY;
1028  break;
1029  }
1030 
1031  rtn = rt_raster_geopoint_to_cell(
1032  raster,
1033  _xy[0], _xy[1],
1034  &(_r[0]), &(_r[1]),
1035  _igt
1036  );
1037  if (rtn != ES_NONE) {
1038  rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
1039  rt_raster_destroy(raster);
1040  return NULL;
1041  }
1042 
1043  RASTER_DEBUGF(4, "Point %d at cell %d x %d", j, (int) _r[0], (int) _r[1]);
1044 
1045  /* raster doesn't cover point */
1046  if ((int) _r[i] < 0) {
1047  RASTER_DEBUGF(4, "Point outside of skewed extent: %d", j);
1048  covers = 0;
1049 
1050  if (_dlastpos != j) {
1051  _dlast = (int) _r[i];
1052  _dlastpos = j;
1053  }
1054  else if ((int) _r[i] < _dlast) {
1055  RASTER_DEBUG(4, "Point going in wrong direction. Reversing direction");
1056  _d[i] *= -1;
1057  _dlastpos = -1;
1058  run = 0;
1059  }
1060 
1061  break;
1062  }
1063 
1064  covers++;
1065  }
1066 
1067  if (!covers) {
1068  x = 0;
1069  y = 0;
1070  if (i < 1)
1071  x = _d[i] * fabs(_r[i]);
1072  else
1073  y = _d[i] * fabs(_r[i]);
1074 
1075  rtn = rt_raster_cell_to_geopoint(
1076  raster,
1077  x, y,
1078  &(_w[0]), &(_w[1]),
1079  _gt
1080  );
1081  if (rtn != ES_NONE) {
1082  rterror("rt_raster_compute_skewed_raster: Could not compute spatial coordinates for raster pixel");
1083  rt_raster_destroy(raster);
1084  return NULL;
1085  }
1086 
1087  /* adjust ul */
1088  if (i < 1)
1089  _gt[0] = _w[i];
1090  else
1091  _gt[3] = _w[i];
1092  rt_raster_set_geotransform_matrix(raster, _gt);
1093  RASTER_DEBUGF(4, "Shifted geotransform: %f, %f, %f, %f, %f, %f",
1094  _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
1095  );
1096 
1097  /* get inverse geotransform matrix */
1098  if (!GDALInvGeoTransform(_gt, _igt)) {
1099  rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
1100  rt_raster_destroy(raster);
1101  return NULL;
1102  }
1103  RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
1104  _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
1105  );
1106  }
1107 
1108  run++;
1109  }
1110  while (!covers);
1111  }
1112 
1113  /* covers test */
1114  rtn = rt_raster_geopoint_to_cell(
1115  raster,
1116  extent.MaxX, extent.MinY,
1117  &(_r[0]), &(_r[1]),
1118  _igt
1119  );
1120  if (rtn != ES_NONE) {
1121  rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
1122  rt_raster_destroy(raster);
1123  return NULL;
1124  }
1125 
1126  RASTER_DEBUGF(4, "geopoint %f x %f at cell %d x %d", extent.MaxX, extent.MinY, (int) _r[0], (int) _r[1]);
1127 
1128  raster->width = _r[0];
1129  raster->height = _r[1];
1130 
1131  /* initialize GEOS */
1132  initGEOS(rtinfo, lwgeom_geos_error);
1133 
1134  /* create reference LWPOLY */
1135  {
1136  LWPOLY *npoly = rt_util_envelope_to_lwpoly(extent);
1137  if (npoly == NULL) {
1138  rterror("rt_raster_compute_skewed_raster: Could not build extent's geometry for covers test");
1139  rt_raster_destroy(raster);
1140  return NULL;
1141  }
1142 
1143  ngeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(npoly), 0);
1144  lwpoly_free(npoly);
1145  }
1146 
1147  do {
1148  covers = 0;
1149 
1150  /* construct sgeom from raster */
1151  if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
1152  rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for covers test");
1153  GEOSGeom_destroy(ngeom);
1154  rt_raster_destroy(raster);
1155  return NULL;
1156  }
1157 
1158  sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
1159  lwgeom_free(geom);
1160 
1161  covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
1162  GEOSGeom_destroy(sgeom);
1163 
1164  if (covers == 2) {
1165  rterror("rt_raster_compute_skewed_raster: Could not run covers test");
1166  GEOSGeom_destroy(ngeom);
1167  rt_raster_destroy(raster);
1168  return NULL;
1169  }
1170 
1171  if (!covers)
1172  {
1173  raster->width++;
1174  raster->height++;
1175  }
1176  }
1177  while (!covers);
1178 
1179  RASTER_DEBUGF(4, "Skewed extent does cover normal extent with dimensions %d x %d", raster->width, raster->height);
1180 
1181  raster->width = (int) ((((double) raster->width) * fabs(_gt[1]) + fabs(scale[0] / 2.)) / fabs(scale[0]));
1182  raster->height = (int) ((((double) raster->height) * fabs(_gt[5]) + fabs(scale[1] / 2.)) / fabs(scale[1]));
1183  _gt[1] = fabs(scale[0]);
1184  _gt[5] = -1 * fabs(scale[1]);
1185  _gt[2] = skew[0];
1186  _gt[4] = skew[1];
1187  rt_raster_set_geotransform_matrix(raster, _gt);
1188 
1189  /* minimize width/height */
1190  for (i = 0; i < 2; i++) {
1191  covers = 1;
1192  do {
1193  if (i < 1)
1194  raster->width--;
1195  else
1196  raster->height--;
1197 
1198  /* construct sgeom from raster */
1199  if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
1200  rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for minimizing dimensions");
1201  GEOSGeom_destroy(ngeom);
1202  rt_raster_destroy(raster);
1203  return NULL;
1204  }
1205 
1206  sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
1207  lwgeom_free(geom);
1208 
1209  covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
1210  GEOSGeom_destroy(sgeom);
1211 
1212  if (covers == 2) {
1213  rterror("rt_raster_compute_skewed_raster: Could not run covers test for minimizing dimensions");
1214  GEOSGeom_destroy(ngeom);
1215  rt_raster_destroy(raster);
1216  return NULL;
1217  }
1218  } while (covers);
1219 
1220  if (i < 1)
1221  raster->width++;
1222  else
1223  raster->height++;
1224 
1225  }
1226 
1227  GEOSGeom_destroy(ngeom);
1228 
1229  return raster;
1230 }
LWGEOM2GEOS
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
Definition: liblwgeom/lwgeom_geos.c:400
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:1218
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,...) __attribute__((format(printf
Wrappers used for reporting errors and info.
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:302
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:306
rtinfo
void void rtinfo(const char *fmt,...) __attribute__((format(printf
FLT_EQ
#define FLT_EQ(x, y)
Definition: librtcore.h:2424
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:491
covers
Datum covers(PG_FUNCTION_ARGS)
Definition: lwgeom_geos_predicates.c:722
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
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:637
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:609
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:686
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:2474

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: