【Python实例】获取指定经纬度范围内的栅格数据（tif）

在地理信息系统（GIS）中，经纬度是一种常用的坐标系统，用于在地球上定位位置。在许多应用中，我们可能需要根据给定的经纬度范围获取特定区域内的栅格位置。本文将介绍如何使用Python来实现这一目标（类似于GIS中的掩膜提取）。

数据准备

确定经纬度范围

首先，我们需要确定所需区域的经纬度范围。经度表示东西方向上的位置，介于-180到180之间，其中0代表本初子午线。纬度表示南北方向上的位置，介于-90到90之间，其中0代表赤道。
例如，我们希望获取北纬30度到40度，东经100度到110度范围内的栅格位置。

安装必要的库

在开始编写代码之前，我们需要安装一些Python库来处理地理空间数据。其中最常用的是GeoPandas和Shapely。

可以使用以下命令来安装它们：

pip install geopandas shapely

基于Python获取指定区域内栅格数据（Mask）

指定经纬度范围内

代码功能：根据全球不透水面积裁剪得到指定经纬度范围内的栅格数据

图形如下：

Python代码如下：

# 设置要裁剪的经纬度范围
min_lon = 102.0  # 最小经度
max_lon = 123.0  # 最大经度
min_lat = 14.0   # 最小纬度
max_lat = 33.0   # 最大纬度

# 导入重采样后的不透水面积数据
ISA_resampled_file = "D:/0 DataBase/5 Imperious surface area/ISA_resampled_500m.tif"
ISA_dataset = gdal.Open(ISA_resampled_file, gdal.GA_ReadOnly)

if ISA_dataset is None:
    raise Exception(f"Unable to open file+{ISA_resampled_file}")

# 获取数据的地理变换参数
geo_transform = ISA_dataset.GetGeoTransform()
top_left_x = geo_transform[0]
pixel_width = geo_transform[1]
top_left_y = geo_transform[3]
pixel_height = geo_transform[5]

# 计算裁剪区域的像素范围
# 计算左上角和右下角的像素位置
ul_x = int((min_lon - top_left_x) / pixel_width)
ul_y = int((top_left_y - max_lat) / abs(pixel_height))
lr_x = int((max_lon - top_left_x) / pixel_width)
lr_y = int((top_left_y - min_lat) / abs(pixel_height))

# 使用 gdal.Translate 裁剪数据
ISA_clipped_file = "D:/0 DataBase/5 Imperious surface area/ISA_clipped.tif"
gdal.Translate(
    ISA_clipped_file,
    ISA_dataset,
    srcWin=[ul_x, ul_y, lr_x - ul_x, lr_y - ul_y]  # 裁剪区域的像素范围
)

# 读取裁剪后的数据
ISA_clipped_dataset = gdal.Open(ISA_clipped_file, gdal.GA_ReadOnly)
band = ISA_clipped_dataset.GetRasterBand(1)
ISA_clipped_data = band.ReadAsArray()

# 打印裁剪后的数据大小
print("裁剪后的不透水面积数据大小:", ISA_clipped_data.shape)

指定研究区（其他tif文件）

Python代码如下：

# 导入重采样后的不透水面积数据
ISA_resampled_file = "D:/0 DataBase/5 Imperious surface area/ISA_resampled_500m.tif"
ISA_dataset = gdal.Open(ISA_resampled_file, gdal.GA_ReadOnly)

if ISA_dataset is None:
    raise Exception(f"Unable to open file+{ISA_resampled_file}")

# 计算不透水面积数据的地理变换参数
ISA_geo_transform = ISA_dataset.GetGeoTransform()
ISA_top_left_x = ISA_geo_transform[0]
ISA_pixel_width = ISA_geo_transform[1]
ISA_top_left_y = ISA_geo_transform[3]
ISA_pixel_height = ISA_geo_transform[5]

# 计算裁剪区域的左上角和右下角经纬度
min_lon = top_left_x
max_lon = top_left_x + (LUCC_width * pixel_width)
max_lat = top_left_y
min_lat = top_left_y + (LUCC_height * pixel_height)

# 计算裁剪区域的像素范围
ul_x = int((min_lon - ISA_top_left_x) / ISA_pixel_width)
ul_y = int((ISA_top_left_y - max_lat) / abs(ISA_pixel_height))
lr_x = int((max_lon - ISA_top_left_x) / ISA_pixel_width)
lr_y = int((ISA_top_left_y - min_lat) / abs(ISA_pixel_height))

# 计算对齐的左上角坐标
aligned_top_left_x = top_left_x
aligned_top_left_y = top_left_y

# 进行裁剪和对齐
ISA_clipped_file = "D:/0 DataBase/5 Imperious surface area/ISA_clipped_aligned.tif"

gdal.Warp(
    ISA_clipped_file,
    ISA_dataset,
    width=LUCC_width,  # 输出与 LUCC 数据相同的宽度
    height=LUCC_height,  # 输出与 LUCC 数据相同的高度
    outputType=gdal.GDT_Float32,  # 输出数据类型
    xRes=pixel_width,  # 与 LUCC 相同的分辨率
    yRes=abs(pixel_height),  # 与 LUCC 相同的分辨率
    dstSRS=LUCC_dataset.GetProjection(),  # 使用 LUCC 数据的投影
    outputBounds=[aligned_top_left_x, aligned_top_left_y + LUCC_height * pixel_height,
                  aligned_top_left_x + LUCC_width * pixel_width, aligned_top_left_y],  # 输出边界
    resampleAlg=gdal.GRA_Bilinear  # 使用双线性插值
)

# 读取裁剪和对齐后的数据
ISA_clipped_dataset = gdal.Open(ISA_clipped_file, gdal.GA_ReadOnly)
band = ISA_clipped_dataset.GetRasterBand(1)
ISA_clipped_data = band.ReadAsArray()

# 打印裁剪和对齐后的数据大小
print("裁剪和对齐后的不透水面积数据大小:", ISA_clipped_data.shape)

参考