Fixed no data value mismanagement

Author: mauroalberti

main_dialog.py

@@ -391,6 +391,7 @@ def write_and_load_result(self, result_ga: GeoArray, result_fpath: str, load_in_
             return
 
         # add required layer to the map canvas - modified after RasterCalc module
+
         if load_in_TOC:
             newLayer = QgsRasterLayer(result_fpath, QFileInfo(result_fpath).baseName())
             QgsProject.instance().addMapLayer(newLayer)

pygsf/defaults/constants.py

@@ -1,3 +1,4 @@
 
 
-GRID_NULL_VALUE = -99999.99
+GRID_NULL_VALUE = -99999
+

pygsf/examples/load_band.py

@@ -15,7 +15,6 @@
 print("Data type: {}".format(band_params["dataType"]))
 print("Unit type: {}".format(band_params["unitType"]))
 print("Statistics: {}".format(band_params["stats"]))
-print("No data value: {}".format(band_params["noData"]))
 
 print("Number of overviews: {}".format(band_params["numOverviews"]))
 print("Number of color table entries: {}".format(band_params["numColorTableEntries"]))

pygsf/libs_utils/gdal/gdal.py

@@ -3,13 +3,14 @@
 
 from typing import Any, Tuple, Dict, Optional, Union
 
+from math import isfinite
+
 import numpy as np
 
 import gdal
 
 from .exceptions import *
 
-from ...defaults.constants import *
 from ...spatial.rasters.geotransform import GeoTransform
 
 
@@ -99,24 +100,19 @@ def read_band(dataset: gdal.Dataset, bnd_ndx: int=1) -> Tuple[dict, 'np.array']:
     if grid_values is None:
         raise RasterIOException("Unable to read data from rasters")
 
-    # transform data into numpy array
-
-    data = np.asarray(grid_values)
-
     # if nodatavalue exists, set null values to NaN in numpy array
 
-    if noDataVal is not None:
-        data = np.where(abs(data - noDataVal) > 1e-10, data, np.NaN)
+    if noDataVal is not None and isfinite(noDataVal):
+        grid_values = np.where(abs(grid_values - noDataVal) > 1e-10, grid_values, np.NaN)
 
     band_params = dict(
         dataType=data_type,
         unitType=unit_type,
         stats=dStats,
-        noData=noDataVal,
         numOverviews=nOverviews,
         numColorTableEntries=nColTableEntries)
 
-    return band_params, data
+    return band_params, grid_values
 
 
 def try_read_raster_band(raster_source: str, bnd_ndx: int=1) -> Tuple[bool, Union[str, Tuple[GeoTransform, str, Dict, 'np.array']]]:

pygsf/spatial/rasters/fields.py

@@ -168,7 +168,9 @@ def magnitude(
     Examples:
     """
 
-    return np.sqrt(fld_x ** 2 + fld_y ** 2)
+    magn = np.sqrt(fld_x ** 2 + fld_y ** 2)
+
+    return magn
 
 
 def orients_r(

pygsf/spatial/rasters/io.py

@@ -59,7 +59,7 @@ def try_write_esrigrid(geoarray: GeoArray, outgrid_flpth: str, esri_nullvalue: N
 
     arr = geoarray.level(level_ndx)
     if arr is None:
-        return False, "Array with index {} does not exist".format((level_ndx))
+        return False, "Array with index {} does not exist".format(level_ndx)
 
     num_rows, num_cols = arr.shape
     llc_x, llc_y = geoarray.level_llc(level_ndx)
@@ -71,7 +71,7 @@ def try_write_esrigrid(geoarray: GeoArray, outgrid_flpth: str, esri_nullvalue: N
     outputgrid.write("XLLCORNER %.8f\n" % llc_x)
     outputgrid.write("YLLCORNER %.8f\n" % llc_y)
     outputgrid.write("CELLSIZE %.8f\n" % cell_size_x)
-    outputgrid.write("NODATA_VALUE %f\n" % esri_nullvalue)
+    outputgrid.write("NODATA_VALUE %.8f\n" % esri_nullvalue)
 
     esrigrid_outvalues = np.where(np.isnan(arr), esri_nullvalue, arr)