Advertisement

GEE python —— Landsat 8图像的分类实例

阅读量：

Classification example for Landsat 8 imagery

土地分类作为地学科学的重要领域之一，在Google Earth Engine平台中应用广泛。它不仅有助于我们更深入地认识和有效地管理地球表面的各类土地类型，并且能够利用其强大的遥感数据处理与分析功能实现精准的土地划分与识别。具体而言，在GEE平台上完成土地分类的基本流程包括以下几个步骤：首先获取并导入所需范围内的Landsat影像数据；其次对原始影像图像进行必要的预处理工作；然后运用监督学习算法对影像图像进行系统化的分类；最后对分类结果进行精确评估并生成相应的地理信息系统（GIS）地图。

  1. 获取遥感数据：在Google Earth Engine平台内可访问并利用这些卫星影像的数据集。其中包含常见的Landsat和Sentinel系列卫星影像。此外还可以上传自已的遥感影像数据集进行分析和研究。

  2. 数据预处理：完成遥感数据的预处理工作。具体而言，在该过程中需要完成包括云覆盖抑制（云去除）、大气补偿（大气校正）、辐射平衡校正（辐射校正）、图像裁剪以及几何投影（投影）等多个关键步骤。

3.特征提取：使用遥感数据来提取地物的特征，如植被指数、地表温度等。

  1. 训练分类器：采用机器学习技术进行分类器训练，在GEE平台中支持应用Supervised Classification和Random Forest等算法。

5.分类结果验证：使用验证样本集来验证分类结果的准确性和可靠性。

6.分类结果输出：生成并展示分类结果的空间数据表示形式，并完成后续的数据处理与分析。

此乃基础的土地分类流程，在GEE平台实施土地分类需运用遥感技术和机器学习方法。建议开展操作前深入学习相关技术原理并实际操作训练。

请参考以下关于GEE JavaScript分类的文章：

以及其中的土地分类相关内容：

土地分类案例

土地分类结果

This classification case study is derived from a scientific study investigating MAD-MEX, an automated system designed for wall-to-wall land cover monitoring within the Mexican REDD-MRV program, which utilizes all Landsat data, as demonstrated by GEBHARDT et al. (2014). For further details, please visit this link: https://www.mdpi.com/2072-4292/6/5/3923.

复制代码 
    from IPython.display import Image
    import ee, folium
    ee.Initialize()
    %matplotlib inline

Get and visualize the Landsat 8 input data
获取指定区域的Landsat影像

复制代码 
    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 &copy; <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')

Through derivation and subsequent graphical representation, the Tasseled Cap Transformation can be effectively visualized.

复制代码 
    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 &copy; <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

Derive indices, spectral angles. Build and visualize image stack

复制代码 
    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 &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='NDVI',
      ).add_to(map)
    
    map.add_child(folium.LayerControl())
    map

Define classification function
定义分类函数

复制代码 
    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

Derive classification function
推导出分类函数

复制代码 
    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}

Display training data of classification
展示分类后的结果

复制代码 
    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 &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='Training Data',
      ).add_to(map)
    
    map.add_child(folium.LayerControl())
    map

Display classification output
Please be patient. It may take a few moments. You might have to run this cell several times.
显示分类输出
请耐心等待。它可能需要一些时间。你可能要运行这行代码数次。

复制代码 
    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 &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    overlay=True,
    name='Classification Output',
      ).add_to(map)
    
    map.add_child(folium.LayerControl())
    map

全部评论 (0)

还没有任何评论哟~

相关文章推荐

GEE python —— Landsat 8图像的分类实例

ClassificationexampleforLandsat8imagery 土地分类是地学科学中的一个重要领域，可以帮助我们更好地了解和管理地球表面的不同土地类型。

基于GEE使用Landsat 8和Landsat 5影像进行分类

2018年Landsat8影像分类 //数据集筛选 vardataset1=ee.ImageCollection'LANDSAT/LC08/C02/T1L2' .filterDate'2018xx',...

GEE python——利用Landsat 8 卫星进行土地分类案例

Landsat8图像的分类示例 在联合国降低因森林砍伐和退化所产生的排放（REDD）计划中，估算全国范围内的森林面积主要基于使用遥感技术的土地覆盖信息。对于墨西哥这样一个幅员辽阔的国家来说，只有通过自...

基于GEE使用Landsat 8和Landsat 5影像计算RSEI

一、使用2018年landsat8影像数据计算遥感生态指数RSEI //couldmasking，GEE上有现有的去云函数 functionmaskL8srimage //Bit0Fill //Bit...

Google Earth Engine（GEE）—Julia Landsat 8 去云案例

代码 usingEarthEngine Initialize 这里的去云波段和之前的sr去云没有什么区别 FunctiontocloudmaskfromthepixelqabandofLandsat8...

GEE案例——Landsat系列（4、5、7、8和9）影像的物候监测分析

简介 Landsat物候分析是一种利用Landsat系列卫星数据进行植被生长和变化的研究方法。物候是指植物生长和发育的季节性变化，包括植物生长、开花、结果和凋落等过程。 Landsat系列卫星是美国国...

GEE Python：利用Python API实现Sentinel-1、2和Landsat-8、9的可用影像数量

目录 简介 参数 Python代码 Python内置函数 简介 生成以下卫星图像集合的可用性图表：1.哨兵12.哨兵23.Landsat84.Landsat9 参数 startDate（开始日期图像集...

GEE 利用Landsat-8影像计算年内最大ndvi

1\.导入内蒙古自治区矢量边界 直接上传边界到assets,import 2\.编译代码 varl8=ee.ImageCollectionLANDSAT/LC08/C01/T1RTTOA; //计算2...

GEE计算遥感生态指数(RSEI)--Landsat 8为例

目录 RSEI原理 湿度指标（Wet） 绿度指标（NDVI） 热度指标（LST） 干度指标（NDBSI） Landsat8波段 归一化（normalization） 主成分分析PCA PCARSEI ...

GEE土地分类——利用landsat 8 和随机森林方法进行土地分类

目录 简介 代码解释 代码 函数 ee.Classifier.smileRandomForestnumberOfTrees,variablesPerSplit,minLeafPopulation,ba...