Python把遥感影像的投影坐标系转为地理坐标系

  本文介绍基于Python语言中的ArcPy模块，批量将多个遥感影像由投影坐标系转为地理坐标系的方法。

  在之前的文章中，我们介绍过将单独1景遥感影像的投影坐标系转为地理坐标系的方法，大家可以参考文章GDAL一行代码将投影坐标系变为地理坐标系；但是，这个方法对于少量遥感影像的场景比较有优势，而对于具有大量待处理遥感影像的批处理场景而言并不是很友好（主要是我不知道怎么在命令行写批量操作的脚本）。而在本文中，我们就介绍一下基于Python的ArcPy模块，实现上述需求的批量操作的方法；此外，如果大家需要反过来，将地理坐标系转为投影坐标系，也可以参考本文的方法，亦可以参考文章地理坐标系变换至投影坐标系：GDAL命令行，或者是文章ArcGIS图层投影与地理坐标系转为投影坐标系的方法。

  首先，明确一下本文的需求。现有一个文件夹，其中含有大量不同格式的文件，如下图所示。其中，我们希望将所有.tif格式的遥感影像文件，由原本的投影坐标系转为地理坐标系。

  如果我们用ArcMap等软件将上述遥感影像文件打开，可以看到其坐标系为投影坐标系WGS_1984_UTM_Zone_48N，且单位为米；如下图所示。我们希望将这个投影坐标系转为地理坐标系，本文就以转为WGS_1984这个地理坐标系为例来介绍。

  本文所需要的代码如下所示。

# -*- coding: utf-8 -*-
"""
Created on Thu May  2 13:05:43 2024

@author: fkxxgis
"""

import os
import arcpy

input_folder = r"E:\04_Reconstruction\99_MODIS\new_data\GF"
output_folder = r"E:\04_Reconstruction\99_MODIS\new_data\GF_WGS84"
arcpy.env.workspace = input_folder
target_spatial_reference = arcpy.SpatialReference(4326)

input_files = arcpy.ListRasters("*.tif")

for input_file in input_files:
    output_file = os.path.join(output_folder, os.path.splitext(input_file)[0] + "_WGS84.tif")

    arcpy.ProjectRaster_management(input_file, output_file, target_spatial_reference)
    print input_file + " finished!"

  首先，我们导入所需使用的库。其中，os库用于文件和路径操作，arcpy则用于处理GIS数据。

  随后，我们定义输入与输出文件的路径。其中，使用r前缀是为了将字符串视为原始字符串，从而避免转义字符的处理。同时，使用arcpy.env.workspace将工作空间设置为输入文件夹，这样在后续处理数据时可直接使用文件名而不需使用完整的路径。随后，通过arcpy.SpatialReference()创建一个代表WGS84地理坐标系的空间参考对象。大家在使用代码时，将这里的4326坐标系编号修改为自己所需目标坐标系的EPSG编号即可；其中，可以在这个网站中，找到自己所需坐标系的EPSG编号。

  接下来，首先使用arcpy.ListRasters()获取输入文件夹中所有的.tif格式文件，随后使用for循环遍历每一个输入文件。在循环中，通过os.path.join()和os.path.splitext()构建输出文件的路径，将输入文件的文件名作为基础，并在结尾加上_WGS84作为标识。然后，使用arcpy.ProjectRaster_management()函数将输入文件投影到WGS84地理坐标系，并将结果保存到输出文件中。最后，通过print语句打印处理完成的文件名。

  执行上述代码，即可得到如下图所示的界面。我这里在print输出内容时，文件名称和finished!之间忘记加空格了，大家理解即可。

  此外，可以在结果文件夹中看到如下图所示的结果文件。

  用ArcMap等软件将上述遥感影像文件打开，可以看到其坐标系为地理坐标系WGS_1984，且单位为度；如下图所示。

  至此，大功告成。