PostGIS 3.7.0dev-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 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:1246
lwpoly_free
void lwpoly_free(LWPOLY *poly)
Definition lwpoly.c:175
lwpoly_as_lwgeom
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition lwgeom.c:357
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:304
rt_util_envelope_to_lwpoly
LWPOLY * rt_util_envelope_to_lwpoly(rt_envelope ext)
Definition rt_util.c:588
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:308
rtinfo
void void rtinfo(const char *fmt,...) __attribute__((format(printf
FLT_EQ
#define FLT_EQ(x, y)
Definition librtcore.h:2436
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:839
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:125
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:2486

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: