Taking 4 variables as an input (shapefile of our region of interest, basename of the ROI for Landsat image (e.g.LC8009061), year, first doy (day of year) and number of images to analyze) this program creates a result layer in which the frequency of cloud appearance for each pixel is shown. It downloads the Landsat Quality band files (if not pres…

clouds.py

import os
import subprocess
import urllib2
import sys
import numpy as np
from osgeo import gdal

GDAL_OPTS = ["COMPRESS=LZW","INTERLEAVE=PIXEL", "TILED=YES",\
        "SPARSE_OK=TRUE", "BIGTIFF=YES" ]
"""Taking 4 variables as an input (shapefile of our region of interest, basename of the ROI for Landsat image (e.g.LC8009061), year, first doy (day of year)
and number of images to analyze) this program creates a result layer in which the frequency of cloud appearance for each pixel is shown. 
It downloads the Landsat Quality band files (if not present), clips them using the shapefile, analyzes cloudiness using a decimal-binary converter, 
and generates the result layer.
Autor: dizlois@gmail.com
"""
        
def fetch_image(name,doy,year,band=None,ground=None,version=None):
    """ Landsat filename convention: e.g.LC80390222000076EDC00
    L: Landsat
    C/O/T: Instrument; C=Combined, O=OLI, T=TIRS
    8: Satellite
    039:	Path (WRS-2)
    022:	Row (WRS-2)
    2000:	Year
    076:	Julian Day (Day 76 of the calendar year = 17 March)
    EDC:	Ground Station where the data was received
    00:	        Archive version number
    """

    if not band: band="BQA" 
    if not ground: ground="LGN"
    if not version: version="00"
    

    ###Name= Mission and region e.g. LC8010060
    ###Download file and save it in the current directory
        
    rootname=name+str(year)+str(doy)+ground+str(version)
    
    url="http://landsat-pds.s3.amazonaws.com/"+name[0]+name[2]+"/"+name[3:6]+"/"+name[6:9]+"/"+rootname+"/"+rootname+"_"+band+".TIF"
    file_name=url.split('/')[-1]
    if os.path.exists ( file_name ):
        print "File found in current directory"
        return file_name
    fetch = urllib2.urlopen(url)
    output = open(file_name,'wb')

    ###Progress bar adapted from : http://blog.radevic.com/2012/07/python-download-url-to-file-with.html
    meta = fetch.info()
    file_size = int(meta.getheaders("Content-Length")[0])
    print("Downloading: {0} Bytes: {1}".format(url, file_size))
    file_size_dl = 0
    block_sz = 8192
    while True:
        buffer = fetch.read(block_sz)
        if not buffer:
            break
        file_size_dl += len(buffer)
        output.write(buffer)
        p = float(file_size_dl) / file_size
        status = "\r{0}  [{1:.2%}]".format(file_size_dl, p)
        status = status + chr(8)*(len(status)+1)
        sys.stdout.write(status) 
    output.write(fetch.read())
    output.close() 
    return file_name
    
##Clipping rasters....       
def cut_shape (filename,shapefile):
    print "Cutting %s raster..."%(filename)
    destination_file=filename.replace ( ".TIF", "_clipped.tif" )
    if os.path.exists ( destination_file ):
        print "Clipped file already in disk..."
        return destination_file
        
    retval = subprocess.call (["gdalwarp","-co","compress=lzw","-srcnodata","1","-dstnodata","-9999","-cutline",shapefile,filename,destination_file], shell=True )
    if retval != 0:
        raise IOError ("Problem selecting input files... ")
    print "Done"
    return destination_file

def Binary_from_decimal(number,bitcount=None):
    ###bitcount variable, when used, set the result to a specific number of bits (larger than minimum)
    Binary=[]
    while number>1:
        if number%2==0:
            Binary.insert(0,0)
        if number%2==1:
            Binary.insert(0,1)
        number=number/2
    Binary.insert(0,1)
    if bitcount:
        dif=bitcount-len(Binary)
        if dif!=0:
            z=dif*[0]
            z.extend(Binary)
            Binary=z
    return Binary

def process_data(clipped_image,results_array,Nodatavalue=None):
    if not Nodatavalue:
        Nodatavalue=-9999
    g=gdal.Open(clipped_image)
    g_array=g.GetRasterBand(1).ReadAsArray()
    tracer=0
    for i in xrange(g_array.shape[0]):
        for j in xrange(g_array.shape[1]):
            if g_array[i,j] !=Nodatavalue:
                if Binary_from_decimal(g_array[i,j],16)[1]+Binary_from_decimal(g_array[i,j],16)[0]==2:
                    results_array[i,j]+=1 
                    tracer+=1
    
    print tracer, " clouded pixels detected"
    return results_array
                    


def get_results (working_directory,shapefile,base_name,first_doy,year,number_images):
    
    first_time=True
    for i in xrange(number_images):
        
        doy=first_doy+(i*16)
        if doy<10:
            doy="00"+str(doy)
        if (10 < doy and doy <100):
            doy="0"+str(doy)
        print "doy ->",doy      
        file_name=working_directory+fetch_image(base_name,doy,year)
        clipped_file=cut_shape (file_name,shapefile)
        g=gdal.Open(clipped_file)
        array=g.GetRasterBand(1).ReadAsArray()
        if first_time:
            results=array*0.
            results[np.where(array==-9999)]=-9999
            first_time=False
        if i+1<number_images:
            results=process_data(clipped_file,results)
            
        if i+1==number_images:
            
            results=(process_data(clipped_file,results))/(i+1)
            drv = gdal.GetDriverByName ( "GTiff" )
            dst_ds = drv.Create ( "results.TIF", g.RasterXSize,\
                g.RasterYSize, 1, gdal.GDT_Float32,options=GDAL_OPTS )
            dst_ds.SetGeoTransform ( g.GetGeoTransform() )
            dst_ds.SetProjection ( g.GetProjectionRef() )
            dst_ds.GetRasterBand(1).WriteArray ( results)
            dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
            dst_ds.GetRasterBand(1).WriteArray ( results)
            dst_ds = None
#Change paths accordingly...
working_directory= "C:/Users/usuario/Desktop/Varios/Scripts/"
os.chdir(working_directory)
shapefile="C:/Users/usuario/Desktop/Varios/Fundacion matrix/RBCC/COLONSO_CHALUPAS.shp"

base_name= "LC8009061"
year=2015
first_doy=11
number_of_images=23

get_results (working_directory,shapefile,base_name,first_doy,year,number_of_images)