遥感云平台-GEE获取CLCD土地利用数据（python+土地利用转移矩阵+桑基图）

内容介绍

本期内容包括调用CLCD数据集并显示，对土地覆盖类型进行重分类，统计每一类型的面积，求解较小区域范围的土地利用转移矩阵并绘制相应的桑基图。
环境配置：Vscode+Jupyter notebook+gee+geemap+python3.10

导入所需要的包并进行初始化

import ee
import os
import geemap
import pandas as pd
import plotly.graph_objects as go
import numpy as np
from matplotlib import rc
from geemap import cartoee
import matplotlib.pyplot as plt
from proplot import rc
import matplotlib.patches as mpatches
from sklearn.metrics import confusion_matrix
ee.Initialize()

Map = geemap.Map()
Map.add_tile_layer(
    url="https://mt1.google.com/vt/lyrs=m&x={x}&y={y}&z={z}&scale=2",
    name="HD Google Basemap",
)
roi = ee.FeatureCollection("")# 自行更改研究区域内的矢量文件
Map.addLayer(roi,{'color':'grey'},'Beijing_GB');
Map.centerObject(roi,10);

vis_CLCD = {
    'min': 1,
    'max': 9,
    'palette': ['FAE39C', '446F33', '33A02C',
        'ABD37B', '1E69B4', 'A6CEE3',
        'CFBDA3', 'E24290', '289BE8',]
}
vis_reCLCD = {
    'min': 1,
    'max': 6,
    'palette': ['FAE39C', '446F33',
        'ABD37B', '1E69B4','CFBDA3', 'E24290']
}
style = {"color": "dff2ff", "width": 1.5, "lineType": "solid", "fillColor": "00000000"}

调用CLCD数据集并重分类

此处参考GEE之分区统计与重分类-以CLCD数据集以及武汉市各区为例

srcFolder = 'projects/lulc-datase/assets/LULC_HuangXin/'

imgList = ee.List([])
for year in range(2000, 2020):  
    tmpImg = ee.Image(srcFolder + 'CLCD_v01_' + str(year)).clip(roi).set('year', year);
    imgList = imgList.add(tmpImg)

imgCollection = ee.ImageCollection.fromImages(imgList)
Map.addLayer(imgCollection.first(),vis_CLCD, "test")

reclassList = ee.List([])
for year in range(2000, 2020):  
    tmpImg = ee.Image(srcFolder + 'CLCD_v01_' + str(year)).clip(roi)\
        .remap(
      [1, 2, 3, 4, 5, 6, 7, 8, 9], # 原始类别
      [1, 2, 2, 3, 4, 4, 5, 6, 4]  # 重分类后的类别(1 耕地 2 林地 3 草地 4 水域 5 未利用地 6 建设用地) 可自行重分类
    ).set('year', year);
    reclassList = reclassList.add(tmpImg)
reclassimgCollection = ee.ImageCollection.fromImages(reclassList)
Map.addLayer(reclassimgCollection .first(),vis_reCLCD, "test")

制图与统计

def get_landcover_by_year(year, image_collection):
    yearstr = ee.String(str(year)).getInfo()
    landcover = image_collection.filter(ee.Filter.eq('year', year)).first()
    return yearstr, landcover
def plot_landcover(bbox, landcover, yearstr, roi, vis_reCLCD, style):
    plt.rcParams['font.family'] = 'Times New Roman'
    rc['tick.labelsize'] = 20
    rc["axes.labelsize"] = 20
    rc["axes.labelweight"] = "bold"
    rc["tick.labelweight"] = "bold"

    landcovervis = landcover.visualize(**vis_reCLCD)
    imgBlend = landcovervis.blend(roi.style(**style))

    fig = plt.figure(figsize=(20, 22),facecolor='white')
    ax = cartoee.get_map(imgBlend, region=bbox)

    gl = cartoee.add_gridlines(ax, interval=[1, 1], linestyle="--", color='lightgray', zorder=0)

    ax.set_title(label='Land Cover '+ yearstr, fontsize=20, pad=5)

    labels = ['Cropland', 'Forest', 'Grassland', 'Water', 'Unused land', 'Construction land']
    colors = ['#FAE39C', '#446F33', '#ABD37B', '#1E69B4', '#CFBDA3', '#E24290']
    patches = [mpatches.Patch(color=color, label=label) for color, label in zip(colors, labels)]
    plt.legend(handles=patches, loc='lower right', fontsize=15)
    legend = plt.legend(handles=patches, loc='lower right', fontsize=15, handleheight=2.5, handlelength=2.5)
    
    plt.show()
    output_path = "C:/Users/UWPK/Desktop/CLCD"+yearstr+".png"
    fig.savefig(output_path, dpi=600, bbox_inches='tight', pad_inches=0.1)
