利用 Landsat-8 数据进行地表温度反演

初始化环境

In [1]:
import aie

aie.Authenticate()
aie.Initialize()

Landsat-8 数据检索

指定区域、时间、云量检索 Landsat-8 ，并对数据进行去云处理。

In [2]:
region = aie.FeatureCollection('China_Province') \
            .filter(aie.Filter.eq('province', '上海市')) 
            
dataset = aie.ImageCollection('LANDSAT_LC08_C02_T1_L2') \
             .filterBounds(region) \
             .filterDate('2019-06-01', '2019-08-31') \
             .filter(aie.Filter.lte('eo:cloud_cover', 10.0))

print(dataset.size().getInfo())

image = dataset.median()   
map = aie.Map(
    center=image.getCenter(),
    height=800,
    zoom=7
)
rgb_params = {
    'bands': ['SR_B4', 'SR_B3', 'SR_B2'],
    'min': 8000,
    'max': 13000
}
map.addLayer(
    image,
    rgb_params,
    'raw_img',
    bounds = image.getBounds()
)
map
Out[2]:

计算 NDVI

In [3]:
ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename(['NDVI'])
ndvi_params = {
    'min': -1.0,
    'max': 1.0,
    'palette': [
        '#FFFFFF', '#CE7E45', '#DF923D', '#F1B555', '#FCD163', '#99B718', '#74A901',
        '#66A000', '#529400', '#3E8601', '#207401', '#056201', '#004C00', '#023B01',
        '#012E01', '#011D01', '#011301'
    ]
}

map.addLayer(
    ndvi,
    ndvi_params,
    'NDVI',
    bounds = image.getBounds()
)
map

计算 Fractional Vegetation

In [4]:
min = aie.Image(-0.38)
max = aie.Image(0.68)
fv = (ndvi.subtract(min).divide(max.subtract(min))).pow(aie.Image(2)).rename(['FV'])
In [5]:
fv_params = {
    'min': 0,
    'max': 1
}
map.addLayer(
    fv,
    fv_params,
    'fv',
    bounds = image.getBounds()
)
map

计算 Thermal

In [6]:
th = image.select('ST_B10').multiply(aie.Image(0.00341802)).add(aie.Image(149)).rename(['TH'])
In [7]:
th_params = {
    'min': 260,
    'max': 320,
    'palette': [
        '#0000FF', '#FFFFFF', '#008000'
    ]
}
map.addLayer(
    th,
    th_params,
    'thermal',
    bounds = image.getBounds()
)
map

计算比辐射率 Emissivity

In [8]:
a = aie.Image(0.004)
b = aie.Image(0.986)
em = fv.multiply(a).add(b).rename(['EM'])
In [9]:
em_params = {
    'min': 0.9865619146722164,
    'max': 0.989699971371314
}
map.addLayer(
    em,
    em_params,
    'emissivity',
    bounds = image.getBounds()
)
map

计算地表温度 ( LST )

In [10]:
tb = th.select(['TH'])
eb = em.select(['EM'])
lst = tb.divide(tb.multiply(aie.Image(0.00115)).divide(aie.Image(1.4388)).multiply(eb.log()).add(aie.Image(1))).rename(['LST'])
In [11]:
lst_params = {
    'min': 291,
    'max': 330,
    'palette': ['#040274', '#040281', '#0502a3', '#0502b8', '#0502ce', '#0502e6',
                '#0602ff', '#235cb1', '#307ef3', '#269db1', '#30c8e2', '#32d3ef',
                '#3be285', '#3ff38f', '#86e26f', '#3ae237', '#b5e22e', '#d6e21f',
                '#fff705', '#ffd611', '#ffb613', '#ff8b13', '#ff6e08', '#ff500d',
                '#ff0000', '#de0101', '#c21301', '#a71001', '#911003']
}

map.addLayer(
    lst,
    lst_params,
    'lst',
    bounds = image.getBounds()
)
map