cloudsen12v2 qgis viz

qgismap.py

import matplotlib.pyplot as plt
from ee_plugin import Map
import numpy as np
import ee
import re

try:
    ee.Image(0)
except:
    ee.Initialize()


def get_percentile(pname = "p05"):
    # split viz
    p1 = float(re.split(r"(\d+)", pname)[1])
    try:
        p2 = float(re.split(r"(\d+)", pname)[3])
    except:
        p2 = p1
    return p1, p2


def getviz(id, rgb = "RGB", viz = "minmax"):
    # S2 copernicus
    s2id_ee = ee.Image("COPERNICUS/S2/%s" % id)    
    
    # Get histograms
    if rgb == "SNOW":
        S11hist  = s2id_ee.select("B11").reduceRegion(ee.Reducer.histogram(), s2id_ee.geometry(), scale = 2000).getInfo()
                
        if viz == "minmax":
            Bmin, Bmax = np.min(S11hist["B11"]["bucketMeans"]), np.max(S11hist["B11"]["bucketMeans"])
        elif viz[0]  == "p":
            p1, p2 = get_percentile(viz)
                        
            histogram_b11 = S11hist["B11"]["bucketMeans"]
            Bmin, Bmax = histogram_b11[int(p1)], histogram_b11[len(histogram_b11) - int(p2) - 1]
        else:
            raise ValueError("viz must be minmax or pxx")
        print("SNOW: [min: %s, max: %s]" % (round(Bmin/10000, 3), round(Bmax/10000, 3)))
        Map.addLayer(s2id_ee, {'bands': ['B12', 'B11', 'B1'], 'min': Bmin, 'max': Bmax}, 'SNOW', True)

    elif rgb == "RGB":
        
        S4hist  = s2id_ee.select("B4").reduceRegion(ee.Reducer.histogram(), s2id_ee.geometry(), scale = 2000).getInfo()
        
        if viz == "minmax":
            Bmin, Bmax = np.min(S4hist["B4"]["bucketMeans"]), np.max(S4hist["B4"]["bucketMeans"])
        elif viz[0]  == "p":
            p1, p2 = get_percentile(viz)            
            histogram_b4 = S4hist["B4"]["bucketMeans"]
            Bmin, Bmax = histogram_b4[int(p1)], histogram_b4[len(histogram_b4) - int(p2) - 1]
        else:
            raise ValueError("viz must be minmax or pxx")
        print("RGB: [min: %s, max: %s]" % (round(Bmin/10000, 3), round(Bmax/10000, 3)))
        Map.addLayer(s2id_ee, {'bands': ['B4', 'B3', 'B2'], 'min': Bmin, 'max': Bmax}, 'RGB', True)
                
    elif rgb == "CIRRUS":
        cirrus = s2id_ee.select("B10").divide(10000)
        S10hist  = cirrus.reduceRegion(ee.Reducer.histogram(), s2id_ee.geometry(), scale = 2000).getInfo()
        if viz == "minmax":
            B10min, B10max = np.min(S10hist["B10"]["bucketMeans"]), np.max(S10hist["B10"]["bucketMeans"])
        elif viz[0]  == "p":
            p1, p2 = get_percentile(viz)
            histogram_b10 = S10hist["B10"]["bucketMeans"]
            B10min, B10max = histogram_b10[int(p1)], histogram_b10[len(histogram_b10) - int(p2) - 1]
        else:
            raise ValueError("viz must be minmax or pxx")
        print("CIRRUS: [min: %s, max: %s]" % (round(B10min, 3), round(B10max, 3)))
        Map.addLayer(cirrus, {'palette': ["#040ed8", "#86d9ff", "#fffe47", "#ff4800", "#9e0000"], 'min': B10min, 'max': B10max}, 'CIRRUS', True)

    return None



# GLOBAL PARAMS -----
id = "20210929T181039_20210929T181920_T12TWP"
getviz(id, rgb = "RGB", viz = "p00-95") # p00-95
getviz(id, rgb = "SNOW", viz = "p00-95") # p05-05
getviz(id, rgb = "CIRRUS", viz = "p00-00") # p05-05