【Python实例】Python读取并绘制tif数据

【Python实例】Python读取并绘制tiff数据

Python实例-以全球不透水面积数据为例
- 数据准备：全球不透水面积数据
- 基于gdal库绘制tif图
- 基于Rasterio库绘制tif图
参考

GeoTIff 是一个标准的.tif 文件或是一个图像文件格式，它包含了一些额外的空间信息，这些信息被当成附属信息(tag），集成在.tif文件内。这些附属信息包含了空间范围、地理参考系统（CRS）、分辨率，以及每一个像素的值。基于此，GeoTiff是一种非常理想的遥感影像和航空相片分发格式。

Python实例-以全球不透水面积数据为例

数据准备：全球不透水面积数据

以全球不透水面积数据为例，数据下载-全球1985-2022逐年不透水层数据（Version 2024）
在这里插入图片描述

资源详情
GAIA (1985-2024)为GAIA(1985-2018)的更新版本。代表性年份(1985年、1990年、1995年、2000年、2005年、2010年和2015年)的GAIA数据总体精度超过90%。GAIA的时间趋势与其他局域、地区和全球尺度的数据集一致。更多细节可以在相关论文中找到(Gong et al., 2020)。
元数据信息
参考坐标系: WGS84 - EPSG:4326
空间分辨率: 30米
空间覆盖范围: (xmin, xmax, ymin, ymax) - (-180, 180, -60, 80)
时间分辨率: 1985-2024 (逐年)
数据格式: GeoTiff
像元值：像元值表示城市出现的频率（0-40）
0: 非城市区域
1: 2024年新增城市区域
2: 2023年新增城市区域
…
40: 1985年及之前存在的城市区域
数据格网
GAIA数据基于渔网格网以5°×5°的数据图幅存储。格网文件以“fishnet_shp”目录提供，可供查找感兴趣区域。
命名规则
每个数据图幅以“GAIA_1985_2024_{longitude}_{latitude}.tif”格式命名，其中，longitude-latitude 指格网左上角经纬度值。另外，每个图幅包含了30米缓冲区范围以便图幅间可以无缝镶嵌。

在GIS中打开界面如下：
在这里插入图片描述

基于gdal库绘制tif图

gdal是最流行的GeoTiff处理库，拥有由C++编写的方法和类。绝大部分的库，诸如georaster等，也是在运用GDAL的基础上，开发出符合Python风格的接口。

from osgeo import gdal
import matplotlib.pyplot as pltdataset = gdal.Open("D:/6 Python Codes/GBA_Data_Process/GAIA/GAIA_1985_2024_110_25.tif", gdal.GA_ReadOnly)
band = dataset.GetRasterBand(1) # 波段序号从1开始，而不是0
plt.figure(figsize=(10, 10))
plt.imshow(band.ReadAsArray())
plt.show()

绘制图形如下：
在这里插入图片描述

设置经纬度等细节，代码如下：

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
from matplotlib.colors import LinearSegmentedColormap# 设置字体为 Times New Roman
plt.rcParams['font.family'] = 'Times New Roman'# TIFF文件路径
tiff_file = "D:/6 Python Codes/GBA_Data_Process/GAIA/GAIA_1985_2024_110_25.tif"
dataset = gdal.Open(tiff_file, gdal.GA_ReadOnly)# 检查数据集是否成功打开
if dataset is None:raise Exception("Unable to open file")# 读取图像数据
band = dataset.GetRasterBand(1)  # 获取第一波段(波段序号从1开始，而不是0)
image_data = band.ReadAsArray()# 获取仿射变换参数
geo_transform = dataset.GetGeoTransform()# 提取仿射变换参数
top_left_x = geo_transform[0]
pixel_width = geo_transform[1]
top_left_y = geo_transform[3]
pixel_height = geo_transform[5]# 计算经纬度范围
height, width = image_data.shape
lon = np.linspace(top_left_x, top_left_x + pixel_width * width, width)
lat = np.linspace(top_left_y, top_left_y + pixel_height * height, height)# 创建橙色渐变色
cmapOrange = LinearSegmentedColormap.from_list("orange_gradient", ["white", "orange"])# 使用Matplotlib绘制地图
plt.figure(figsize=(10, 8))
img = plt.imshow(image_data, extent=(lon.min(), lon.max(), lat.min(), lat.max()), cmap= cmapOrange)
cbar = plt.colorbar(img, label='Pixel values')
cbar.set_label('Pixel values', fontsize=14)plt.xlabel('Longitude (°)', fontsize=16)
plt.ylabel('Latitude (°)', fontsize=16)
plt.title('GAIA 1985 Year', fontsize=16)# 设置坐标轴刻度值字体大小
plt.tick_params(axis='both', labelsize=14)# 自定义坐标轴刻度值，后面加上°
def format_ticks(x, pos):return f'{x:.2f}°'# 应用自定义格式化函数
plt.gca().xaxis.set_major_formatter(FuncFormatter(format_ticks))
plt.gca().yaxis.set_major_formatter(FuncFormatter(format_ticks))plt.show()

绘制图形如下：
在这里插入图片描述

基于Rasterio库绘制tif图

Rasterio由mapbox团队开发，它提供了一系列用于读取地理空间数据的Python接口，可配合Matplotlib库使用。

代码如下：

import rasterio
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.ticker import FuncFormatter# 设置字体为 Times New Roman
plt.rcParams['font.family'] = 'Times New Roman'# 打开TIFF文件
tiff_file = "D:/6 Python Codes/GBA_Data_Process/GAIA/GAIA_1985_2024_110_25.tif"
with rasterio.open(tiff_file) as dataset:# 读取图像数据image_data = dataset.read(1)# 获取仿射变换参数transform = dataset.transform# 计算经纬度范围
height, width = image_data.shape
lon = np.linspace(transform.c, transform.c + transform.a * width, width)
lat = np.linspace(transform.f, transform.f + transform.e * height, height)# 创建橙色渐变色
cmap = LinearSegmentedColormap.from_list("orange_gradient", ["white", "orange"])# 绘制地图
plt.figure(figsize=(10, 8))
img = plt.imshow(image_data, extent=(lon.min(), lon.max(), lat.min(), lat.max()), cmap=cmap)
cbar = plt.colorbar(img, label='Pixel values')
cbar.set_label('Pixel values', fontsize=12)# 设置坐标轴标签
plt.xlabel('Longitude (°)', fontsize=12)
plt.ylabel('Latitude (°)', fontsize=12)
plt.title('TIFF Map', fontsize=14)# 设置坐标轴刻度值字体大小
plt.tick_params(axis='both', labelsize=10)# 自定义坐标轴刻度值，后面加上°
def format_ticks(x, pos):return f'{x:.2f}°'# 应用自定义格式化函数
plt.gca().xaxis.set_major_formatter(FuncFormatter(format_ticks))
plt.gca().yaxis.set_major_formatter(FuncFormatter(format_ticks))plt.show()