Rasterio and OpenCV reads JPEG differently

Very interesting issue...

I also checked on my side with (Python 3.6.9 on Ubuntu 18.04) and I can reproduce it:

rio.__version__: 1.1.8
cv2.__version__: 4.4.0
skimage.__version__: 0.17.2

Let me share my investigations.

First, two tiny practical differences: I reshaped rasterio image using their built-in tool as explained here:
https://rasterio.readthedocs.io/en/latest/topics/image_processing.html#imageorder

from rasterio.plot import reshape_as_image
im1 = reshape_as_image(im1)

Second, I also used an OpenCV parameter to convert the image to RGB as stated here: https://docs.opencv.org/master/d8/d01/group__imgproc__color__conversions.html

im2 = cv2.cvtColor(im2, cv2.COLOR_BGR2RGB)

Please, also note that using the cv2.IMREAD_UNCHANGED flag while reading the image with OpenCV doesn't change anything.

I finally double checked the dtypes:

ims = [im1, im2, im3]
for i in range(len(ims)):
    print("shape: {} and dtype: {}".format(ims[i].shape, ims[i].dtype))

The output is consistent:

shape: (1000, 1000, 3) and dtype: uint8
shape: (1000, 1000, 3) and dtype: uint8
shape: (1000, 1000, 3) and dtype: uint8

From here, I dig into the image using a plotting function:

def plotlist(lst, cmap=None):
    figsize = (12,12)
    fig, axs = plt.subplots(
        nrows=1,
        ncols=len(lst),
        figsize=figsize,
        dpi=100,
        sharex=True,
        sharey=True, )
    for i, ax in enumerate(axs):
        ax.imshow(lst[i], interpolation='none', cmap=cmap)

    plt.tight_layout()
    plt.show()

And used a super tiny subset of only 4 pixels (the four last on the right of the bottom row of the image):

subimgs = [im[-1:,-4:,:] for im in ims] # bottom right
plotlist(subimgs)

Which rendered:

Rendering of the 4 last pixels of the last row.

Pretty much the same right? Let's plot the differences:

deltas = []
deltas.append(np.subtract(subimgs[0],subimgs[1])) # rio - cv2
deltas.append(np.subtract(subimgs[0],subimgs[2])) # rio - skimage
deltas.append(np.subtract(subimgs[2],subimgs[1])) # skimage - cv2
plotlist(deltas)

Which rendered:

Rendering of the differences between the images on the 4 last pixels of the last row.

Whoops, something wrong with OpenCV? Yes, maybe, but from there you cannot be sure, it may also be the case that both rasterio and skimage are wrong in the same amount (tiny chances but...)

Let's make this clear by printing the actual values of these 4 last pixels:

>>> subimgs[0] # rio
array([[[58, 40, 38],
        [57, 37, 36],
        [63, 42, 39],
        [54, 30, 28]]], dtype=uint8)

>>> subimgs[2] # cv2
array([[[58, 40, 38],
        [56, 38, 36],
        [62, 42, 41],
        [51, 31, 30]]], dtype=uint8)

>>> subimgs[2] # skimage

array([[[58, 40, 38],
        [57, 37, 36],
        [63, 42, 39],
        [54, 30, 28]]], dtype=uint8)

Again, they look pretty much the same right? But wait... there actually are tiny little differences... (carefully look at each individual values...)

So, what is the truth?

This may (or may not) help, but using a third party tool can help decide. So I opened the image using Gimp and used the pipette on these 4 pixels, and here are the results:

Gimp
pixel 1 RGB: [58, 40, 38]
pixel 2 RGB: [56, 38, 36]
pixel 3 RGB: [62, 42, 41]
pixel 4 RGB: [51, 31, 30]

Yes, these values are the same than the one from OpenCV!
As a GIS person, I also decided to load the image into QGIS and query for the values on each of these pixels; results are the same than in Gimp and OpenCV.
As a last test, I also see what happens with Python GDAL:

from osgeo import gdal
print(gdal.__version__) # 3.1.0
ds = gdal.Open(path)
for i in range(3):
    i+=1
    print("{}".format(np.array(ds.GetRasterBand(i).ReadAsArray())[-1:,-4:]))

# Output (read pixel RGB values as columns here as each row represents a single band):
[[58 56 62 51]]
[[40 38 42 31]]
[[38 36 41 30]]

Again, the same than Gimp and OpenCV...

So now we have balanced the problem;
OpenCV == Gimp == QGIS == GDAL
rasterio == skimage

As I said, there are 2 possibilities from here;

1. maybe rasterio and skimage use the same lib under the hood and the artifact is hidden here?
2. maybe Gimp, QGIS, GDAL and OpenCV are based on the same lib (which I doubt) and the artifact is hidden here?

If someone know what is the dependency tree of all the image libraries used here, it could help fixing it.

I would intuitively say that the chances are higher that there is a glitch in skimage and that rasterio is, somehow, based on it.

My advice; consider also investigating in this direction instead of only OpenCV.