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:

enter image description here

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.

1 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
    #uniques values are already sorted
    for value in uniques:
        sumPartial += value*counts[index]    
        if sumPartial>=sumTot:
            DNm = value
    return DNm


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.