New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GRASS 7.0 NoData handling #4918

Closed
HuidaeCho opened this Issue May 1, 2014 · 5 comments

Comments

Projects
None yet
4 participants
@HuidaeCho
Copy link

HuidaeCho commented May 1, 2014

GRASS 7 uses NaN as a no data value for floating-point rasters and msGetGDALNoDataValue correctly returns NaN. When comparing pixel values with the NaN no data value, mapdrawgdal.c failes to find no data pixels because it simply uses the equal comparison. Because of this issue, dfScaleMin and dfScaleMax are set to a NaN and, as a result, COLORRANGE & DATARANGE don't work properly and an empty image is returned. A simple example is:

MAP
  SIZE 600 600
  EXTENT 2029205.71251009 1483554.15327001 2036686.24774584 1494117.98023943
  LAYER
    TYPE RASTER
    STATUS DEFAULT
    DATA "/home/grassdata/epsg102667/tmp/cellhd/contourareas"
    CLASS
      STYLE
        COLORRANGE 0 0 0 255 255 255 # black to white
        DATARANGE 702 768 # min value to max value
      END
    END
  END
END
@HuidaeCho

This comment has been minimized.

Copy link

HuidaeCho commented May 1, 2014

So here comes the patch:

--- mapdrawgdal.c.orig  2014-01-02 07:41:49.000000000 -0500
+++ mapdrawgdal.c       2014-04-30 23:50:31.632994038 -0400
@@ -49,6 +49,9 @@
 #include "gdal_alg.h"

 static int
+IsNoDataValue( double dfValue, double dfNoDataValue );
+
+static int
 LoadGDALImages( GDALDatasetH hDS, int band_numbers[4], int band_count,
                 layerObj *layer,
                 int src_xoff, int src_yoff, int src_xsize, int src_ysize,
@@ -1498,7 +1501,7 @@
       float fNoDataValue = (float) dfNoDataValue;

       for( i = 0; i < nPixelCount; i++ ) {
-        if( bGotNoData && pafRawData[i] == fNoDataValue )
+        if( bGotNoData && IsNoDataValue(pafRawData[i], fNoDataValue) )
           continue;

         if( !bMinMaxSet ) {
@@ -2138,7 +2141,7 @@
   bGotFirstValue = FALSE;

   for( i = 1; i < nPixelCount; i++ ) {
-    if( bGotNoData && pafRawData[i] == fNoDataValue )
+    if( bGotNoData && IsNoDataValue(pafRawData[i], fNoDataValue) )
       continue;

     if( !bGotFirstValue ) {
@@ -2301,7 +2304,7 @@
       /*
        * Skip nodata pixels ... no processing.
        */
-      if( bGotNoData && fRawValue == fNoDataValue ) {
+      if( bGotNoData && IsNoDataValue(fRawValue, fNoDataValue) ) {
         continue;
       }

@@ -2355,6 +2358,15 @@
 }

 /************************************************************************/
+/*                          IsNoDataValue()                             */
+/************************************************************************/
+static int
+IsNoDataValue( double dfValue, double dfNoDataValue )
+{
+  return (isnan(dfNoDataValue) && isnan(dfValue)) || dfValue == dfNoDataValue;
+}
+
+/************************************************************************/
 /*                          msGetGDALNoDataValue()                      */
 /************************************************************************/
 double
@tbonfort

This comment has been minimized.

Copy link
Member

tbonfort commented May 3, 2014

Patch seems correct, however Is NoDataValue() is called for each pixel in the input image which might cause some overhead. I would prefer if you could unroll those loops so the isNan calls only happen when dfNoDataValue is effectively NaN

@tbonfort tbonfort self-assigned this May 3, 2014

@tbonfort tbonfort added this to the 6.4.2 release milestone May 3, 2014

@HuidaeCho

This comment has been minimized.

Copy link

HuidaeCho commented May 3, 2014

You're right. I reduced isnan calls in half only for checking pixel values. Cheers!

--- mapdrawgdal.c.orig  2014-01-02 07:41:49.000000000 -0500
+++ mapdrawgdal.c       2014-05-03 07:30:47.730150463 -0400
@@ -49,6 +49,9 @@
 #include "gdal_alg.h"

 static int
+IsNoData( double dfValue, double dfNoDataValue, int bIsNoDataNan );
+
+static int
 LoadGDALImages( GDALDatasetH hDS, int band_numbers[4], int band_count,
                 layerObj *layer,
                 int src_xoff, int src_yoff, int src_xsize, int src_ysize,
@@ -1449,6 +1452,7 @@
     const char *pszScaleInfo;
     float *pafRawData;
     int   nPixelCount = dst_xsize * dst_ysize, i, bGotNoData = FALSE;
+    int bIsNoDataNan;
     GDALRasterBandH hBand =GDALGetRasterBand(hDS,band_numbers[iColorIndex]);
     pszScaleInfo = CSLFetchNameValue( layer->processing, "SCALE" );
     if( pszScaleInfo == NULL ) {
@@ -1489,6 +1493,7 @@
     pafRawData = pafWholeRawData + iColorIndex * dst_xsize * dst_ysize;

     dfNoDataValue = msGetGDALNoDataValue( layer, hBand, &bGotNoData );
+    bIsNoDataNan = isnan(dfNoDataValue);

     if( dfScaleMin == dfScaleMax ) {
       int bMinMaxSet = 0;
@@ -1498,7 +1503,7 @@
       float fNoDataValue = (float) dfNoDataValue;

       for( i = 0; i < nPixelCount; i++ ) {
-        if( bGotNoData && pafRawData[i] == fNoDataValue )
+        if( bGotNoData && IsNoData(pafRawData[i], fNoDataValue, bIsNoDataNan) )
           continue;

         if( !bMinMaxSet ) {
@@ -2088,6 +2093,7 @@
   const char *pszScaleInfo;
   const char *pszBuckets;
   int  *cmap, c, j, k, bGotNoData = FALSE, bGotFirstValue;
+  int bIsNoDataNan;
   unsigned char *rb_cmap[4];
   CPLErr eErr;
   rasterBufferObj *mask_rb = NULL;
@@ -2126,6 +2132,7 @@
   }

   fNoDataValue = (float) msGetGDALNoDataValue( layer, hBand, &bGotNoData );
+  bIsNoDataNan = isnan(fNoDataValue);

   /* ==================================================================== */
   /*      Determine scaling.                                              */
@@ -2138,7 +2145,7 @@
   bGotFirstValue = FALSE;

   for( i = 1; i < nPixelCount; i++ ) {
-    if( bGotNoData && pafRawData[i] == fNoDataValue )
+    if( bGotNoData && IsNoData(pafRawData[i], fNoDataValue, bIsNoDataNan) )
       continue;

     if( !bGotFirstValue ) {
@@ -2301,7 +2308,7 @@
       /*
        * Skip nodata pixels ... no processing.
        */
-      if( bGotNoData && fRawValue == fNoDataValue ) {
+      if( bGotNoData && IsNoData(fRawValue, fNoDataValue, bIsNoDataNan) ) {
         continue;
       }

@@ -2355,6 +2362,16 @@
 }

 /************************************************************************/
+/*                          IsNoData()                                  */
+/************************************************************************/
+static int
+IsNoData( double dfValue, double dfNoDataValue, int bIsNoDataNan )
+{
+  return (bIsNoDataNan && isnan(dfValue)) ||
+        (!bIsNoDataNan && dfValue == dfNoDataValue);
+}
+
+/************************************************************************/
 /*                          msGetGDALNoDataValue()                      */
 /************************************************************************/
 double
@sdlime

This comment has been minimized.

Copy link
Member

sdlime commented Jun 10, 2014

@rouault, any chance to you could review this patch and #4927 if you have a moment? This one seems more straight forward than the other.

@mapserver-bot

This comment has been minimized.

Copy link

mapserver-bot commented Feb 23, 2016

This is an automated comment

This issue has been closed due to lack of activity. This doesn't mean the issue is invalid, it simply got no attention within the last year. Please reopen with missing/relevant information if still valid.

Typically, issues fall in this state for one of the following reasons:

  • Hard, impossible or not enough information to reproduce
  • Missing test case
  • Lack of a champion with interest and/or funding to address the issue
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment