GEE python——利用Landsat 8 卫星进行土地分类案例
Landsat 8 图像的分类示例
在联合国实施的 REDD 计划中,在估算全国范围内的森林面积时主要依据的是遥感技术获取的土地覆盖数据。考虑到墨西哥作为幅员辽阔的国家的特点,在这一过程中仅凭自动化图像分类技术便能够以标准化且具成本效益的方式及时提供所需信息。本文介绍并详细探讨了墨西哥实用的土地覆被监测系统(Landsat Land Cover Monitoring System)。该系统通过整合国家尺度地图基准数据、所有可用 Landsat 卫星图像以及实地调查数据进行了内部验证工作。研究发现该系统自 1993 年至 2008 年间每年都能生成一张年度全国土地覆被图共计七幅图谱,并采用层次化分类方法将地形地貌划分为 9 个类别及 12 个子类别。经过评估后发现总体准确率达到 76% ,其中热带雨林区域的表现尤为突出其分类准确率高达 78% 而温带森林地区的准确率则更高达 82% 。尽管该监测系统是专为解决墨西哥 REDD+ 计划中的具体需求而设计但在遵循通用流程的基础上也可以推广至其他参与国以确保 REDD+ 计划中的方法学标准得到遵守并实现结果的一致性与可比性

该研究工作源自2014年发表的相关研究论文"MAD-MEX:使用所有 Landsat 数据对墨西哥 REDD-MRV 计划进行自动化墙到墙土地覆盖监测"的Landsat8遥感影像分类示例。如需进一步了解,请参考相关研究论文链接:

Folium包
这是一个用于生成交互式地图的Python库。它以Leayter.js为基础开发而成,并能够将地理数据以直观的形式展示并成功嵌入至Web应用程序中。Folium提供了一个完整功能模块集合来生成多样化的地图图形包括点标记、弹窗信息、热力图、线状物、面状区域以及分组控制等元素。该库以其强大的功能和高度可定制性著称,在使用时能够帮助开发者轻松构建出专业水准的地图界面
Folium的主要特点如下:
简单易用:Folium API高度易于使用。该API使其能够迅速掌握并生成美观的地图。用户只需导入Folium,在地图中添加要素后导出为HTML文件即可。
动态交互:Folium具备动态交互功能,在地图上可配置点击事件、悬停事件等多种互动类型。人机之间得以实现互动交流,并进而显著提升了用户体验水平。
- 多种地图显示方式:Folium提供多种地图显示方式, 包括OpenStreetMap/StamenTerrain/StamenToner等多种类型. 根据需求选择合适的显示方式以呈现不同类型的地理信息.
4. 多种类型的地图功能组合:Folium提供了这些功能组合,并且它们之间还可以灵活地进行叠加或替代使用。具体来说, 它们包括标记位置, 实时弹出信息框, 基于浓度分布的热力图, 连线或轮廓线以及面状区域的可视化等多种表现形式。在实际应用中, 用户可以根据具体情况选择最适合的数据展示方式以达到预期效果
可扩展性:基于Jinja2的模板引擎构建了Folium的基础架构,在这一设计下,开发者能够通过定制化地图展示模块的开发设计来显著提升其功能性。在现有功能框架的基础上,允许开发者编写更加丰富的地图展示逻辑,并通过相应的配置参数实现特定业务需求下的个性化地图呈现效果。
6. 兼容性:Folium支持与其他主流Python工具集良好集成,包括但不限于Pandas、GeoPandas和Matplotlib。通过与这些工具集的集成,用户能够方便地处理地理数据并将其可视化在Folium生成的地图上。
使用Folium创建地图的基本步骤如下:
1. 导入Folium包:首先需要导入Folium包,以便可以使用其中的API。
2. 创建地图:借助Folium的Map类功能模块,在应用中能够生成一个空白的地图界面。根据需求配置地图的中心坐标位置和缩放比例。
通过调用add_to方法可以在地图上增添标记点、弹窗以及热力图等多种地图要素,并可以根据需要设置这些要素的位置、样式以及属性等信息。
的地图显示:通过使用save方法将地图保存为HTML格式后,通过浏览器打开即可查看地图.
在基础的地图创建和要素配置之外,Folium提供了一系列高级功能组件。包括分组控制台、图例面板以及 click 和 hover 事件处理接口等。这些高级功能将显著提升地图的交互体验与视觉效果。
作为一个强大的Python库,Folium被广泛应用于地图交互与地理数据可视化领域;它不仅提供了丰富的功能,还具备高度的可定制性,能够帮助用户轻松构建专业级的地图.无论是从事地理信息科学、数据可视化还是Web开发,这一工具都能满足您的需求.对于初学者而言,Folium凭借其直观易懂的学习曲线,能让您快速掌握地图制作技能;而对于资深开发者来说,其强大的API体系则能为企业级项目提供有力支持.
代码
from IPython.display import Image
import ee, folium
ee.Initialize()
%matplotlib inline
获取 Landsat 8 输入数据并将其可视化
area_of_interest = ee.Geometry.Rectangle([-98.75, 19.15, -98.15,18.75])
mexico_landcover_2010_landsat = ee.Image("users/renekope/MEX_LC_2010_Landsat_v43").clip(area_of_interest)
landsat8_collection = ee.ImageCollection('LANDSAT/LC8_L1T_TOA').filterDate('2016-01-01', '2018-04-19').min()
landsat8_collection = landsat8_collection.slice(0,9)
vis = {
'bands': ['B6', 'B5', 'B2'],
'min': 0,
'max': 0.5,
'gamma': [0.95, 1.1, 1],
'region':area_of_interest}
image = landsat8_collection.clip(area_of_interest)
mapid = image.getMapId(vis)
map = folium.Map(location=[19.15,-98.75],zoom_start=9, height=500,width=700)
folium.TileLayer(
tiles=mapid['tile_fetcher'].url_format,
attr='Map Data © <a href="https://earthengine.google.com/">Google Earth Engine</a>',
overlay=True,
name='Landsat 8 ',
).add_to(map)
map.add_child(folium.LayerControl())
map
推导植被指数和其他栅格操作的功能
def NDVI(image):
return image.normalizedDifference(['B5', 'B4'])
def SAM(image):
band1 = image.select("B1")
bandn = image.select("B2","B3","B4","B5","B6","B7","B8","B9");
maxObjSize = 256;
b = band1.divide(bandn);
spectralAngleMap = b.atan();
spectralAngleMap_sin = spectralAngleMap.sin();
spectralAngleMap_cos = spectralAngleMap.cos();
sum_cos = spectralAngleMap_cos.reduce(ee.call("Reducer.sum"));
sum_sin = spectralAngleMap_sin.reduce(ee.call("Reducer.sum"));
return ee.Image.cat(sum_sin, sum_cos, spectralAngleMap_sin, spectralAngleMap_cos);
#Enhanced Vegetation Index
def EVI(image):
# L(Canopy background)
# C1,C2(Coefficients of aerosol resistance term)
# GainFactor(Gain or scaling factor)
gain_factor = ee.Image(2.5);
coefficient_1 = ee.Image(6);
coefficient_2 = ee.Image(7.5);
l = ee.Image(1);
nir = image.select("B5");
red = image.select("B4");
blue = image.select("B2");
evi = image.expression(
"Gain_Factor*((NIR-RED)/(NIR+C1*RED-C2*BLUE+L))",
{
"Gain_Factor":gain_factor,
"NIR":nir,
"RED":red,
"C1":coefficient_1,
"C2":coefficient_2,
"BLUE":blue,
"L":l
}
)
return evi
#Atmospherically Resistant Vegetation Index
def ARVI(image):
red = image.select("B4")
blue = image.select("B2")
nir = image.select("B5")
red_square = red.multiply(red)
arvi = image.expression(
"NIR - (REDsq - BLUE)/(NIR+(REDsq-BLUE))",{
"NIR": nir,
"REDsq": red_square,
"BLUE": blue
}
)
return arvi
#Leaf Area Index
def LAI(image):
nir = image.select("B5")
red = image.select("B4")
coeff1 = ee.Image(0.0305);
coeff2 = ee.Image(1.2640);
lai = image.expression(
"(((NIR/RED)*COEFF1)+COEFF2)",
{
"NIR":nir,
"RED":red,
"COEFF1":coeff1,
"COEFF2":coeff2
}
)
return lai
def tasseled_cap_transformation(image):
#Tasseled Cap Transformation for Landsat 8 based on the
#scientfic work "Derivation of a tasselled cap transformation based on Landsat 8 at-satellite reflectance"
#by M.Baigab, L.Zhang, T.Shuai & Q.Tong (2014). The bands of the output image are the brightness index,
#greenness index and wetness index.
b = image.select("B2", "B3", "B4", "B5", "B6", "B7");
#Coefficients are only for Landsat 8 TOA
brightness_coefficents= ee.Image([0.3029, 0.2786, 0.4733, 0.5599, 0.508, 0.1872])
greenness_coefficents= ee.Image([-0.2941, -0.243, -0.5424, 0.7276, 0.0713, -0.1608]);
wetness_coefficents= ee.Image([0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559]);
fourth_coefficents= ee.Image([-0.8239, 0.0849, 0.4396, -0.058, 0.2013, -0.2773]);
fifth_coefficents= ee.Image([-0.3294, 0.0557, 0.1056, 0.1855, -0.4349, 0.8085]);
sixth_coefficents= ee.Image([0.1079, -0.9023, 0.4119, 0.0575, -0.0259, 0.0252]);
#Calculate tasseled cap transformation
brightness = image.expression(
'(B * BRIGHTNESS)',
{
'B':b,
'BRIGHTNESS': brightness_coefficents
})
greenness = image.expression(
'(B * GREENNESS)',
{
'B':b,
'GREENNESS': greenness_coefficents
})
wetness = image.expression(
'(B * WETNESS)',
{
'B':b,
'WETNESS': wetness_coefficents
})
fourth = image.expression(
'(B * FOURTH)',
{
'B':b,
'FOURTH': fourth_coefficents
})
fifth = image.expression(
'(B * FIFTH)',
{
'B':b,
'FIFTH': fifth_coefficents
})
sixth = image.expression(
'(B * SIXTH)',
{
'B':b,
'SIXTH': sixth_coefficents
})
bright = brightness.reduce(ee.call("Reducer.sum"));
green = greenness.reduce(ee.call("Reducer.sum"));
wet = wetness.reduce(ee.call("Reducer.sum"));
four = fourth.reduce(ee.call("Reducer.sum"));
five = fifth.reduce(ee.call("Reducer.sum"));
six = sixth.reduce(ee.call("Reducer.sum"));
tasseled_cap = ee.Image(bright).addBands(green).addBands(wet).addBands(four).addBands(five).addBands(six)
return tasseled_cap.rename('brightness','greenness','wetness','fourth','fifth','sixth')
进行缨帽变换
tct = tasseled_cap_transformation(landsat8_collection)
image = tct.clip(area_of_interest)
vis_tct = {'min':-1,'max':2,'size':'800',
'bands':['brightness','greenness','wetness'],
'region':area_of_interest}
mapid = image.getMapId(vis_tct)
map = folium.Map(location=[19.15,-98.75],zoom_start=9, height=500,width=700)
folium.TileLayer(
tiles=mapid['tile_fetcher'].url_format,
attr='Map Data © <a href="https://earthengine.google.com/">Google Earth Engine</a>',
overlay=True,
name='Tasseled Cap Transformation',
).add_to(map)
map.add_child(folium.LayerControl())
map
推导指数和光谱角 建立图像堆栈并将其可视化
ndvi = NDVI(landsat8_collection)
sam = SAM(landsat8_collection)
evi = EVI(landsat8_collection)
arvi = ARVI(landsat8_collection)
lai = LAI(landsat8_collection)
spectral_indices_stack = ee.Image(ndvi).addBands(lai).addBands(sam).addBands(arvi).addBands(evi).addBands(tct).addBands(landsat8_collection)
image = ndvi.clip(area_of_interest)
vis_ndvi = {'min':-1,'max':1,'size':'800',
'region':area_of_interest}
mapid = image.getMapId(vis_ndvi)
map = folium.Map(location=[19.15,-98.75],zoom_start=9, height=500,width=700)
folium.TileLayer(
tiles=mapid['tile_fetcher'].url_format,
attr='Map Data © <a href="https://earthengine.google.com/">Google Earth Engine</a>',
overlay=True,
name='NDVI',
).add_to(map)
map.add_child(folium.LayerControl())
map
定义分类功能
def classification(raster_input, training_dataset,number_of_training_points, region, classification_algorithm):
bands = raster_input.bandNames()
points = ee.FeatureCollection.randomPoints(region, number_of_training_points, number_of_training_points, 1)
training = training_dataset.addBands(raster_input).reduceToVectors(
reducer='mean',
geometry=points,
geometryType='centroid',
scale=30,
crs='EPSG:4326'
)
classifier = ee.Classifier.randomForest().train(
features=training,
classProperty='label',
inputProperties=raster_input.bandNames(),
)
out = raster_input.classify(classifier)
return out
该文主要对基于函数的分类方法进行系统阐述,并对多种分类模型在不同应用场景下的应用方式进行深入探讨。具体而言,在讨论过程中首先阐述了几类典型模型的基本原理;然后深入探讨其在实际应用中所涉及的具体实现方式,并通过理论分析与案例研究相结合的方式展开推导过程
output = classification(spectral_indices_stack, mexico_landcover_2010_landsat, 10000, area_of_interest, 'Cart')
palette = ['5d9cd4','007e00','003c00','aaaa00','aa8000','8baa00','ffb265','00d900','aa007f','ff55ff','ff557f','ff007f','ff55ff','aaffff','00ffff','55aaff','e29700','bd7e00','966400','a2ecb1','c46200','aa5500','6d3600','00aa7f','008a65','005941','e9e9af','faff98',
'00007f','c7c8bc','4d1009','000000','fef7ff','6daa50','3a7500','0b5923','ffaaff','ffd1fa']
palette = ','.join(palette)
# make a visualizing variable
vis_classification = {'min': 0, 'max': len(palette), 'palette': palette, 'region':area_of_interest}
显示分类训练数据
image = mexico_landcover_2010_landsat.clip(area_of_interest)
mapid = image.getMapId(vis_classification)
map = folium.Map(location=[19.15,-98.75],zoom_start=9, height=500,width=700)
folium.TileLayer(
tiles=mapid['tile_fetcher'].url_format,
attr='Map Data © <a href="https://earthengine.google.com/">Google Earth Engine</a>',
overlay=True,
name='Training Data',
).add_to(map)
map.add_child(folium.LayerControl())
map
显示分类结果
请耐心等待。可能需要一些时间。您可能需要多次运行此单元。
image = output.clip(area_of_interest)
mapid = image.getMapId(vis_classification)
map = folium.Map(location=[19.15,-98.75],zoom_start=9, height=500,width=700)
folium.TileLayer(
tiles=mapid['tile_fetcher'].url_format,
attr='Map Data © <a href="https://earthengine.google.com/">Google Earth Engine</a>',
overlay=True,
name='Classification Output',
).add_to(map)
map.add_child(folium.LayerControl())
map

