n-raster iterator.
Returns a raster with one band. The raster returned should be freed by the caller
838 {
839
841
843
844
846
848 int allnull = 0;
849 int allempty = 0;
850 int aligned = 0;
851 double offset[4] = {0.};
853
854 int i = 0;
855 int status = 0;
856 int inextent = 0;
859 int _x = 0;
860 int _y = 0;
861
862 int _width = 0;
863 int _height = 0;
864
865 double minval;
867 int isnodata;
868 int nodata;
869
871
872 assert(itrset != NULL && itrcount > 0);
873 assert(rtnraster != NULL);
874
875
876 *rtnraster = NULL;
877
878
879 if (callback == NULL) {
880 rterror(
"rt_raster_iterator: Callback function not provided");
882 }
883
884
886 rterror(
"rt_raster_iterator: Custom extent cannot be empty if extent type is ET_CUSTOM");
888 }
889
890
892 rterror(
"rt_raster_iterator: Pixel type cannot be PT_END");
894 }
895
896
898 rterror(
"rt_raster_iterator: Could not initialize internal variables");
900 }
901
902
904 rterror(
"rt_raster_iterator: Could not populate for internal variables");
907 }
908
909
910 if (allnull == itrcount) {
911 RASTER_DEBUG(3,
"all rasters are NULL, returning NULL");
912
914
916 }
917
918 else if (allempty == itrcount) {
919 RASTER_DEBUG(3,
"all rasters are empty, returning empty raster");
920
922
924 if (rtnrast == NULL) {
925 rterror(
"rt_raster_iterator: Could not create empty raster");
927 }
929
930 *rtnraster = rtnrast;
932 }
933
934
937
938
939
941 RASTER_DEBUG(4,
"using custom extent as reference raster");
943 }
944
945 else {
946 for (i = 0; i < itrcount; i++) {
948 RASTER_DEBUGF(4,
"using raster at index %d as reference raster", i);
950 break;
951 }
952 }
953 }
954
955
956 if (rast == NULL) {
957 rterror(
"rt_raster_iterator: Could not find reference raster to use for alignment tests");
958
960
962 }
963
964 do {
965 aligned = 1;
966
967
968 if (extenttype ==
ET_CUSTOM && rast != customextent) {
970 rterror(
"rt_raster_iterator: Could not test for alignment between reference raster and custom extent");
971
973
975 }
976
978 if (!aligned)
979 break;
980 }
981
982 for (i = 0; i < itrcount; i++) {
983
985 continue;
986
988 rterror(
"rt_raster_iterator: Could not test for alignment between reference raster and raster %d", i);
989
991
993 }
994 RASTER_DEBUGF(5,
"raster at index %d alignment: %d", i, aligned);
995
996
997 if (!aligned)
998 break;
999 }
1000 }
1001 while (0);
1002
1003
1004 if (!aligned) {
1005 rterror(
"rt_raster_iterator: The set of rasters provided (custom extent included, if appropriate) do not have the same alignment");
1006
1008
1010 }
1011
1012
1013 i = -1;
1014 switch (extenttype) {
1017
1019 if (rtnrast == NULL) {
1020 rterror(
"rt_raster_iterator: Could not allocate memory for output raster");
1021
1023
1025 }
1026
1027 for (i = 0; i < itrcount; i++) {
1030 break;
1031 }
1032 }
1034 rtnrast->
bands = NULL;
1035
1036
1038 for (i = i + 1; i < itrcount; i++) {
1040 continue;
1041
1044
1045 if (rast == NULL || status !=
ES_NONE) {
1046 rterror(
"rt_raster_iterator: Could not compute %s extent of rasters",
1047 extenttype ==
ET_UNION ?
"union" :
"intersection"
1048 );
1049
1051
1053 }
1055 rtinfo(
"rt_raster_iterator: Computed raster for %s extent is empty",
1056 extenttype ==
ET_UNION ?
"union" :
"intersection"
1057 );
1058
1060
1063 }
1064
1067 }
1068
1069 break;
1070
1071
1072
1073
1075 i = 0;
1076
1078 if (i < 0) {
1079 if (itrcount < 2)
1080 i = 0;
1081 else
1082 i = 1;
1083 }
1084
1086 if (i < 0) i = itrcount - 1;
1087
1088
1089 if (_param->
raster[i] == NULL) {
1090 RASTER_DEBUGF(3,
"returning NULL as %s raster is NULL and extent type is ET_%s",
1091 (i == 0 ? "first" : (i == 1 ? "second" : "last")),
1092 (i == 0 ? "FIRST" : (i == 1 ? "SECOND" : "LAST"))
1093 );
1094
1096
1098 }
1099
1100 else if (_param->
isempty[i]) {
1101 RASTER_DEBUGF(3,
"returning empty raster as %s raster is empty and extent type is ET_%s",
1102 (i == 0 ? "first" : (i == 1 ? "second" : "last")),
1103 (i == 0 ? "FIRST" : (i == 1 ? "SECOND" : "LAST"))
1104 );
1105
1107
1109 if (rtnrast == NULL) {
1110 rterror(
"rt_raster_iterator: Could not create empty raster");
1112 }
1114
1115 *rtnraster = rtnrast;
1117 }
1118
1119
1122 if (rtnrast == NULL) {
1123 rterror(
"rt_raster_iterator: Could not allocate memory for output raster");
1124
1126
1128 }
1129
1130 switch (extenttype) {
1132 memcpy(rtnrast, customextent,
sizeof(
struct rt_raster_t));
1133 break;
1134
1135 default:
1137 break;
1138 }
1140 rtnrast->
bands = NULL;
1141 break;
1142 }
1143
1146
1147 RASTER_DEBUGF(4,
"rtnrast (width, height, ulx, uly, scalex, scaley, skewx, skewy, srid) = (%d, %d, %f, %f, %f, %f, %f, %f, %d)",
1148 _width,
1149 _height,
1157 );
1158
1159
1161 rterror(
"rt_raster_iterator: Could not initialize empty values and NODATA");
1162
1165
1167 }
1168
1169
1171 rtnrast,
1172 pixtype,
1173 nodataval,
1174 hasnodata, nodataval,
1175 0
1176 ) < 0) {
1177 rterror(
"rt_raster_iterator: Could not add new band to output raster");
1178
1181
1183 }
1184
1185
1187 if (rtnband == NULL) {
1188 rterror(
"rt_raster_iterator: Could not get new band from output raster");
1189
1192
1194 }
1195
1196
1198
1199
1201 rterror(
"rt_raster_iterator: Could not initialize callback function argument");
1202
1206
1208 }
1209
1210
1211 for (i = 0; i < itrcount; i++) {
1213 continue;
1214
1218 rterror(
"rt_raster_iterator: Could not compute raster offsets");
1219
1223
1225 }
1226
1227 _param->
offset[i][0] = offset[2];
1228 _param->
offset[i][1] = offset[3];
1229 RASTER_DEBUGF(4,
"rast %d offset: %f %f", i, offset[2], offset[3]);
1230 }
1231
1232
1233
1234
1235 for (_y = 0; _y < _height; _y++) {
1236 for (_x = 0; _x < _width; _x++) {
1237 RASTER_DEBUGF(4,
"iterating output pixel (x, y) = (%d, %d)", _x, _y);
1240
1241
1242 for (i = 0; i < itrcount; i++) {
1244
1245
1246
1247
1248
1249
1250 if (
1252 (_param->
band.
rtband[i] == NULL && itrset[i].nbnodata) ||
1254 ) {
1255 RASTER_DEBUG(4,
"empty raster, band does not exist or band is NODATA. using empty values and NODATA");
1256
1259
1262
1263 continue;
1264 }
1265
1266
1267 x = _x - (int) _param->
offset[i][0];
1268 y = _y - (int) _param->
offset[i][1];
1270
1273
1274
1275 npixels = NULL;
1276 status = 0;
1277 if (distancex > 0 && distancey > 0) {
1279
1282 x, y,
1283 distancex, distancey,
1284 1,
1285 &npixels
1286 );
1287 if (status < 0) {
1288 rterror(
"rt_raster_iterator: Could not get pixel neighborhood");
1289
1293
1295 }
1296 }
1297
1298
1299
1300 if (
1301 (x >= 0 && x < _param->width[i]) &&
1302 (y >= 0 && y < _param->height[i])
1303 ) {
1307 x, y,
1308 &value,
1309 &isnodata
1311 rterror(
"rt_raster_iterator: Could not get the pixel value of band");
1312
1316
1318 }
1319 inextent = 1;
1320 }
1321
1322 else {
1323 RASTER_DEBUG(4,
"Outside band extent, setting value to NODATA");
1324
1327
1328 else
1330
1331 inextent = 0;
1332 isnodata = 1;
1333 }
1334
1335
1336 status++;
1337 if (status > 1)
1339 else
1341
1342 if (npixels == NULL) {
1343 rterror(
"rt_raster_iterator: Could not reallocate memory for neighborhood");
1344
1348
1350 }
1351
1352 npixels[status - 1].
x =
x;
1353 npixels[status - 1].
y =
y;
1354 npixels[status - 1].
nodata = 1;
1356
1357
1358 if ((!_param->
band.
hasnodata[i] && inextent) || !isnodata) {
1359 npixels[status - 1].
nodata = 0;
1360 }
1361 RASTER_DEBUGF(4,
"value, nodata: %f, %d", value, npixels[status - 1].nodata);
1362
1363
1365 npixels, status, mask,
1366 x, y,
1367 distancex, distancey,
1370 NULL, NULL
1371 );
1374 rterror(
"rt_raster_iterator: Could not create 2D array of neighborhood");
1375
1379
1381 }
1382 }
1383
1384
1387 nodata = 0;
1388 status = callback(_param->
arg, userarg, &value, &nodata);
1389
1390
1392
1393
1394 if (status == 0) {
1395 rterror(
"rt_raster_iterator: Callback function returned an error");
1396
1400
1402 }
1403
1404
1405 status = 0;
1406 if (!nodata) {
1408 RASTER_DEBUGF(4,
"burning pixel (%d, %d) with value: %f", _x, _y, value);
1409 }
1410 else if (!hasnodata) {
1412 RASTER_DEBUGF(4,
"burning pixel (%d, %d) with minval: %f", _x, _y, minval);
1413 }
1414 else {
1416 }
1418 rterror(
"rt_raster_iterator: Could not set pixel value");
1419
1423
1425 }
1426 }
1427 }
1428
1429
1431
1432 *rtnraster = rtnrast;
1434}
void rterror(const char *fmt,...) __attribute__((format(printf
Wrappers used for reporting errors and info.
void * rtalloc(size_t size)
Wrappers used for managing memory.
#define RASTER_DEBUG(level, msg)
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
#define RASTER_DEBUGF(level, msg,...)
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
void void rtinfo(const char *fmt,...) __attribute__((format(printf
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
struct rt_pixel_t * rt_pixel
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
void rt_band_destroy(rt_band band)
Destroy a raster band.
void * rtrealloc(void *mem, size_t size)
uint16_t rt_raster_get_height(rt_raster raster)
rt_errorstate rt_pixel_set_to_array(rt_pixel npixel, uint32_t count, rt_mask mask, int x, int y, uint16_t distancex, uint16_t distancey, double ***value, int ***nodata, int *dimx, int *dimy)
uint16_t rt_raster_get_width(rt_raster raster)
uint32_t rt_band_get_nearest_pixel(rt_band band, int x, int y, uint16_t distancex, uint16_t distancey, int exclude_nodata_value, rt_pixel *npixels)
Get nearest pixel(s) with value (not NODATA) to specified pixel.
void rtdealloc(void *mem)
rt_errorstate rt_raster_from_two_rasters(rt_raster rast1, rt_raster rast2, rt_extenttype extenttype, rt_raster *rtnraster, double *offset)
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
static int _rti_iterator_arg_callback_init(_rti_iterator_arg _param)
static int _rti_iterator_arg_empty_init(_rti_iterator_arg _param)
static _rti_iterator_arg _rti_iterator_arg_init()
static int _rti_iterator_arg_populate(_rti_iterator_arg _param, rt_iterator itrset, uint16_t itrcount, uint16_t distancex, uint16_t distancey, int *allnull, int *allempty)
static void _rti_iterator_arg_callback_clean(_rti_iterator_arg _param)
static void _rti_iterator_arg_destroy(_rti_iterator_arg _param)
struct _rti_iterator_arg_t::@16 band
struct _rti_iterator_arg_t::@19 empty