def download_landcover(yearstr,landcover,roi):
    out_dir = 'E:\\Download'
    filename = os.path.join(out_dir, "landcover"+yearstr+".tif")
    geemap.download_ee_image(
        landcover, filename=filename, scale=30, region=roi.geometry(), crs='EPSG:4326'
    )
def statistical(landcover,roi):
    areaImage = ee.Image.pixelArea().addBands(landcover)
    clcdArea = areaImage.reduceRegion(**{
        'reducer': ee.Reducer.sum().group(**{
            'groupField': 1,
            'groupName': 'remapped'
        }),
        'geometry': roi.geometry(),
        'scale': 30,
        'maxPixels': 1e13 
    })
    land_use_mapping = {
        1: 'Cropland',
        2: 'Forest',
        3: 'Grassland',
        4: 'Water',
        5: 'Unused land',
        6: 'Construction land'
    }
    clcdarea = clcdArea.getInfo()
    clcddata = clcdarea['groups']
    df = pd.DataFrame(clcddata)
    df.columns = ['Category', 'Total(m²)']
    df['Category'] = df['Category'].map(land_use_mapping)
    pd.set_option('display.float_format', '{:.4f}'.format)
    return df

# 从上面的reclassimgCollection中获取数据，yearstr是为了方便命名
yearstr2000, landcover2000 = get_landcover_by_year(2000, reclassimgCollection)
yearstr2019, landcover2019 = get_landcover_by_year(2019, reclassimgCollection)

# 制图
bbox = [115.203,39.389,117.532,41.09] # 改为研究区域的经纬度范围
plot_landcover(bbox, landcover2019, yearstr2019, roi, vis_reCLCD, style)

# 下载
# download_landcover(yearstr2000, landcover2000)

# 统计各类型的面积
df2000 = statistical(landcover2000,roi)
print('landcover2000',df2000)
df2019 = statistical(landcover2019,roi)
print('landcover2019',df2019)

制图结果：

统计结果：

小范围的土地转移矩阵和桑基图

由于geemap中的ee_to_numpy函数有数据大小限制，需要注意。

def numpydata(landcover,roi):
    lcdata = geemap.ee_to_numpy(landcover, region=roi,scale=30)
    lcdata_array = np.asarray(lcdata).flatten()
    return lcdata_array
roi2 = ''# 定义研究区域
lcdata2000_array = numpydata(landcover2000,roi2)
lcdata2019_array = numpydata(landcover2019,roi2)
transition_area = confusion_matrix(lcdata2000_array, lcdata2019_array)
print(transition_area)

# 绘制桑基图
node_labels = ['Cropland', 'Forest', 'Grassland', 'Water', 'Unused land', 'Construction land',
               'Cropland', 'Forest', 'Grassland', 'Water', 'Unused land', 'Construction land']
df = pd.DataFrame(transition_area)
sources = []
targets = []
values = []
for i in range(len(df)):
    for j in range(len(df.columns)):
        sources.append(i)
        targets.append(j + len(node_labels) // 2)
        values.append(df.iloc[i, j])
fig = go.Figure(data=[go.Sankey(
    node = dict(
      pad = 15,
      thickness = 20,
      line = dict(color = "black", width = 0.5),
      label = node_labels
    ),
    link = dict(
      source = sources,
      target = targets,
      value = values
    )
)])
fig.show()

土地转移矩阵结果：

桑基图：