Skip to content

Commit cf5f99b

Browse files
KomzpaKomzpa
committed
raster: preserve map algebra nodata semantics
Keep map algebra outputs from reserving a NODATA value when source NODATA values are handled by the expression or callback path, while still reserving NODATA for synthesized extent gaps or unhandled source NODATA. Closes #2807 Closes #1100
1 parent c1e91f8 commit cf5f99b

6 files changed

Lines changed: 743 additions & 11 deletions

File tree

‎raster/rt_pg/rtpg_mapalgebra.c‎

Lines changed: 274 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -211,6 +211,10 @@ struct rtpg_nmapalgebra_arg_t {
211211
rtpg_nmapalgebra_callback_arg callback;
212212
};
213213

214+
static int rtpg_nmapalgebra_missing_band_can_emit_nodata(uint8_t *hasband, int numraster);
215+
static int
216+
rtpg_nmapalgebra_source_has_nodata(rt_raster *rasters, uint8_t *hasband, int *nband, int numraster, double *nodataval);
217+
214218
static rtpg_nmapalgebra_arg rtpg_nmapalgebra_arg_init(void) {
215219
rtpg_nmapalgebra_arg arg = NULL;
216220

@@ -278,6 +282,121 @@ static void rtpg_nmapalgebra_arg_destroy(rtpg_nmapalgebra_arg arg) {
278282
pfree(arg);
279283
}
280284

285+
static int
286+
rtpg_nmapalgebra_same_grid_extent(rt_raster rast1, rt_raster rast2)
287+
{
288+
int aligned = 0;
289+
290+
if (rast1 == NULL || rast2 == NULL || rt_raster_is_empty(rast1) || rt_raster_is_empty(rast2))
291+
return 0;
292+
293+
if (rt_raster_same_alignment(rast1, rast2, &aligned, NULL) != ES_NONE || !aligned)
294+
return 0;
295+
296+
return (rt_raster_get_width(rast1) == rt_raster_get_width(rast2) &&
297+
rt_raster_get_height(rast1) == rt_raster_get_height(rast2) &&
298+
FLT_EQ(rt_raster_get_x_offset(rast1), rt_raster_get_x_offset(rast2)) &&
299+
FLT_EQ(rt_raster_get_y_offset(rast1), rt_raster_get_y_offset(rast2)));
300+
}
301+
302+
static int
303+
rtpg_nmapalgebra_covers_grid_extent(rt_raster rast, rt_raster extent)
304+
{
305+
int aligned = 0;
306+
int covers = 0;
307+
308+
if (rast == NULL || extent == NULL || rt_raster_is_empty(rast) || rt_raster_is_empty(extent))
309+
return 0;
310+
311+
if (rt_raster_same_alignment(rast, extent, &aligned, NULL) != ES_NONE || !aligned)
312+
return 0;
313+
314+
if (rt_raster_covers(rast, -1, extent, -1, &covers) != ES_NONE)
315+
return 0;
316+
317+
return covers;
318+
}
319+
320+
static int
321+
rtpg_nmapalgebra_extent_can_synthesize_gaps(rt_raster *rasters,
322+
uint8_t *hasband,
323+
int numraster,
324+
rt_extenttype extenttype,
325+
rt_raster customextent)
326+
{
327+
rt_raster ref = NULL;
328+
int i = 0;
329+
330+
if (numraster < 2 || extenttype == ET_INTERSECTION)
331+
return 0;
332+
333+
switch (extenttype)
334+
{
335+
case ET_CUSTOM:
336+
ref = customextent;
337+
break;
338+
case ET_SECOND:
339+
ref = rasters[(numraster > 1) ? 1 : 0];
340+
break;
341+
case ET_LAST:
342+
ref = rasters[numraster - 1];
343+
break;
344+
case ET_UNION:
345+
for (i = 0; i < numraster; i++)
346+
{
347+
if (hasband[i])
348+
{
349+
ref = rasters[i];
350+
break;
351+
}
352+
}
353+
break;
354+
default:
355+
ref = rasters[0];
356+
break;
357+
}
358+
359+
if (ref == NULL)
360+
return 0;
361+
362+
for (i = 0; i < numraster; i++)
363+
{
364+
if (!hasband[i])
365+
continue;
366+
if (extenttype == ET_UNION)
367+
{
368+
if (!rtpg_nmapalgebra_same_grid_extent(ref, rasters[i]))
369+
return 1;
370+
}
371+
else if (!rtpg_nmapalgebra_covers_grid_extent(rasters[i], ref))
372+
return 1;
373+
}
374+
375+
return 0;
376+
}
377+
378+
static int
379+
rtpg_nmapalgebra_source_has_nodata(rt_raster *rasters, uint8_t *hasband, int *nband, int numraster, double *nodataval)
380+
{
381+
for (int i = 0; i < numraster; i++)
382+
{
383+
rt_band band;
384+
385+
if (!hasband[i])
386+
continue;
387+
388+
band = rt_raster_get_band(rasters[i], nband[i]);
389+
if (!rt_band_get_hasnodata_flag(band))
390+
continue;
391+
392+
if (nodataval)
393+
rt_band_get_nodata(band, nodataval);
394+
return 1;
395+
}
396+
397+
return 0;
398+
}
399+
281400
static int rtpg_nmapalgebra_rastbandarg_process(rtpg_nmapalgebra_arg arg, ArrayType *array, int *allnull, int *allempty, int *noband) {
282401
Oid etype;
283402
Datum *e;
@@ -642,6 +761,8 @@ Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS)
642761
int allnull = 0;
643762
int allempty = 0;
644763
int noband = 0;
764+
int source_has_nodata = 0;
765+
double source_nodataval = 0;
645766

646767
rt_raster raster = NULL;
647768
rt_band band = NULL;
@@ -1020,11 +1141,20 @@ Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS)
10201141
if (arg->pixtype == PT_END)
10211142
arg->pixtype = rt_band_get_pixtype(band);
10221143

1023-
/* set hasnodata and nodataval */
1024-
arg->hasnodata = rt_band_get_hasnodata_flag(band);
1144+
/* set output hasnodata and nodataval */
1145+
source_has_nodata = rtpg_nmapalgebra_source_has_nodata(
1146+
arg->raster, arg->hasband, arg->nband, arg->numraster, &source_nodataval);
1147+
arg->hasnodata = source_has_nodata;
1148+
if (!arg->hasnodata && rtpg_nmapalgebra_missing_band_can_emit_nodata(arg->hasband, arg->numraster))
1149+
arg->hasnodata = 1;
1150+
if (!arg->hasnodata && rtpg_nmapalgebra_extent_can_synthesize_gaps(
1151+
arg->raster, arg->hasband, arg->numraster, arg->extenttype, arg->cextent))
1152+
arg->hasnodata = 1;
10251153
arg->callback.hasnodata = arg->hasnodata;
1026-
if (arg->hasnodata)
1154+
if (rt_band_get_hasnodata_flag(band))
10271155
rt_band_get_nodata(band, &(arg->nodataval));
1156+
else if (source_has_nodata)
1157+
arg->nodataval = source_nodataval;
10281158
else
10291159
arg->nodataval = rt_band_get_min_value(band);
10301160

@@ -1173,6 +1303,128 @@ static void rtpg_nmapalgebraexpr_arg_destroy(rtpg_nmapalgebraexpr_arg arg) {
11731303
pfree(arg);
11741304
}
11751305

