PostGIS 3.0.6dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ 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_free
void lwpoly_free(LWPOLY *poly)
Definition lwpoly.c:175
lwpoly_as_lwgeom
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition lwgeom.c:311
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
rt_util_envelope_to_lwpoly
LWPOLY * rt_util_envelope_to_lwpoly(rt_envelope ext)
Definition rt_util.c:439
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
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, 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, and rt_envelope::UpperLeftY.

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: