欢迎访问宙启技术站
智能推送

使用Python和GDAL进行地理数据的重投影

发布时间:2023-12-27 16:04:10

以下是使用Python和GDAL库进行地理数据的重投影的示例代码:

首先,确保已经安装了GDAL库。可以使用pip安装,命令为:pip install GDAL

接下来,导入所需的库和模块:

from osgeo import gdal, osr

# 输入文件路径
input_file = "input.tif"

# 输出文件路径
output_file = "output.tif"

# 输入投影
input_projection = "+proj=latlong +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

# 输出投影
output_projection = "+proj=utm +zone=51 +datum=WGS84 +units=m +no_defs"

接下来,读取输入文件的投影和地理转换信息:

# 打开输入文件
input_dataset = gdal.Open(input_file)

# 获取输入文件的投影信息
input_projection_ref = input_dataset.GetProjection()

# 获取输入文件的地理转换信息
input_geotransform = input_dataset.GetGeoTransform()

然后,创建一个新的输出数据集,并设置输出投影和地理转换信息:

# 创建输出数据集
output_dataset = gdal.GetDriverByName("GTiff").Create(output_file, input_dataset.RasterXSize, input_dataset.RasterYSize, 1, gdal.GDT_Float32)

# 设置输出数据集的投影信息
output_dataset.SetProjection(output_projection)

# 设置输出数据集的地理转换信息
output_dataset.SetGeoTransform(input_geotransform)

接下来,进行地理数据的重投影:

# 创建输入和输出的空间参考对象
input_spatial_ref = osr.SpatialReference()
input_spatial_ref.ImportFromProj4(input_projection)

output_spatial_ref = osr.SpatialReference()
output_spatial_ref.ImportFromProj4(output_projection)

# 创建一个坐标转换对象,用于进行地理数据的重投影
coord_transform = osr.CoordinateTransformation(input_spatial_ref, output_spatial_ref)

# 读取输入数据的地理坐标,并进行重投影
for y in range(input_dataset.RasterYSize):
    for x in range(input_dataset.RasterXSize):
        # 读取输入数据的地理坐标
        lon, lat, _ = gdal.ApplyGeoTransform(input_geotransform, x, y)

        # 进行地理数据的重投影
        new_lon, new_lat, _ = coord_transform.TransformPoint(lon, lat)

        # 将地理坐标转换为像素坐标
        new_x, new_y = gdal.ApplyGeoTransform(output_dataset.GetGeoTransform(), new_lon, new_lat)

        # 读取输入数据的像素值
        pixel_value = input_dataset.GetRasterBand(1).ReadAsArray(x, y, 1, 1)[0][0]

        # 将像素值写入输出数据集
        output_dataset.GetRasterBand(1).WriteArray(pixel_value, int(new_x), int(new_y))

最后,保存输出数据集,并关闭输入和输出文件:

# 保存输出数据集
output_dataset.FlushCache()

# 关闭输入和输出文件
input_dataset = None
output_dataset = None

这是一个简单的示例,演示了如何在Python中使用GDAL库进行地理数据的重投影。根据具体的需求,你可能需要调整代码以适应不同的数据格式和投影信息。希望这个例子可以帮助到你!