ValueError: zero-size array to reduction operation maximum which has no identity

[tobkug]

Hi I receive following error. The DEM is without any doubt of very bad quality but I'm interest what causes the error.

*** Iteration 1 ***

`Warping` all inputs to the following:
Resolution: 1.0
Extent: [780221.0, 2052737.0, 780622.0, 2053138.0]
Projection: '+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs '
Resampling alg: cubic

1 of 2: /media/ref.tif
2 of 2: [memory]
Elevation difference stats for uncorrected input DEMs (src - ref)
Removing outliers
Initial pixel count:
107630
Absolute dz filter: 100.00
107630
Excluding values outside of range: nan +/- 3*nan
Excluding values outside of range: nan to nan
107630
Computing slope
Slope filter: 0.10 - 40.00
Initial count: 65729
Excluding values outside of range: 0.100000 to 40.000000
37142
Computing aspect
Filtered difference map
/home/pygeotools/pygeotools/lib/malib.py:1665: UserWarning: Warning: converting a masked element to nan.
  print("count: %i min: %0.2f max: %0.2f mean: %0.2f std: %0.2f med: %0.2f mad: %0.2f q1: %0.2f q2: %0.2f iqr: %0.2f mode: %0.2f p16: %0.2f p84: %0.2f spread: %0.2f" % stats)
count: 37142 min: nan max: nan mean: nan std: nan med: nan mad: nan q1: nan q2: nan iqr: nan mode: -20.60 p16: nan p84: nan spread: nan
Computing sub-pixel offset between DEMs using mode: nuth
Computing common mask
Initial sample count:
37142
Removing outliers
3-sigma filter: nan - nan
0
Computing 1-degree bin statistics: count
Computing 1-degree bin statistics: median
Traceback (most recent call last):
  File "/home/demcoreg/demcoreg/dem_align.py", line 602, in <module>
    main()
  File "/home/demcoreg/demcoreg/dem_align.py", line 289, in main
    mask_list=mask_list, max_dz=max_dz, slope_lim=slope_lim, plot=True)
  File "/home/demcoreg/demcoreg/dem_align.py", line 147, in compute_offset
    fit_param, fig = coreglib.compute_offset_nuth(diff, slope, aspect, plot=plot)
  File "/home/demcoreg/demcoreg/coreglib.py", line 297, in compute_offset_nuth
    bin_ptp = np.cos(np.radians(bin_centers)).ptp()
  File "/home/.local/lib/python3.7/site-packages/numpy/core/_methods.py", line 247, in _ptp
    umr_maximum(a, axis, None, out, keepdims),
ValueError: zero-size array to reduction operation maximum which has no identity

`

[ShashankBice]

Check if no-data is properly defined and masked in your DEMs. The program fails to compute any valid stats (all nans), hence the minimizer/fit program complains. Can use replace_odv.py

[tobkug]

I used replace_ndv.py but the error remains. So the no-data should be properly defined.

[dshean]

Can you paste output of gdalinfo -stats for both input DEMs? It's possible that one uses nan as nodata, or that there is a nodata value defined in geotiff header, like 0 or -9999, but there are nan values in the geotiff.

There was a recent commit to replace_ndv.py to deal with this case (did you pull latest?). We should probably just add a check for nan in iolib functions in pygeotools, which should fix the issues you're experiencing with dem_align.py.

[tobkug]

after using replace_ndv.py with -9999 as new nodata value. I pulled the lastet.

Reference Dem:

Driver: GTiff/GeoTIFF
Files: /media/RefDem.tif
Size is 601, 601
Coordinate System is:
PROJCS["WGS_1984_UTM_Zone_18N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32618"]]
Origin = (780121.000000000000000,2053238.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  780121.000, 2053238.000) ( 72d20'46.72"W, 18d33' 4.42"N)
Lower Left  (  780121.000, 2052637.000) ( 72d20'47.03"W, 18d32'44.89"N)
Upper Right (  780722.000, 2053238.000) ( 72d20'26.24"W, 18d33' 4.14"N)
Lower Right (  780722.000, 2052637.000) ( 72d20'26.55"W, 18d32'44.60"N)
Center      (  780421.500, 2052937.500) ( 72d20'36.63"W, 18d32'54.51"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  Min=-43.781 Max=7.619 
  Minimum=-43.781, Maximum=7.619, Mean=-14.493, StdDev=4.976
  NoData Value=nan
  Metadata:
    RepresentationType=ATHEMATIC
    STATISTICS_COVARIANCES=24.75758567854982
    STATISTICS_MAXIMUM=7.6189703941345
    STATISTICS_MEAN=-14.493315708945
    STATISTICS_MINIMUM=-43.781005859375
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=4.9756995165052

Source Dem:

Driver: GTiff/GeoTIFF
Files: SrcDem.tif
Size is 622, 617
Coordinate System is:
PROJCS["WGS 84 / UTM zone 18N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32618"]]
Origin = (780117.000000000000000,2053254.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  780117.000, 2053254.000) ( 72d20'46.85"W, 18d33' 4.95"N)
Lower Left  (  780117.000, 2052637.000) ( 72d20'47.16"W, 18d32'44.89"N)
Upper Right (  780739.000, 2053254.000) ( 72d20'25.66"W, 18d33' 4.65"N)
Lower Right (  780739.000, 2052637.000) ( 72d20'25.97"W, 18d32'44.59"N)
Center      (  780428.000, 2052945.500) ( 72d20'36.41"W, 18d32'54.77"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Min=-143.445 Max=104.745 
  Minimum=-143.445, Maximum=104.745, Mean=-25.119, StdDev=11.048
  NoData Value=nan
  Metadata:
    STATISTICS_MAXIMUM=104.74529266357
    STATISTICS_MEAN=-25.119146955916
    STATISTICS_MINIMUM=-143.44519042969
    STATISTICS_STDDEV=11.047686537924
    STATISTICS_VALID_PERCENT=64.11

[dshean]

Hmm.

The above shows you still have nan as your nodata value in both DEMs. The latest commit on pygeotools will handle this case and replace with the input value: https://github.com/dshean/pygeotools/commit/9c8d7d5cff5f4fd1c08cf2d46dd1566d229ae9ad

When you run replace_ndv.py it writes out a new file with *_ndv.tif appended. That is the file with updated nodata values, not the original. Maybe you can try using those?

[tobkug]

ok, embarrassing. I didn't realize that a new tif. file was generated.

Thanks for the help and patience.

[dshean]

OK! Glad to hear we figured this out and you got something to work. The doc is not great, and I've made this same replace_ndv.py output file mistake myself from time to time. Seemed safer to create a new tif, to avoid risks of making some irreversible change to the original file.