介绍两种按照shp行政边界切割tif文件的方式
1. gdal.Warp
from osgeo import gdal
shp_path="shp文件路径"
path="tif数据路径"
#读取shp
shp_dataset = gdal.OpenEx(shp_path, gdal.OF_VECTOR)
shp_layer = shp_dataset.GetLayer()
shp_srs = shp_layer.GetSpatialRef().ExportToWkt()
# 读取tif
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(path, len(new_im_data[0]), len(new_im_data), 1, gdal.GDT_Float32)
# 我的shp文件中有多个城市 所以我这里是循环切割每个城市
for city_layer in shp_layer:
print(city_layer.name)
# 执行裁剪 output_dataset裁剪后的数据 获取数据第一个参数写'',需要保存到文件直接写保存路径format可不写
output_dataset = gdal.Warp('',
# 原始数据
dataset,
# 保存到内存
format='MEM',
cropToCutline=True,
# srs可以自定义或者使用shp文件中srs
# dstSRS='EPSG:4326',
dstSRS=shp_srs,
#shp文件地址
cutlineDSName=shp_path,
# 查询条件 Name是shp文件中的字段 查询Name=XXX shp中只有一个city可以不使用这个参数
cutlineWhere=f'"Name"=\'{city_layer.name}\'',
#指定重采样的方法
resampleAlg='bilinear',
# 原始数据无效值
srcNodata=-99,
# 新数据设置无效值
dstNodata=np.nan
)
output_data = output_dataset.ReadAsArray()
2. rasterio.mask
import rasterio
from rasterio.mask import mask
import geopandas as gpd
# 打开 shapefile 文件
shapefile_path = 'path/to/shapefile.shp'
shapefile = gpd.read_file(shapefile_path)
# 打开要裁剪的 GeoTIFF 文件
tiff_path = 'path/to/input.tif'
with rasterio.open(tiff_path) as src:
# 使用 mask 函数进行裁剪
out_image, out_transform = mask(src, shapefile.geometry, crop=True)
out_meta = src.meta.copy()
# 更新裁剪后图像的元数据
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
# 将裁剪后的图像保存到新文件中
output_path = 'path/to/output.tif'
with rasterio.open(output_path, "w", **out_meta) as dest:
dest.write(out_image)