如何使用Python读取CropScape的DBF文件并绘制作物分布可视化图

解决方案 #

你当前读取的CDL_2020_clip_xxx.tif.vat.dbf是栅格数据的属性对照表，仅存储栅格数值对应的作物类型、配色规则，没有实际的空间坐标信息，因此无法单独用这个DBF文件绘图，必须搭配同目录下前缀完全一致的TIF栅格文件完成可视化，具体实现代码如下：

依赖导入

import geopandas as gpd
import rasterio
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch
import numpy as np

步骤1：基于DBF构建自定义颜色映射表

# 读取DBF属性表
filename = 'CDL_2020_clip_20210826173123_1609691609.tif.vat.dbf'
dbf = gpd.read_file(filename)

# 生成颜色映射表，VALUE对应栅格值，RGBA对应显示颜色透明度
max_val = dbf['VALUE'].max()
color_list = np.zeros((max_val + 1, 4))
for _, row in dbf.iterrows():
    color_list[int(row['VALUE'])] = [row['RED'], row['GREEN'], row['BLUE'], row['OPACITY']]
cdl_cmap = ListedColormap(color_list)

# 提前获取玉米、大豆对应的分类值
corn_val = dbf[dbf['CLASS_NAME']=='Corn']['VALUE'].iloc[0]
soybean_val = dbf[dbf['CLASS_NAME']=='Soybean']['VALUE'].iloc[0]

步骤2：读取配套的TIF栅格数据

# 读取同前缀的TIF栅格文件
tif_path = filename.replace('.vat.dbf', '')
with rasterio.open(tif_path) as src:
    img = src.read(1) # CDL为单波段分类栅格，读取第一波段即可
    extent = [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top]
    crs = src.crs

步骤3：绘制可视化图件

fig, ax = plt.subplots(figsize=(10, 10), dpi=100)

# 可选：仅保留玉米、大豆两类显示，其余像素设为背景
img[~np.isin(img, [corn_val, soybean_val])] = 0

# 绘制栅格，使用官方配色
im = ax.imshow(img, extent=extent, cmap=cdl_cmap, vmin=0, vmax=max_val)

# 可选：叠加飞行路径等矢量数据，注意要转成和栅格一致的坐标系
# flight_path = gpd.read_file('飞行路径文件路径.shp').to_crs(crs)
# flight_path.plot(ax=ax, color='red', linewidth=1.5, label='飞行路径')

# 添加自定义图例
legend_elements = [
    Patch(facecolor=color_list[corn_val][:3], label='玉米'),
    Patch(facecolor=color_list[soybean_val][:3], label='大豆')
]
ax.legend(handles=legend_elements, loc='upper right', fontsize=12)

ax.set_xlabel('经度', fontsize=12)
ax.set_ylabel('纬度', fontsize=12)
ax.set_title('2020年艾奥瓦州玉米大豆分布', fontsize=14)
plt.show()

内容的提问来源于stack exchange，提问作者Natasha