# Difference in NDVI Results Computed Via SNAP Tool and Rasterio Also there is difference when reading band values in SNAP and Rasterio

Hello Everyone,
I am getting difference in ndvi results computed via rasterio and SNAP.I am doing this on RAW imagery that I downloaded from https://dataspace.copernicus.eu/ .I am using SNAP 9.0 and I have used NDVI Processor to calculate NDVI on SNAP Tool.I have also noticed that SNAP reads bands values between 0 and 1 while rasterio uses UINT16 data type. I am attaching the screenshots below

Polygon = POLYGON((73.24819412786306 32.521843540039484,
73.25105892991965 32.521843540039484,
73.25105892991965 32.52284573593144,
73.24819412786306 32.52284573593144,
73.24819412786306 32.521843540039484))

Tile Name: S2A_MSIL2A_20230224T054821_N0509_R048_T43SCS_20230224T090324.SAFE

Rasterio code:

import rasterio
import numpy as np
from shapely.geometry import Polygon
import pyproj
from shapely.geometry import box

def calculate_ndvi(red_band, nir_band):
ndvi = (nir_band.astype(np.float64) - red_band.astype(np.float64)) / (nir_band.astype(np.float64) + red_band.astype(np.float64))
return ndvi

def get_polygon(coordinates, file_path):
with rasterio.open(file_path) as src:
original_polygon = Polygon([(coord[“g”], coord[“f”]) for coord in coordinates])

``````    # Create a transformer to convert polygon CRS to raster CRS
transformer = pyproj.Transformer.from_crs(pyproj.CRS('EPSG:4326'), src.crs, always_xy=True)

# Transform the original polygon to the raster CRS
transformed_polygon = []
for x, y in original_polygon.exterior.coords:
xx, yy = transformer.transform(x, y)
transformed_polygon.append((xx, yy))
transformed_polygon = Polygon(transformed_polygon)

# Mask the raster with the transformed polygon
out_img, out_transform = mask(src, [transformed_polygon], all_touched=True, crop=True)
out_img = out_img.astype(np.float64)/
out_img[out_img == 0.0] = np.nan

return out_img
``````

# Define the file paths for the red and NIR bands

red_band_file = r’D:\43SCS\43SCS_2023-02-24\S2A_MSIL2A_20230224T054821_N0509_R048_T43SCS_20230224T090324.SAFE\GRANULE\L2A_T43SCS_A040085_20230224T054823\IMG_DATA\R10m\T43SCS_20230224T054821_B04_10m.jp2’
nir_band_file = r’D:\43SCS\43SCS_2023-02-24\S2A_MSIL2A_20230224T054821_N0509_R048_T43SCS_20230224T090324.SAFE\GRANULE\L2A_T43SCS_A040085_20230224T054823\IMG_DATA\R10m\T43SCS_20230224T054821_B08_10m.jp2’

# Get the region of interest for both bands

coordinates = [{‘g’: 73.24819412786306, ‘f’: 32.521843540039484}, {‘g’: 73.24819412786306, ‘f’: 32.521843540039484}, {‘g’: 73.24819412786306, ‘f’: 32.521843540039484}, {‘g’: 73.24819412786306, ‘f’: 32.521843540039484}, {‘g’: 73.25105892991965, ‘f’: 32.521843540039484}, {‘g’: 73.25105892991965, ‘f’: 32.521843540039484}, {‘g’: 73.25105892991965, ‘f’: 32.521843540039484}, {‘g’: 73.25105892991965, ‘f’: 32.521843540039484}, {‘g’: 73.25105892991965, ‘f’: 32.52284573593144}, {‘g’: 73.25105892991965, ‘f’: 32.52284573593144}, {‘g’: 73.25105892991965, ‘f’: 32.52284573593144}, {‘g’: 73.25105892991965, ‘f’: 32.52284573593144}, {‘g’: 73.24819412786306, ‘f’: 32.52284573593144}, {‘g’: 73.24819412786306, ‘f’: 32.52284573593144}, {‘g’: 73.24819412786306, ‘f’: 32.52284573593144}, {‘g’: 73.24819412786306, ‘f’: 32.52284573593144}, {‘g’: 73.24819412786306, ‘f’: 32.521843540039484}, {‘g’: 73.24819412786306, ‘f’: 32.521843540039484}, {‘g’: 73.24819412786306, ‘f’: 32.521843540039484}, {‘g’: 73.24819412786306, ‘f’: 32.521843540039484}]
red_band_roi = get_polygon(coordinates, red_band_file)
nir_band_roi = get_polygon(coordinates, nir_band_file)

# Calculate NDVI

ndvi = calculate_ndvi(red_band_roi, nir_band_roi)

# Print NDVI statistics

print(“NDVI Statistics:”)
print(f"Max NDVI value Rasterio: {np.nanmax(ndvi)}“)
print(f"Min NDVI value Rasterio: {np.nanmin(ndvi)}”)

The results are computed on same polygon:
NDVI Statistics Rasterio:
Max NDVI value Rasterio: 0.6266440390326686
Min NDVI value Rasterio: 0.1268224485719992

NDVI Statistics SNAP:
Max NDVI value SNAP:0.8740970492362976
Min NDVI value SNAP:0.19548484683036804

Different data readings on Same Raster:
SNAP reading B4 data between 0 and 1 and there is some DL written there I donot know what that is:

And here is the rasterio reading same band 4 .It shows data type of UINT16

Unfortunately, no one answered this. I calculated NDVI in SNAP software with an NDVI processor and with BandMath both of them gave the same results. However, when I calculate NDVI in Python (or Matlab) after exporting the 12 bands (or only NIR and RED) I get underestimated values. Certainly, SNAP adds a correction during the calculation but I couldn’t find what it is. Poor documentation.

Probably the differences are caused by the case that Rasterio doesn’t account the offset and scaling. The data is provided as Digital Numbers which need to be scaled to reflectance values.

This can be found in the SentiWiki: S2 Products (copernicus.eu)