Formula Dark Object Subtraction 1 QGIS (Semi-automatic Classification Plugin)

Question

I am using QGIS with the plugin Semi-Automatic Classification Plugin to apply DOS1 atmospheric correction to bands 3 and 8 of Sentinel-2. Example:

I want to automate this task with a python script but I don't know which formula is using this plugin.

I've got similar results by getting the smallest pixel (ignoring no-data value) for every band and substract this pixel to all pixels from that band:

inFile = 'band3.jp2'
with rasterio.open(inFile) as r:
    values = r.read(1)
    values[values==0]=65535 #Ignore no-data values
    minValue = np.min(values)
    values = values - minValue
    values[values==(65535-minValue)]=0 #Set no-data values again to 0

But the result is not exactly the same. I have found the source code of the plugin but I don't understand it at all I just want to apply the same DOS1 formula with a simple script in python.

Answer 1

The most similar formula I found was applying this:

DN = DN - (LDNm – 100)

DN: the value of the pixel

LDNm: radiance that corresponds to a digital count value for which the sum of all the pixels with digital counts lower or equal to this value is equal to the 0.01% of all the pixels from the image considered

Python code to calculate LDNm:

def findDNmin(inputRaster):
    uniques, counts = np.unique(inputRaster, return_counts=True)
    sumTot = np.sum(inputRaster)*0.0001
    sumPartial=0
    index=0
    #uniques values are already sorted
    for value in uniques:
        sumPartial += value*counts[index]    
        if sumPartial>=sumTot:
            DNm = value
            break
        index+=1 
    return DNm

Documentation

The difference between my result and the result from QGIS plugin is almost 0 but still not the same (+-3 difference pixel value)