Computing vegetation indices using scikit-image: "Divide by zero" warning prevents Red Green Index calculation

Question

I am attempting to locate all green pixels in an image by converting vegetative indices using the scikit-image module in python. However, when I try to compute the Red-Green index, I get the following warning:

C:\Users\AppData\Local\Temp\ipykernel_1508\3901298948.py:3: RuntimeWarning: divide by zero encountered in divide
  RG_ratio = R/G
C:\Users\AppData\Local\Temp\ipykernel_1508\3901298948.py:3: RuntimeWarning: invalid value encountered in divide
  RG_ratio = R/G

The code below is what I have tried. Here is the image I am working with. It must be saved as a .jpg for this code to work:

Here is what I have done so far:

# Import modules
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from skimage.color import rgb2gray, label2rgb
from skimage.filters import threshold_otsu

# Read example image
filename = 'wheat.jpg'
RGB = mpimg.imread(filename)

# Display image
plt.axis('off')
plt.imshow(RGB)
plt.show()

# Convert image from float in the range 0-1 to unsigned integer in the range 0-255
RGB = (RGB*255).astype('uint8')
RGB.dtype

# Separate RGB channels

R = RGB[:, :, 0] #All the rows and all the columns for the first band (R)
G = RGB[:, :, 1]
B = RGB[:, :, 2]

# Display bands and histograms for each band
plt.figure(figsize=(12,8))
plt.subplot(2, 3, 1)
plt.axis('off')
plt.imshow(R, cmap='gray')
plt.title('Red')
plt.subplot(2, 3, 2)
plt.axis('off')
plt.imshow(G, cmap='gray')
plt.title('Green')
plt.subplot(2, 3, 3)
plt.axis('off')
plt.imshow(B, cmap='gray')
plt.title('Blue')
plt.subplot(2, 3, 4)
plt.axis('off')
plt.hist(R.flatten())
plt.subplot(2, 3, 5)
plt.axis('off')
plt.hist(G.flatten())
plt.subplot(2, 3, 6)
plt.axis('off')
plt.hist(B.flatten())
plt.show()

# Compute Red-Green Ratio
RG_ratio = R/G

# Compute Excess Green Index
ExG = 2*G - R - B

The warning occurs after running the RG_Ratio code

C:\Users\AppData\Local\Temp\ipykernel_1508\3901298948.py:3: RuntimeWarning: divide by zero encountered in divide
  RG_ratio = R/G
C:\Users\AppData\Local\Temp\ipykernel_1508\3901298948.py:3: RuntimeWarning: invalid value encountered in divide
  RG_ratio = R/G

Now, to try to produce a histogram of the RG_Ratio index:

# Plot histogram of RG_ratio
plt.figure()
plt.hist(RG_ratio.flatten(), bins='scott')
plt.show()

The following error is thrown:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [37], in <cell line: 3>()
      1 # Plot histogram of RG_ratio
      2 plt.figure()
----> 3 plt.hist(RG_ratio.flatten(), bins='scott')
      4 plt.show()

File ~\anaconda3\envs\agron893\lib\site-packages\matplotlib\pyplot.py:2600, in hist(x, bins, range, density, weights, cumulative, bottom, histtype, align, orientation, rwidth, log, color, label, stacked, data, **kwargs)
   2594 @_copy_docstring_and_deprecators(Axes.hist)
   2595 def hist(
   2596         x, bins=None, range=None, density=False, weights=None,
   2597         cumulative=False, bottom=None, histtype='bar', align='mid',
   2598         orientation='vertical', rwidth=None, log=False, color=None,
   2599         label=None, stacked=False, *, data=None, **kwargs):
-> 2600     return gca().hist(
   2601         x, bins=bins, range=range, density=density, weights=weights,
   2602         cumulative=cumulative, bottom=bottom, histtype=histtype,
   2603         align=align, orientation=orientation, rwidth=rwidth, log=log,
   2604         color=color, label=label, stacked=stacked,
   2605         **({"data": data} if data is not None else {}), **kwargs)

File ~\anaconda3\envs\agron893\lib\site-packages\matplotlib\__init__.py:1414, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1411 @functools.wraps(func)
   1412 def inner(ax, *args, data=None, **kwargs):
   1413     if data is None:
-> 1414         return func(ax, *map(sanitize_sequence, args), **kwargs)
   1416     bound = new_sig.bind(ax, *args, **kwargs)
   1417     auto_label = (bound.arguments.get(label_namer)
   1418                   or bound.kwargs.get(label_namer))

File ~\anaconda3\envs\agron893\lib\site-packages\matplotlib\axes\_axes.py:6641, in Axes.hist(self, x, bins, range, density, weights, cumulative, bottom, histtype, align, orientation, rwidth, log, color, label, stacked, **kwargs)
   6637 # Loop through datasets
   6638 for i in range(nx):
   6639     # this will automatically overwrite bins,
   6640     # so that each histogram uses the same bins
-> 6641     m, bins = np.histogram(x[i], bins, weights=w[i], **hist_kwargs)
   6642     tops.append(m)
   6643 tops = np.array(tops, float)  # causes problems later if it's an int

File <__array_function__ internals>:180, in histogram(*args, **kwargs)

File ~\anaconda3\envs\agron893\lib\site-packages\numpy\lib\histograms.py:793, in histogram(a, bins, range, normed, weights, density)
    681 r"""
    682 Compute the histogram of a dataset.
    683 
   (...)
    789 
    790 """
    791 a, weights = _ravel_and_check_weights(a, weights)
--> 793 bin_edges, uniform_bins = _get_bin_edges(a, bins, range, weights)
    795 # Histogram is an integer or a float array depending on the weights.
    796 if weights is None:

File ~\anaconda3\envs\agron893\lib\site-packages\numpy\lib\histograms.py:396, in _get_bin_edges(a, bins, range, weights)
    392 if weights is not None:
    393     raise TypeError("Automated estimation of the number of "
    394                     "bins is not supported for weighted data")
--> 396 first_edge, last_edge = _get_outer_edges(a, range)
    398 # truncate the range if needed
    399 if range is not None:

File ~\anaconda3\envs\agron893\lib\site-packages\numpy\lib\histograms.py:315, in _get_outer_edges(a, range)
    312         raise ValueError(
    313             'max must be larger than min in range parameter.')
    314     if not (np.isfinite(first_edge) and np.isfinite(last_edge)):
--> 315         raise ValueError(
    316             "supplied range of [{}, {}] is not finite".format(first_edge, last_edge))
    317 elif a.size == 0:
    318     # handle empty arrays. Can't determine range, so use 0-1.
    319     first_edge, last_edge = 0, 1

ValueError: supplied range of [0.0, inf] is not finite

I do want to compute the RG_ratio index, but I'm not sure what is causing the "divide by zero" warning. How can I fix this issue?

Answer 1

This link contains the solution.

The problem occurs in Numpy and not in scikit-learn.

You can disable the warning with numpy.seterr. Put this before the possible division by zero:

np.seterr(divide='ignore')

That'll disable zero division warnings globally. If you just want to disable them for a little bit, you can use numpy.errstate in a with clause:

with np.errstate(divide='ignore'):
    # some code here

For a zero by zero division (undetermined, results in a NaN), the error behaviour has changed with numpy version 1.12.0: this is now considered "invalid", while previously it was "divide".

Thus, if there is a chance you your numerator could be zero as well, use

np.seterr(divide='ignore', invalid='ignore')

or

with np.errstate(divide='ignore', invalid='ignore'):
    # some code here

Answer 2

Those warnings arise when some values in G are 0. This means that some values in R/G will be either inf (infinity) (if R is not 0 at that pixel) or nan (not a number) if R is also 0 at that pixel. You have several options:

• Make sure you have no 0 pixels before division. You can do this with R = np.clip(R, 1, 255) and G = np.clip(R, 1, 255). This would be my recommended approach.
• Ignore the inf and nan pixels later, when computing the histogram: plt.hist(RG_ratio[np.isfinite(RG_ratio)]). This will compute only the histogram of pixels that did not go to infinity or to nan.

This is a fundamental problem with the R/G ratio so I presume there is a standard approach in the GIS literature, you should find out what it is.