编写一个批量改投影并按照矢量范围裁剪栅格数据的代码

在Python中，你可以使用GDAL库来进行栅格数据的批量投影转换和矢量范围裁剪。GDAL是一个强大的开源库，用于地理空间处理。下面是一个基本的示例，假设你已经有了一个包含多个栅格文件的目录路径，并且有一个预设的新投影文件：

from osgeo import gdal, ogr

def batch_project_and_crop(input_dir, output_dir, new_projection_file):
    # 创建一个新的临时工作空间
    temp_dir = 'temp'
    if not os.path.exists(temp_dir):
        os.makedirs(temp_dir)

    # 获取新投影的信息
    driver = ogr.GetDriverByName('ESRI Shapefile')
    with driver.Open(new_projection_file, 0) as shp:
        src_proj = shp.GetSpatialRef()

    for file in os.listdir(input_dir):
        if file.endswith('.tif'):  # 检查文件是否是栅格文件
            input_path = os.path.join(input_dir, file)
            output_path = os.path.join(output_dir, file[:-4] + '_cropped.tif')

            # 打开原始栅格文件
            src_ds = gdal.Open(input_path)
            src_band = src_ds.GetRasterBand(1)

            # 创建目标文件，指定新的投影
            dst_driver = gdal.GetDriverByName('GTiff')
            dst_ds = dst_driver.CreateCopy(os.path.join(temp_dir, file), src_ds, options=['COMPRESS=LZW', 'TILED=YES'])
            dst_ds.SetProjection(src_proj.ExportToWkt())

            # 创建几何范围裁剪操作
            transform = dst_ds.GetGeoTransform()
            with ogr.Open(new_projection_file, update=True) as shp:
                geom = shp.GetLayer().GetSpatialRef().ExportToWkt()  # 使用shapefile的坐标系统
                clip_window = ogr.Geometry(ogr.wkbBox2D)
                clip_window.AddPoint(transform[0], transform[3])
                clip_window.AddPoint(transform[0] + transform[1]*src_ds.RasterXSize, transform[3] + transform[5]*src_ds.RasterYSize)
                gdal.Warp(output_path, dst_ds, cutlineDSName=new_projection_file, cropToCutline=True, warpMemoryLimit=2**30)

            # 删除临时文件
            src_ds = None
            dst_ds = None

    # 清理临时目录
    shutil.rmtree(temp_dir)

# 调用函数
batch_project_and_crop('/path/to/input/directory', '/path/to/output/directory', '/path/to/new_projection.shp')

这个代码会将输入目录下的每个栅格文件转换到新投影，并按照提供的矢量范围裁剪。