GEE python：Landsat进行土地分类（全代码）

Classification example for Landsat 8 imagery

土地分类是地学科学中的一个重要领域，可以帮助我们更好地了解和管理地球表面的不同土地类型。在Google Earth Engine（GEE）中，在使用Landsat影像进行土地分类时，进行土地分类可以利用其强大的遥感数据处理和分析能力。以下是使用GEE进行土地分类的基本步骤：

1.收集遥感数据：使用GEE中提供的遥感数据集，如Landsat、Sentinel等，或者上传自己的遥感数据集。

2.数据预处理：对遥感数据进行预处理，如云去除、大气校正、辐射校正、裁剪、投影等。

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

4.训练分类器：使用机器学习算法训练分类器，在GEE中可以使用Supervised Classification、Random Forest等算法。

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

6.分类结果输出：将分类结果输出为栅格图层或矢量要素图层，并进行后续的分析和应用。

以上是基本的土地分类步骤，具体的操作流程需要根据数据和研究目的进行适当的调整和修改。在GEE中进行土地分类需要一定的遥感和机器学习算法知识，建议在进行操作前进行相关学习和实践。

有关GEE JavaScript的分类请点击：
https://blog.csdn.net/qq_31988139/article/details/129870610
土地分类

土地分类案例

引用: Classification example for Landsat 8 imagery based on the scientfic work “MAD-MEX: Automatic Wall-to-Wall Land Cover Monitoring for the Mexican REDD-MRV Program Using All Landsat Data” by S.Gebhardt et. al 2014. Please find the link to the paper here: 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

Functions to derive vegetation indices and other raster operations
添加其它波段作为参与变脸给计算的变量

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')

Derive and visualize Tasseled Cap Transformation
推导并可视化缨帽变换

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