1306+
static int
1307+
rtpg_nmapalgebra_missing_band_can_emit_nodata(uint8_t *hasband, int numraster)
1308+
{
1309+
/* Missing requested bands are passed to the iterator as NODATA sources,
1310+
* even though no source band metadata can advertise that up front. */
1311+
for (int i = 0; i < numraster; i++)
1312+
{
1313+
if (!hasband[i])
1314+
return 1;
1315+
}
1316+
1317+
return 0;
1318+
}
1319+
1320+
static int
1321+
rtpg_nmapalgebraexpr_nodata_flags_can_emit(rtpg_nmapalgebraexpr_callback_arg *callback,
1322+
int *hasnodata,
1323+
int numraster,
1324+
int both_can_occur)
1325+
{
1326+
if (numraster > 1)
1327+
{
1328+
if (both_can_occur && hasnodata[0] && hasnodata[1] && !callback->nodatanodata.hasval)
1329+
return 1;
1330+
if (hasnodata[0] && !callback->expr[1].hasval && !callback->expr[1].spi_plan)
1331+
return 1;
1332+
if (hasnodata[1] && !callback->expr[2].hasval && !callback->expr[2].spi_plan)
1333+
return 1;
1334+
}
1335+
else if (hasnodata[0] && !callback->expr[1].hasval && !callback->expr[1].spi_plan)
1336+
return 1;
1337+
1338+
return 0;
1339+
}
1340+
1341+
static int
1342+
rtpg_nmapalgebraexpr_source_can_emit_nodata(rtpg_nmapalgebraexpr_callback_arg *callback,
1343+
rt_raster *rasters,
1344+
uint8_t *hasband,
1345+
int *nband,
1346+
int numraster)
1347+
{
1348+
int hasnodata[2] = {0, 0};
1349+
1350+
for (int i = 0; i < numraster && i < 2; i++)
1351+
{
1352+
rt_band band;
1353+
1354+
if (!hasband[i])
1355+
{
1356+
hasnodata[i] = 1;
1357+
continue;
1358+
}
1359+
1360+
band = rt_raster_get_band(rasters[i], nband[i]);
1361+
hasnodata[i] = rt_band_get_hasnodata_flag(band);
1362+
}
1363+
1364+
return rtpg_nmapalgebraexpr_nodata_flags_can_emit(callback, hasnodata, numraster, 1);
1365+
}
1366+
1367+
static int
1368+
rtpg_nmapalgebraexpr_extent_can_emit_nodata(rtpg_nmapalgebraexpr_callback_arg *callback,
1369+
rt_raster *rasters,
1370+
uint8_t *hasband,
1371+
int numraster,
1372+
rt_extenttype extenttype,
1373+
rt_raster customextent)
1374+
{
1375+
int gap[2] = {0, 0};
1376+
int both_can_occur = 0;
1377+
int refindex = 0;
1378+
1379+
if (numraster < 2 || extenttype == ET_INTERSECTION)
1380+
return 0;
1381+
1382+
switch (extenttype)
1383+
{
1384+
case ET_SECOND:
1385+
refindex = (numraster > 1) ? 1 : 0;
1386+
break;
1387+
case ET_LAST:
1388+
refindex = numraster - 1;
1389+
break;
1390+
default:
1391+
refindex = 0;
1392+
break;
1393+
}
1394+
1395+
for (int i = 0; i < numraster && i < 2; i++)
1396+
{
1397+
if (!hasband[i])
1398+
continue;
1399+
1400+
if (extenttype == ET_UNION)
1401+
{
1402+
for (int j = 0; j < numraster; j++)
1403+
{
1404+
if (i == j || !hasband[j])
1405+
continue;
1406+
if (!rtpg_nmapalgebra_covers_grid_extent(rasters[i], rasters[j]))
1407+
{
1408+
gap[i] = 1;
1409+
break;
1410+
}
1411+
}
1412+
}
1413+
else if (extenttype == ET_CUSTOM)
1414+
{
1415+
if (!rtpg_nmapalgebra_covers_grid_extent(rasters[i], customextent))
1416+
gap[i] = 1;
1417+
}
1418+
else if (i != refindex && !rtpg_nmapalgebra_covers_grid_extent(rasters[i], rasters[refindex]))
1419+
gap[i] = 1;
1420+
}
1421+
1422+
if (extenttype == ET_CUSTOM)
1423+
both_can_occur = 1;
1424+
1425+
return rtpg_nmapalgebraexpr_nodata_flags_can_emit(callback, gap, numraster, both_can_occur);
1426+
}
1427+
11761428
static int rtpg_nmapalgebraexpr_callback(
11771429
rt_iterator_arg arg, void *userarg,
11781430
double *value, int *nodata
@@ -1427,6 +1679,9 @@ Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS)
14271679
int allempty = 0;
14281680
int noband = 0;
14291681
int len = 0;
1682+
int source_can_emit_nodata = 0;
1683+
int source_has_nodata = 0;
1684+
double source_nodataval = 0;
14301685

14311686
TupleDesc tupdesc;
14321687
SPITupleTable *tuptable = NULL;
@@ -1725,11 +1980,24 @@ Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS)
17251980
if (arg->bandarg->pixtype == PT_END)
17261981
arg->bandarg->pixtype = rt_band_get_pixtype(band);
17271982

1728-
/* set hasnodata and nodataval */
1729-
arg->bandarg->hasnodata = rt_band_get_hasnodata_flag(band);
1983+
/* set output hasnodata and nodataval */
1984+
source_has_nodata = rtpg_nmapalgebra_source_has_nodata(
1985+
arg->bandarg->raster, arg->bandarg->hasband, arg->bandarg->nband, numraster, &source_nodataval);
1986+
source_can_emit_nodata = rtpg_nmapalgebraexpr_source_can_emit_nodata(
1987+
&(arg->callback), arg->bandarg->raster, arg->bandarg->hasband, arg->bandarg->nband, numraster);
1988+
arg->bandarg->hasnodata = source_can_emit_nodata;
1989+
if (!arg->bandarg->hasnodata && rtpg_nmapalgebraexpr_extent_can_emit_nodata(&(arg->callback),
1990+
arg->bandarg->raster,
1991+
arg->bandarg->hasband,
1992+
numraster,
1993+
arg->bandarg->extenttype,
1994+
arg->bandarg->cextent))
1995+
arg->bandarg->hasnodata = 1;
17301996
arg->callback.hasnodata = arg->bandarg->hasnodata;
1731-
if (arg->bandarg->hasnodata)
1997+
if (rt_band_get_hasnodata_flag(band))
17321998
rt_band_get_nodata(band, &(arg->bandarg->nodataval));
1999+
else if (source_can_emit_nodata && source_has_nodata)
2000+
arg->bandarg->nodataval = source_nodataval;
17332001
else
17342002
arg->bandarg->nodataval = rt_band_get_min_value(band);
17352003

0 commit comments

Comments
 (0)