基于 Landsat-8 數(shù)據(jù)進(jìn)行貴陽市冷熱島分析
初始化環(huán)境
import aie
aie.Authenticate()
aie.Initialize()
Landsat-8 數(shù)據(jù)檢索
指定區(qū)域熊响、時間旨别、云量檢索 Landsat-8 ，并對數(shù)據(jù)進(jìn)行去云處理汗茄。
region = aie.FeatureCollection('China_City') \
.filter(aie.Filter.eq('city', '貴陽市')) \
.geometry()
dataset = aie.ImageCollection('LANDSAT_LC08_C02_T1_L2') \
.filterBounds(region) \
.filterDate('2021-05-01', '2021-10-1') \
.filter(aie.Filter.lte('eo:cloud_cover', 30.0))
print(dataset.size().getInfo())
Landsat-8 數(shù)據(jù)去云
def removeLandsatCloud(image):
cloudShadowBitMask = (1 << 4)
cloudsBitMask = (1 << 3)
qa = image.select('QA_PIXEL')
mask = qa.bitwiseAnd(aie.Image(cloudShadowBitMask)).eq(aie.Image(0)).And(qa.bitwiseAnd(aie.Image(cloudsBitMask)).eq(aie.Image(0)))
return image.updateMask(mask)
# median運(yùn)算
image = dataset.median()
# clip裁剪
image = image.clip(region)
# 去云
# 去云會造成數(shù)據(jù)缺失，所以我這里不去云了
# image = removeLandsatCloud(image)
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
計算 NDVI
ndvi計算 NIR-R/NIR+R
ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename(['NDVI'])
ndvi_params = {
'min': 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
貴陽市ndvi分布圖
task = aie.Export.image.toAsset(ndvi,'guiyang_ndvi',50)
task.start()
計算 Fractional Vegetation
植被覆蓋度計算 NDVI-NDVI_min / NDVI_max-NDVI_min
import numpy as np
scale = 1000
histogram = ndvi.reduceRegion(aie.Reducer.histogram(2000), None, scale)
histogram_info = histogram.getInfo()
# print(histogram_info)
bucketKey = histogram_info['NDVI_range']
bucketValue = histogram_info['NDVI_counts']
key = np.array(bucketValue)
accSum = np.cumsum(key)
# print(accSum[20])
# print(accSum[-1])
accPercent = accSum / accSum[-1]
p2 = np.searchsorted(accPercent, 0.2)
min_ndvi = bucketKey[p2 + 1]
print('min_ndvi:%f' % min_ndvi)
p98 = np.searchsorted(accPercent, 0.98)
max_ndvi = bucketKey[p98]
print('max_ndvi:%f' % max_ndvi)
ndvi
min = aie.Image(min_ndvi)
max = aie.Image(max_ndvi)
fv = ndvi.where(ndvi.lt(min),aie.Image(0))\
.where(ndvi.gt(max),aie.Image(1))\
.where(ndvi.gt(min).And(ndvi.lt(max)),\
ndvi.subtract(min).divide(max.subtract(min)))\
.rename(['FV'])
fv_params = {
'min': 0,
'max': 1
}
map.addLayer(
fv,
fv_params,
'fv',
bounds = image.getBounds()
)
map
計算比輻射率 Emissivity
em = ndvi.where(ndvi.lte(aie.Image(0.2)),aie.Image(0.995))\
.where(ndvi.gt(aie.Image(0.2)).And(ndvi.lt(aie.Image(0.4))),\
aie.Image(0.9589).add((aie.Image(0.086)).multiply(fv))\
.subtract(aie.Image(0.0617).multiply(fv).pow(aie.Image(2))))\
.where(ndvi.gte(aie.Image(0.4)),aie.Image(0.9625).add(aie.Image(0.0614))\
.multiply(fv).subtract(aie.Image(0.0461)).multiply(fv).pow(aie.Image(2)))\
.rename(['EM'])
scale = 1000
histogram = em.reduceRegion(aie.Reducer.mean())
histogram_info = histogram.getInfo()
print(histogram_info)
histogram = em.reduceRegion(aie.Reducer.max())
histogram_info = histogram.getInfo()
print(histogram_info)
histogram = em.reduceRegion(aie.Reducer.min())
histogram_info = histogram.getInfo()
print(histogram_info)
EM
計算 Thermal
同溫黑體輻射亮度計算
th = image.select('ST_B10').multiply(aie.Image(0.00341802)).add(aie.Image(149)).rename(['TH'])
scale = 1000
histogram = th.reduceRegion(aie.Reducer.mean())
histogram_info = histogram.getInfo()
print(histogram_info)
histogram = th.reduceRegion(aie.Reducer.max())
histogram_info = histogram.getInfo()
print(histogram_info)
histogram = th.reduceRegion(aie.Reducer.min())
histogram_info = histogram.getInfo()
print(histogram_info)
Thermal范圍
th_params = {
'min': 240,
'max': 330,
'palette': [
'#0000FF', '#FFFFFF', '#008000'
]
}
map.addLayer(
th,
th_params,
'thermal',
bounds = image.getBounds()
)
map
計算地表溫度 ( LST )
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))).subtract(aie.Image(273)).rename(['LST'])
scale = 1000
histogram = lst.reduceRegion(aie.Reducer.mean())
histogram_info = histogram.getInfo()
print(histogram_info)
histogram = lst.reduceRegion(aie.Reducer.max())
histogram_info = histogram.getInfo()
print(histogram_info)
histogram = lst.reduceRegion(aie.Reducer.min())
histogram_info = histogram.getInfo()
print(histogram_info)
histogram = lst.reduceRegion(aie.Reducer.histogram(2000), None, scale)
histogram_info = histogram.getInfo()
# print(histogram_info)
bucketKey = histogram_info['LST_range']
bucketValue = histogram_info['LST_counts']
key = np.array(bucketValue)
accSum = np.cumsum(key)
# print(accSum[20])
# print(accSum[-1])
accPercent = accSum / accSum[-1]
p1 = np.searchsorted(accPercent, 0.01)
min_lst = bucketKey[p1 + 1]
print('min_LST:%f' % min_lst)
p99 = np.searchsorted(accPercent, 0.99)
max_lst = bucketKey[p99]
print('max_LST:%f' % max_lst)
lst范圍
消除lst噪點
用1%和99%消除噪點
# 消除噪點，用1%和99.9%消除噪點
lst = lst.where(lst.lte(aie.Image(10)),aie.Image(10))\
.where(lst.gt(aie.Image(10)).And(ndvi.lt(aie.Image(44))),lst)\
.where(lst.gte(aie.Image(44)),aie.Image(44))
lst_params = {
'min': 10,
'max': 45,
'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
貴陽市地表溫度結(jié)果
task = aie.Export.image.toAsset(lst,'lst',50)
task.start()
城市冷熱島分析
城市冷島分析
scale = 1000
histogram = lst.reduceRegion(aie.Reducer.mean())
histogram_info = histogram.getInfo()
print(histogram_info)
mean_lst = aie.Image(histogram_info.get('LST_mean'))
# subtract運(yùn)算
uci = lst.subtract(mean_lst)
# lst二分類递览，將lst小于0的區(qū)域設(shè)置為1叼屠，將lst大于0的區(qū)域設(shè)置為nodata
uci_extent = lst.where(lst.lte(mean_lst),aie.Image(1))\
.where(lst.gt(mean_lst),aie.Image(0))
uci_final = lst.where(lst.lte(mean_lst),uci)\
.where(lst.gt(mean_lst),aie.Image(0))
task = aie.Export.image.toAsset(uci_final,'guiyang_uci_final',50)
task.start()
mask_vis = {
'min': 0,
'max': 1,
'palette': ['#ffffff', '#0000cd'] # 0:白色, 1:藍(lán)色
}
map.addLayer(uci_extent,mask_vis, 'uci_extent', bounds=region.getBounds()) # 藍(lán)色區(qū)域為冷島
冷島
histogram = uci_final.reduceRegion(aie.Reducer.min())
histogram_info = histogram.getInfo()
print(histogram_info)
冷島強(qiáng)度分布圖
城市熱島分析
tr = lst.subtract(mean_lst).divide(mean_lst)
# lst二分類，將lst小于0的區(qū)域設(shè)置為1绞铃，將lst大于0的區(qū)域設(shè)置為nodata
uhi_extent = tr.where(tr.lt(aie.Image(0)),aie.Image(0))\
.where(tr.gt(aie.Image(0)),aie.Image(1))
mask_vis = {
'min': 0,
'max': 1,
'palette': ['#ffffff', '#ff7070'] # 0:白色, 1:紅色
}
map.addLayer(uhi_extent,mask_vis, 'uhi_extent', bounds=region.getBounds()) # 紅色區(qū)域為熱島
熱島
uhi_final = uhi_extent.multiply(tr)
task = aie.Export.image.toAsset(tr,'guiyang_uhi_final',50)
task.start()
uhi_params = {
'min': 0.1,
'max': 0.8,
'palette': ['#fff705', '#ffd611', '#ffb613', '#ff8b13', '#ff6e08', '#ff500d',
'#ff0000', '#de0101', '#c21301', '#a71001', '#911003']
}
map.addLayer(uhi_final,uhi_params, 'uhi_final', bounds=region.getBounds())
熱島強(qiáng)度圖
總結(jié)
總的來說环鲤，這個不難，難在精準(zhǔn)憎兽，難在精準(zhǔn)的地表溫度反演冷离，比如我遇到的異常值問題，需要去除噪點纯命。還有就是阿里云的小哥哥們很熱心西剥，也多虧他們的幫助，當(dāng)然我這個地表溫度反演做的比較簡單亿汞，也不是很準(zhǔn)確瞭空，希望大家有更加準(zhǔn)確的算法。這個實驗還是比較好的擬合了貴陽市的實際疗我，也說明AIE可以很好的使用咆畏。
本文主要引用了AIE的官方案例，《城市與區(qū)域規(guī)劃空間分析實驗教程》的實驗19吴裤，還有網(wǎng)絡(luò)上使用ENVI利用landsat8數(shù)據(jù)進(jìn)行地表溫度反演的例子旧找。
