[postgis-tickets] r17465 - [raster] Support NODATA=NaN raster ingestion
Darafei
komzpa at gmail.com
Sun Jun 2 07:16:29 PDT 2019
Author: komzpa
Date: 2019-06-02 19:16:28 -0700 (Sun, 02 Jun 2019)
New Revision: 17465
Modified:
trunk/NEWS
trunk/raster/loader/raster2pgsql.c
trunk/raster/rt_core/librtcore.h
trunk/raster/rt_core/rt_band.c
trunk/raster/rt_core/rt_pixel.c
trunk/raster/rt_core/rt_raster.c
trunk/raster/rt_core/rt_spatial_relationship.c
trunk/raster/rt_core/rt_util.c
trunk/raster/rt_core/rt_warp.c
trunk/raster/rt_pg/rtpg_gdal.c
trunk/raster/rt_pg/rtpg_geometry.c
trunk/raster/rt_pg/rtpg_mapalgebra.c
trunk/raster/rt_pg/rtpg_raster_properties.c
trunk/raster/test/regress/tickets.sql
trunk/raster/test/regress/tickets_expected
Log:
[raster] Support NODATA=NaN raster ingestion
Notable example of such raster is Facebook Population Density.
Closes #4412
Closes https://github.com/postgis/postgis/pull/407
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/NEWS 2019-06-03 02:16:28 UTC (rev 17465)
@@ -152,6 +152,7 @@
- #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
Modified: trunk/raster/loader/raster2pgsql.c
===================================================================
--- trunk/raster/loader/raster2pgsql.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/loader/raster2pgsql.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -492,11 +492,8 @@
}
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;
Modified: trunk/raster/rt_core/librtcore.h
===================================================================
--- trunk/raster/rt_core/librtcore.h 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_core/librtcore.h 2019-06-03 02:16:28 UTC (rev 17465)
@@ -2228,12 +2228,13 @@
);
/*
- 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
Modified: trunk/raster/rt_core/rt_band.c
===================================================================
--- trunk/raster/rt_core/rt_band.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_core/rt_band.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -1755,11 +1755,11 @@
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;
}
Modified: trunk/raster/rt_core/rt_pixel.c
===================================================================
--- trunk/raster/rt_core/rt_pixel.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_core/rt_pixel.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -393,7 +393,8 @@
/* 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;
}
Modified: trunk/raster/rt_core/rt_raster.c
===================================================================
--- trunk/raster/rt_core/rt_raster.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_core/rt_raster.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -767,10 +767,8 @@
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 @@
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 @@
_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 @@
_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 @@
}
/* 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 @@
_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 @@
_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;
}
}
Modified: trunk/raster/rt_core/rt_spatial_relationship.c
===================================================================
--- trunk/raster/rt_core/rt_spatial_relationship.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_core/rt_spatial_relationship.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -875,10 +875,9 @@
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 @@
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 @@
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;
}
Modified: trunk/raster/rt_core/rt_util.c
===================================================================
--- trunk/raster/rt_core/rt_util.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_core/rt_util.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -77,7 +77,9 @@
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 @@
#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 @@
#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),
Modified: trunk/raster/rt_core/rt_warp.c
===================================================================
--- trunk/raster/rt_core/rt_warp.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_core/rt_warp.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -286,11 +286,9 @@
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 @@
}
/* 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 @@
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 @@
_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 @@
_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;
}
}
Modified: trunk/raster/rt_pg/rtpg_gdal.c
===================================================================
--- trunk/raster/rt_pg/rtpg_gdal.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_pg/rtpg_gdal.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -542,13 +542,15 @@
/* 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 @@
/* 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 */
Modified: trunk/raster/rt_pg/rtpg_geometry.c
===================================================================
--- trunk/raster/rt_pg/rtpg_geometry.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_pg/rtpg_geometry.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -826,13 +826,15 @@
/* 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 @@
/* 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]);
Modified: trunk/raster/rt_pg/rtpg_mapalgebra.c
===================================================================
--- trunk/raster/rt_pg/rtpg_mapalgebra.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_pg/rtpg_mapalgebra.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -1895,11 +1895,8 @@
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;
}
Modified: trunk/raster/rt_pg/rtpg_raster_properties.c
===================================================================
--- trunk/raster/rt_pg/rtpg_raster_properties.c 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/rt_pg/rtpg_raster_properties.c 2019-06-03 02:16:28 UTC (rev 17465)
@@ -721,8 +721,8 @@
}
/* 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 @@
}
/* 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++) {
Modified: trunk/raster/test/regress/tickets.sql
===================================================================
--- trunk/raster/test/regress/tickets.sql 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/test/regress/tickets.sql 2019-06-03 02:16:28 UTC (rev 17465)
@@ -123,3 +123,5 @@
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;
Modified: trunk/raster/test/regress/tickets_expected
===================================================================
--- trunk/raster/test/regress/tickets_expected 2019-06-02 18:43:17 UTC (rev 17464)
+++ trunk/raster/test/regress/tickets_expected 2019-06-03 02:16:28 UTC (rev 17465)
@@ -17,3 +17,4 @@
#4102.1|-10
#4102.2|-10
#3457|4950
+#4412|t
More information about the postgis-tickets
mailing list