[raster] Support NODATA=NaN raster ingestion

Notable example of such raster is Facebook Population Density.

Closes #4412
Closes #407



git-svn-id: http://svn.osgeo.org/postgis/trunk@17465 b70326c6-7e19-0410-871a-916f4a2858ee

NEWS

@@ -152,6 +152,7 @@ PostGIS 3.0.0
   - #4225, Upgrade tiger to use tiger 2018 census files
   - #4198, Add ST_ConstrainedDelaunayTriangles SFCGAL function (Darafei
            Praliaskouski)
+  - #4412, Support ingesting rasters with NODATA=NaN (Darafei Praliaskouski)
 
 
 PostGIS 2.5.0

raster/loader/raster2pgsql.c

@@ -492,11 +492,8 @@ static void calc_tile_size(
 			}
 			r = r - (double) d;
 
-			if (
-				FLT_EQ(_r, -1) ||
-				(r < _r) ||
-				FLT_EQ(r, _r)
-			) {
+			if (FLT_EQ(_r, -1.0) || (r < _r) || FLT_EQ(r, _r))
+			{
 				/*_d = d;*/
 				_r = r;
 				_i = i;

raster/rt_core/librtcore.h

@@ -2228,12 +2228,13 @@ rt_util_hsv_to_rgb(
 );
 
 /*
-	helper macros for consistent floating point equality checks
+	helper macros for consistent floating point equality checks.
+	NaN equals NaN for NODATA support.
 */
-#define FLT_NEQ(x, y) (fabs(x - y) > FLT_EPSILON)
-#define FLT_EQ(x, y) (!FLT_NEQ(x, y))
-#define DBL_NEQ(x, y) (fabs(x - y) > DBL_EPSILON)
-#define DBL_EQ(x, y) (!DBL_NEQ(x, y))
+#define FLT_NEQ(x, y) ((x != y) && !(isnan(x) && isnan(y)) && (fabs(x - y) > FLT_EPSILON))
+#define FLT_EQ(x, y) ((x == y) || (isnan(x) && isnan(y)) || (fabs(x - y) <= FLT_EPSILON))
+#define DBL_NEQ(x, y) ((x != y) && !(isnan(x) && isnan(y)) && (fabs(x - y) > DBL_EPSILON))
+#define DBL_EQ(x, y) ((x == y) || (isnan(x) && isnan(y)) || (fabs(x - y) <= DBL_EPSILON))
 
 /*
 	helper macro for symmetrical rounding

raster/rt_core/rt_band.c

@@ -1755,11 +1755,11 @@ rt_band_check_is_nodata(rt_band band) {
 	int isnodata = 0;
 
 	assert(NULL != band);
+	band->isnodata = FALSE;
 
 	/* Check if band has nodata value */
 	if (!band->hasnodata) {
 		RASTER_DEBUG(3, "Band has no NODATA value");
-		band->isnodata = FALSE;
 		return FALSE;
 	}
 

raster/rt_core/rt_pixel.c

@@ -393,7 +393,8 @@ rt_errorstate rt_pixel_set_to_array(
 			/* unweighted (boolean) mask */
 			if (mask->weighted == 0) {
 				/* pixel is set to zero or nodata */
-				if (FLT_EQ(mask->values[_y][_x],0) || mask->nodata[_y][_x] == 1) {
+				if (FLT_EQ(mask->values[_y][_x], 0.0) || mask->nodata[_y][_x] == 1)
+				{
 					values[_y][_x] = 0;
 					nodatas[_y][_x] = 1;
 				}

raster/rt_core/rt_raster.c

@@ -767,10 +767,8 @@ rt_raster_cell_to_geopoint(
 		memcpy(_gt, gt, sizeof(double) * 6);
 
 	/* scale of matrix is not set */
-	if (
-		FLT_EQ(_gt[1], 0) ||
-		FLT_EQ(_gt[5], 0)
-	) {
+	if (FLT_EQ(_gt[1], 0.0) || FLT_EQ(_gt[5], 0.0))
+	{
 		rt_raster_get_geotransform_matrix(raster, _gt);
 	}
 
@@ -1006,7 +1004,8 @@ rt_raster_compute_skewed_raster(
 	if (scale == NULL)
 		return NULL;
 	for (i = 0; i < 2; i++) {
-		if (FLT_EQ(scale[i], 0)) {
+		if (FLT_EQ(scale[i], 0.0))
+		{
 			rterror("rt_raster_compute_skewed_raster: Scale cannot be zero");
 			return 0;
 		}
@@ -1020,12 +1019,8 @@ rt_raster_compute_skewed_raster(
 	_gt[5] *= -1;
 
 	/* skew not provided or skew is zero, return raster of correct dim and spatial attributes */
-	if (
-		(skew == NULL) || (
-			FLT_EQ(skew[0], 0) &&
-			FLT_EQ(skew[1], 0)
-		)
-	) {
+	if ((skew == NULL) || (FLT_EQ(skew[0], 0.0) && FLT_EQ(skew[1], 0.0)))
+	{
 		int _dim[2] = {
 			(int) fmax((fabs(extent.MaxX - extent.MinX) + (fabs(scale[0]) / 2.)) / fabs(scale[0]), 1),
 			(int) fmax((fabs(extent.MaxY - extent.MinY) + (fabs(scale[1]) / 2.)) / fabs(scale[1]), 1)
@@ -2641,12 +2636,8 @@ rt_raster_gdal_rasterize(
 		_scale[1] = fabs(*scale_y);
 	}
 	/* user-defined width/height */
-	else if (
-		(NULL != width) &&
-		(NULL != height) &&
-		(FLT_NEQ(*width, 0.0)) &&
-		(FLT_NEQ(*height, 0.0))
-	) {
+	else if ((NULL != width) && (NULL != height) && (*width != 0) && (*height != 0))
+	{
 		_dim[0] = abs(*width);
 		_dim[1] = abs(*height);
 
@@ -2851,10 +2842,8 @@ rt_raster_gdal_rasterize(
 	}
 
 	/* reprocess extent if skewed */
-	if (
-		FLT_NEQ(_skew[0], 0) ||
-		FLT_NEQ(_skew[1], 0)
-	) {
+	if (FLT_NEQ(_skew[0], 0.0) || FLT_NEQ(_skew[1], 0.0))
+	{
 		rt_raster skewedrast;
 
 		RASTER_DEBUG(3, "Computing skewed extent's envelope");
@@ -3118,7 +3107,7 @@ rt_raster_gdal_rasterize(
 			_gt[1] = *scale_x;
 
 			/* check for skew */
-			if (NULL != skew_x && FLT_NEQ(*skew_x, 0))
+			if (NULL != skew_x && FLT_NEQ(*skew_x, 0.0))
 				_gt[2] = *skew_x;
 		}
 		/* positive scale-y */
@@ -3148,7 +3137,7 @@ rt_raster_gdal_rasterize(
 			_gt[5] = *scale_y;
 
 			/* check for skew */
-			if (NULL != skew_y && FLT_NEQ(*skew_y, 0))
+			if (NULL != skew_y && FLT_NEQ(*skew_y, 0.0))
 				_gt[4] = *skew_y;
 		}
 	}

raster/rt_core/rt_spatial_relationship.c

@@ -875,10 +875,9 @@ int rt_raster_intersects_algorithm(
 								noval1 = 1;
 							}
 							/* cell is outside bounds of grid */
-							else if (
-								(Qr[pX] < 0 || Qr[pX] > width1 || FLT_EQ(Qr[pX], width1)) ||
-								(Qr[pY] < 0 || Qr[pY] > height1 || FLT_EQ(Qr[pY], height1))
-							) {
+							else if ((Qr[pX] < 0 || Qr[pX] >= width1) ||
+								 (Qr[pY] < 0 || Qr[pY] >= height1))
+							{
 								noval1 = 1;
 							}
 							else if (hasnodata1 == FALSE)
@@ -898,10 +897,9 @@ int rt_raster_intersects_algorithm(
 								noval2 = 1;
 							}
 							/* cell is outside bounds of grid */
-							else if (
-								(Qr[pX] < 0 || Qr[pX] > width2 || FLT_EQ(Qr[pX], width2)) ||
-								(Qr[pY] < 0 || Qr[pY] > height2 || FLT_EQ(Qr[pY], height2))
-							) {
+							else if ((Qr[pX] < 0 || Qr[pX] >= width2) ||
+								 (Qr[pY] < 0 || Qr[pY] >= height2))
+							{
 								noval2 = 1;
 							}
 							else if (hasnodata2 == FALSE)
@@ -1245,10 +1243,9 @@ rt_raster_intersects(
 								continue;
 							}
 
-							if (
-								(Qr[pX] < 0 || Qr[pX] > *widthL || FLT_EQ(Qr[pX], *widthL)) ||
-								(Qr[pY] < 0 || Qr[pY] > *heightL || FLT_EQ(Qr[pY], *heightL))
-							) {
+							if ((Qr[pX] < 0 || Qr[pX] >= *widthL) ||
+							    (Qr[pY] < 0 || Qr[pY] >= *heightL))
+							{
 								continue;
 							}
 

raster/rt_core/rt_util.c

@@ -77,7 +77,9 @@ rt_util_clamp_to_32BUI(double value) {
 
 float
 rt_util_clamp_to_32F(double value) {
-    return (float)fmin(fmax((value), -FLT_MAX), FLT_MAX);
+	if (isnan(value))
+		return value;
+	return (float)fmin(fmax((value), -FLT_MAX), FLT_MAX);
 }
 
 /**
@@ -650,7 +652,8 @@ rt_util_dbl_trunc_warning(
 #endif
 				result = 1;
 			}
-			else if (FLT_NEQ(checkvalint, initialvalue)) {
+			else if (checkvalint != initialvalue)
+			{
 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
 				rtwarn("Value set for %s band got truncated from %f to %d",
 					rt_pixtype_name(pixtype),
@@ -671,7 +674,8 @@ rt_util_dbl_trunc_warning(
 #endif
 				result = 1;
 			}
-			else if (FLT_NEQ(checkvaluint, initialvalue)) {
+			else if (checkvaluint != initialvalue)
+			{
 #if POSTGIS_RASTER_WARN_ON_TRUNCATION > 0
 				rtwarn("Value set for %s band got truncated from %f to %u",
 					rt_pixtype_name(pixtype),

raster/rt_core/rt_warp.c

@@ -286,11 +286,9 @@ rt_raster rt_raster_gdal_warp(
 			gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
 
 		/* substitute spatial info (lack of) with a real one EPSG:32731 (WGS84/UTM zone 31s) */
-		if (
-			FLT_EQ(gt[0], 0) && FLT_EQ(gt[3], 0) &&
-			FLT_EQ(gt[1], 1) && FLT_EQ(gt[5], -1) &&
-			FLT_EQ(gt[2], 0) && FLT_EQ(gt[4], 0)
-		) {
+		if (FLT_EQ(gt[0], 0.0) && FLT_EQ(gt[3], 0.0) && FLT_EQ(gt[1], 1.0) && FLT_EQ(gt[5], -1.0) &&
+		    FLT_EQ(gt[2], 0.0) && FLT_EQ(gt[4], 0.0))
+		{
 			double ngt[6] = {166021.4431, 0.1, 0, 10000000.0000, 0, -0.1};
 
 			rtwarn("Raster has default geotransform. Adjusting metadata for use of GDAL Warp API");
@@ -432,7 +430,8 @@ rt_raster rt_raster_gdal_warp(
 	}
 
 	/* scale not defined, use suggested */
-	if (FLT_EQ(_scale[0], 0) && FLT_EQ(_scale[1], 0)) {
+	if (FLT_EQ(_scale[0], 0.0) && FLT_EQ(_scale[1], 0.0))
+	{
 		_scale[0] = fabs(_gt[1]);
 		_scale[1] = fabs(_gt[5]);
 	}
@@ -472,10 +471,8 @@ rt_raster rt_raster_gdal_warp(
 	RASTER_DEBUGF(4, "Using skew: %f x %f", _skew[0], _skew[1]);
 
 	/* reprocess extent if skewed */
-	if (
-		FLT_NEQ(_skew[0], 0) ||
-		FLT_NEQ(_skew[1], 0)
-	) {
+	if (FLT_NEQ(_skew[0], 0.0) || FLT_NEQ(_skew[1], 0.0))
+	{
 		rt_raster skewedrast;
 
 		RASTER_DEBUG(3, "Computing skewed extent's envelope");
@@ -706,7 +703,7 @@ rt_raster rt_raster_gdal_warp(
 			_gt[1] = *scale_x;
 
 			/* check for skew */
-			if (NULL != skew_x && FLT_NEQ(*skew_x, 0))
+			if (NULL != skew_x && FLT_NEQ(*skew_x, 0.0))
 				_gt[2] = *skew_x;
 		}
 		/* positive scale-y */
@@ -730,7 +727,7 @@ rt_raster rt_raster_gdal_warp(
 			_gt[5] = *scale_y;
 
 			/* check for skew */
-			if (NULL != skew_y && FLT_NEQ(*skew_y, 0))
+			if (NULL != skew_y && FLT_NEQ(*skew_y, 0.0))
 				_gt[4] = *skew_y;
 		}
 	}

raster/rt_pg/rtpg_gdal.c

@@ -542,13 +542,15 @@ Datum RASTER_GDALWarp(PG_FUNCTION_ARGS)
 	/* scale x */
 	if (!PG_ARGISNULL(4)) {
 		scale[0] = PG_GETARG_FLOAT8(4);
-		if (FLT_NEQ(scale[0], 0)) scale_x = &scale[0];
+		if (FLT_NEQ(scale[0], 0.0))
+			scale_x = &scale[0];
 	}
 
 	/* scale y */
 	if (!PG_ARGISNULL(5)) {
 		scale[1] = PG_GETARG_FLOAT8(5);
-		if (FLT_NEQ(scale[1], 0)) scale_y = &scale[1];
+		if (FLT_NEQ(scale[1], 0.0))
+			scale_y = &scale[1];
 	}
 
 	/* grid alignment x */
@@ -566,13 +568,15 @@ Datum RASTER_GDALWarp(PG_FUNCTION_ARGS)
 	/* skew x */
 	if (!PG_ARGISNULL(8)) {
 		skew[0] = PG_GETARG_FLOAT8(8);
-		if (FLT_NEQ(skew[0], 0)) skew_x = &skew[0];
+		if (FLT_NEQ(skew[0], 0.0))
+			skew_x = &skew[0];
 	}
 
 	/* skew y */
 	if (!PG_ARGISNULL(9)) {
 		skew[1] = PG_GETARG_FLOAT8(9);
-		if (FLT_NEQ(skew[1], 0)) skew_y = &skew[1];
+		if (FLT_NEQ(skew[1], 0.0))
+			skew_y = &skew[1];
 	}
 
 	/* width */

raster/rt_pg/rtpg_geometry.c

@@ -826,13 +826,15 @@ Datum RASTER_asRaster(PG_FUNCTION_ARGS)
 	/* scale x */
 	if (!PG_ARGISNULL(1)) {
 		scale[0] = PG_GETARG_FLOAT8(1);
-		if (FLT_NEQ(scale[0], 0)) scale_x = &scale[0];
+		if (FLT_NEQ(scale[0], 0.0))
+			scale_x = &scale[0];
 	}
 
 	/* scale y */
 	if (!PG_ARGISNULL(2)) {
 		scale[1] = PG_GETARG_FLOAT8(2);
-		if (FLT_NEQ(scale[1], 0)) scale_y = &scale[1];
+		if (FLT_NEQ(scale[1], 0.0))
+			scale_y = &scale[1];
 	}
 	POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: scale (x, y) = %f, %f", scale[0], scale[1]);
 
@@ -1188,13 +1190,15 @@ Datum RASTER_asRaster(PG_FUNCTION_ARGS)
 	/* skewx */
 	if (!PG_ARGISNULL(12)) {
 		skew[0] = PG_GETARG_FLOAT8(12);
-		if (FLT_NEQ(skew[0], 0)) skew_x = &skew[0];
+		if (FLT_NEQ(skew[0], 0.0))
+			skew_x = &skew[0];
 	}
 
 	/* skewy */
 	if (!PG_ARGISNULL(13)) {
 		skew[1] = PG_GETARG_FLOAT8(13);
-		if (FLT_NEQ(skew[1], 0)) skew_y = &skew[1];
+		if (FLT_NEQ(skew[1], 0.0))
+			skew_y = &skew[1];
 	}
 	POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: skew (x, y) = %f, %f", skew[0], skew[1]);
 

raster/rt_pg/rtpg_mapalgebra.c

@@ -1895,11 +1895,8 @@ static int rtpg_union_mean_callback(
 	POSTGIS_RT_DEBUGF(4, "rast0: %f %d", arg->values[0][0][0], arg->nodata[0][0][0]);
 	POSTGIS_RT_DEBUGF(4, "rast1: %f %d", arg->values[1][0][0], arg->nodata[1][0][0]);
 
-	if (
-		!arg->nodata[0][0][0] &&
-		FLT_NEQ(arg->values[0][0][0], 0) &&
-		!arg->nodata[1][0][0]
-	) {
+	if (!arg->nodata[0][0][0] && FLT_NEQ(arg->values[0][0][0], 0.0) && !arg->nodata[1][0][0])
+	{
 		*value = arg->values[1][0][0] / arg->values[0][0][0];
 		*nodata = 0;
 	}

raster/rt_pg/rtpg_raster_properties.c

@@ -721,8 +721,8 @@ Datum RASTER_rasterToWorldCoord(PG_FUNCTION_ARGS)
 	}
 
 	/* raster skewed? */
-	skewed[0] = FLT_NEQ(rt_raster_get_x_skew(raster), 0) ? true : false;
-	skewed[1] = FLT_NEQ(rt_raster_get_y_skew(raster), 0) ? true : false;
+	skewed[0] = FLT_NEQ(rt_raster_get_x_skew(raster), 0.0) ? true : false;
+	skewed[1] = FLT_NEQ(rt_raster_get_y_skew(raster), 0.0) ? true : false;
 
 	/* column and row */
 	for (i = 1; i <= 2; i++) {
@@ -816,9 +816,9 @@ Datum RASTER_worldToRasterCoord(PG_FUNCTION_ARGS)
 	}
 
 	/* raster skewed? */
-	skewed = FLT_NEQ(rt_raster_get_x_skew(raster), 0) ? true : false;
+	skewed = FLT_NEQ(rt_raster_get_x_skew(raster), 0.0) ? true : false;
 	if (!skewed)
-		skewed = FLT_NEQ(rt_raster_get_y_skew(raster), 0) ? true : false;
+		skewed = FLT_NEQ(rt_raster_get_y_skew(raster), 0.0) ? true : false;
 
 	/* longitude and latitude */
 	for (i = 1; i <= 2; i++) {

raster/test/regress/tickets.sql

@@ -123,3 +123,5 @@ SELECT '#4102.1', ST_BandNoDataValue(ST_AddBand(ST_MakeEmptyRaster(2, 2, 0, 0, 1
 SELECT '#4102.2', ST_BandNoDataValue(ST_AddBand(ST_MakeEmptyRaster(2, 2, 0, 0, 1, -1, 0, 0, 0), 1, '32BSI', 0, -10), 1) AS rast;
 
 select '#3457', ST_Area((ST_DumpAsPolygons(ST_Clip(ST_ASRaster(ST_GeomFromText('POLYGON((0 0,100 0,100 100,0 100,0 0))',4326),ST_Addband(ST_MakeEmptyRaster(1,1,0,0,1,-1,0,0,4326),'32BF'::text,0,-1),'32BF'::text,1,-1), ST_GeomFromText('POLYGON((0 0,100 100,100 0,0 0))',4326)))).geom);
+
+SELECT '#4412', exists(select ST_PixelAsPolygons(ST_AddBand(ST_MakeEmptyRaster(2, 2, 0, 0, 1, -1, 0, 0, 0), 1, '64BF', 0, 'NaN'), 1)) AS rast;

raster/test/regress/tickets_expected

@@ -17,3 +17,4 @@ NOTICE:  Invalid band index (must use 1-based). Returning NULL
 #4102.1|-10
 #4102.2|-10
 #3457|4950
+#4412|t