在地理数据处理中，矢量裁剪栅格是一个非常重要的操作，它可以帮助我们提取感兴趣的区域并获得更精确的分析结果。其重要性包括：
- 区域限定：地球科学研究通常需要关注特定的地理区域。通过矢量裁剪栅格，我们可以将栅格数据限定在特定的研究区域内，去除不相关的数据，从而减少计算量和提高分析效率。
- 空间一致性：栅格数据通常具有均匀的空间分辨率，但我们可能只对特定区域感兴趣。通过矢量裁剪栅格，我们可以将栅格数据限定在感兴趣的区域内，确保分析结果在这个区域内具有更高的精度和一致性。
- 数据匹配：在地球科学研究中，我们通常需要将不同来源的数据进行比较和分析。通过矢量裁剪栅格，我们可以将不同数据源的栅格数据裁剪为相同的区域，使它们具有相同的空间范围和分辨率，从而方便进行数据匹配和比较。
- 空间分析：地球科学研究中常常需要进行各种空间分析，如土地利用变化、水文模拟、地质建模等。通过矢量裁剪栅格，我们可以将栅格数据限定在特定区域内，以便更准确地进行空间分析，并获得与实际情况更符合的结果。
简单来说，制作数据的人要尽可能的考虑到多个用户的需求，所以往往数据的范围比较大，比如全国范围、全球范围，但是使用数据的人往往有自己聚焦的研究区，所以其 “研究区” 更明确。所以，在进行科研或者实际工作时，我们一般会裁剪出自己的目标区域数据，减少数据量，也方便分析和挖掘数据。
🌰 假设你获取到了全中国多年的气温栅格数据，而你的研究区是北京市，老板让你求北京市的平均温度序列，后续还要用北京市的数据进行其他分析。这时候我们第一步自然是先把我们需要的数据裁剪出来，既方便处理，也方便出图。
裁剪栅格数据一般有两种情况：
- 依据范围裁剪：有大概的范围，比如有研究区的经纬度范围，需要根据相应的经纬度范围来进行裁剪
- 基于矢量数据裁剪：有研究区的矢量数据，那么我们可以使用矢量数据来进行裁剪
可能用到的包
import rasterio as rio
from rasterio.plot import show
from rasterio.windows import Window
from rasterio.transform import from_origin
import geopandas as gpd
from rasterio.mask import mask读取栅格数据并简单可视化：
dataset = rio.open('/home/mw/input/GeoPy17233/GeoPy1/temp_grid_84/tem2015_84.tif')
show(dataset)我们读取的数据是整个中国范围的，可以通过下面的代码查看数据的范围：
dataset.boundsBoundingBox(left=65.45935893761649, bottom=13.120243022441713, right=137.45935893761649, top=56.01024302244171)
依据经纬度范围裁剪
假设我们的研究区是湖南省，借助搜索引擎，发现其经纬度范围为东经 108.786106 ~ 114.256514，北纬 24.643089 ~ 30.1287。
我们用下面的代码对上面的数据进行裁剪并且保存：
# 定义经纬度范围,湖南省范围
min_lon = 108
max_lon = 115
min_lat = 24
max_lat = 31# 打开栅格数据
with rio.open('/home/mw/input/GeoPy17233/GeoPy1/temp_grid_84/tem2015_84.tif') as src:# 获取栅格数据的元数据和变换信息profile = src.profiletransform = src.transform# 计算裁剪后的栅格数据的行列范围col_start, row_start = ~transform * \(min_lon, max_lat) # ~transform表示坐标转换的操作符号col_end, row_end = ~transform * (max_lon, min_lat)# 计算裁剪后的栅格数据的宽度和高度width = int(col_end - col_start)height = int(row_end - row_start)# 创建裁剪后的栅格数据数组cropped_data = src.read(window=Window(col_start, row_start, width, height))# 更新元数据profile.update(width=width, height=height, transform=from_origin(min_lon, max_lat, profile['transform'][0], -profile['transform'][4]))# 创建新的栅格数据文件并写入裁剪后的数据with rio.open('/home/mw/temp/tem2015_84_hunan.tif', 'w', **profile) as dst:dst.write(cropped_data)# 可视化裁剪后的数据
dataset = rio.open('/home/mw/temp/tem2015_84_hunan.tif')
show(dataset)
从上面结果可以看出我们已经成功裁剪出了我们需要范围的数据。
# 定义经纬度范围,湖南省范围
min_lon = 108
max_lon = 115
min_lat = 24
max_lat = 31# 打开栅格数据
with rio.open('/home/mw/input/GeoPy17233/GeoPy1/temp_grid_84/tem2015_84.tif') as src:# 获取栅格数据的元数据和变换信息profile = src.profiletransform = src.transform# 计算裁剪后的栅格数据的行列范围col_start, row_start = ~transform * \(min_lon, max_lat) # ~transform表示坐标转换的操作符号col_end, row_end = ~transform * (max_lon, min_lat)# 计算裁剪后的栅格数据的宽度和高度width = int(col_end - col_start)height = int(row_end - row_start)# 创建裁剪后的栅格数据数组cropped_data = src.read(window=Window(col_start, row_start, width, height))# 更新元数据profile.update(width=width,height=height,transform=from_origin(min_lon,max_lat,profile['transform'][0],profile['transform'][4])) # 提示：仔细观察这一行~# 创建新的栅格数据文件并写入裁剪后的数据with rio.open('/home/mw/temp/tem2015_84_hunan_wrong.tif', 'w', **profile) as dst:dst.write(cropped_data)# 可视化裁剪后的数据
dataset = rio.open('/home/mw/temp/tem2015_84_hunan_wrong.tif')
show(dataset)
用矢量数据来裁剪
下面我们再来学习如何用矢量文件来进行裁剪。
以北京市为例：
# 打开栅格数据
with rio.open('/home/mw/input/GeoPy17233/GeoPy1/temp_grid_84/tem2015_84.tif') as src:# 打开矢量文件beijing = gpd.read_file('/home/mw/input/GeoPy17233/GeoPy1/beijing.shp')# 对栅格数据进行矢量裁剪clipped_data, clipped_transform = mask(src, beijing.geometry, crop=True)# 获取裁剪后的元数据clipped_profile = src.profileclipped_profile.update(transform=clipped_transform,width=clipped_data.shape[2],height=clipped_data.shape[1])# 创建新的栅格数据文件并写入裁剪后的数据with rio.open('/home/mw/temp/tem2015_84_beijing.tif', 'w', **clipped_profile) as dst:dst.write(clipped_data)# 可视化结果
dataset = rio.open('/home/mw/temp/tem2015_84_beijing.tif')
show(dataset)
可以看出我们已经成功的用北京市的矢量裁剪出了我们想要的结果。
