GDAL命令：一行代码转换坐标系

  本文介绍基于gdal模块，在命令行中通过GDAL命令的方式（不是Python或者C++代码，就是gdal模块自身提供的命令行工具），对栅格遥感影像数据加以投影，即将原本的地理坐标系转为投影坐标系的方法。

  首先明确一下本文的需求。我们现在有一个.tif格式的栅格遥感影像文件，其空间坐标系为GCS_WGS_1984，也就是WGS84，是一个地理坐标系；在ArcMap软件中将其打开，可以看到其空间坐标系及空间分辨率的单位（经纬度），如下图所示。

  我们现在希望，将这一景遥感影像加以投影，即将其坐标系由原本的地理坐标系转换为投影坐标系，目标投影坐标系为WGS_1984_UTM_Zone_48N，也就是一个UTM投影坐标系。在之前的文章中，我们也多次介绍过基于ArcGIS等软件，或者GEE等在线平台，直接或间接地实现矢量、栅格数据投影（或者重投影）的具体方法，大家可以参考文章ArcGIS矢量图层投影与地理坐标系转为投影坐标系——ArcMap，或者文章ArcMap通过模型构建器导出地理与投影坐标系转换的Python代码，再或者文章Google Earth Engine谷歌地球引擎地理坐标系、投影坐标系的变换与重投影加以查看。而本文，我们就介绍基于gdal模块（这个模块可以是大家单独配置的，也可以是在Python、C++等代码语言的环境下配置的），快速、方便地实现空间数据投影的方法。

  首先，我们需要配置好gdal模块。如果大家是用的Anaconda环境，那么就可以基于文章Anaconda环境配置GDAL的方法中介绍的方法，借着Python环境配置一下gdal模块；如果想通过其他方式配置gdal模块，那么参照gdal模块官网的介绍加以操作即可。

  配置gdal模块完毕后，我们打开电脑中的任意命令行工具。如果前期是在Python环境配置的gdal模块，那么就建议用Python环境下的命令行工具——否则，如果直接用操作系统自带的命令行工具，可能会出现由于环境变量配置不当导致的代码执行错误。例如，如果大家前期是在Anaconda环境的Python中配置的gdal模块，那么此时就打开Anaconda下属的Prompt工具即可；如下图所示，这两个Prompt工具选择任意一个均可。

  随后，在弹出的命令行中，我们首先cd进入存储有原文件（也就是待投影的栅格遥感影像文件）的路径下，然后输入如下的代码。

gdalwarp vegetation_type.tif result.tif -t_srs "EPSG:32648"

  其中，vegetation_type.tif就是原文件（待投影的文件）的名称，result.tif就是输出文件的名称；-t_srs表示目标坐标系（或者叫输出坐标系），其后面的参数就是我们期望的投影坐标系，随后的"EPSG:32648"就是WGS_1984_UTM_Zone_48N这个投影坐标系。大家可以在这个网站（https://epsg.io/）中，找到自己所需坐标系的EPSG编号。

  运行上述代码，如下图所示。

  其中，需要注意，我们也可以不cd进入存储有原文件（也就是待投影的栅格遥感影像文件）的路径，但那样就必须在上述代码的前2个参数中，将栅格遥感影像文件的名称用完整的绝对路径来表示；否则就会如上图紫色框上方的那个报错一样，找不到指定的文件。

  此外，需要注意的是，大家执行上述代码后，可能会出现ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db这个错误提示，如下图所示。

  遇到这种情况，我们就需要首先找到配置gdal模块时的路径，并在其中找到proj这个文件夹；因为我这里是在Anaconda环境的Python中配置的，所以就在Anaconda环境的Library文件夹找这个proj文件夹即可。随后，按照Windows环境变量的设置方法提到的方法，在系统变量中，新建一个名叫PROJ_LIB的变量，并将proj这个文件夹的路径作为其值。如下图所示。

  随后，再运行上述代码，即可不再报错。

  此时，如果我们用ArcGIS打开结果文件，可以看到其已经完成了投影，坐标系已经是WGS_1984_UTM_Zone_48N，且空间分辨率的单位为米；如下图所示。

  以上，我们利用了gdal模块提供的一个命令行工具——gdalwarp命令，实现了栅格图像投影的需求。gdal模块提供的这些命令行工具，可以在命令提示符或终端中执行，就不需要我们再写Python、C++等语言的代码了，所以比较方便。这些命令行工具通常作为gdal模块的一部分提供——在正确安装gdal模块后，其会自动添加到系统的环境变量中，以便在任何命令行工具里执行这些命令。

  除了上述命令行工具，按道理我们还可以用Python代码的方式，基于gdal模块提供的Python语言的API——gdal.Warp()函数，或者gdal.Translate()函数等，来实现栅格投影的需求；如以下代码所示。

# -*- coding: utf-8 -*-
"""
Created on Mon Feb 19 22:48:16 2024

@author: fkxxgis
"""
import os
from osgeo import gdal

os.environ["PROJ_LIB"] = r"C:\ProgramData\anaconda3\Library\share\proj";
os.environ["GDAL_DATA"] = r"C:\ProgramData\anaconda3\Library\share"

original_file_path = r"F:\Data_Reflectance_Rec\Type\vegetation_type.tif"
projected_file_path = r"F:\Data_Reflectance_Rec\Type\vegetation_type_pro1.tif"

target_projection = 'EPSG: 32648'

# gdal.Warp(projected_file_path, original_file_path, dstSRS = target_projection, targetAlignedPixels = True)
gdal.Translate(projected_file_path, original_file_path, projWin = None, outputSRS = target_projection)

  但是，经过尝试，这两个函数在我这里都行不通。其中，第一个gdal.Warp()函数在我这里会出现TypeError: in method 'wrapper_GDALWarpDestName', argument 4 of type 'GDALWarpAppOptions *'这样的报错，如下图所示。

  据说出现这个报错的原因是gdal模块自身版本的问题，所以可能还不太好解决。而对于第二个gdal.Translate()函数，其在我这里虽然可以不报错地执行代码，但是得到的栅格遥感影像结果文件还是地理坐标系，依然没有被投影。针对上述这些问题，也加以了多次尝试，但一直得不到正确的结果，只好作罢，最后发现还是用GDAL命令的方式，更加方便、快捷一些。

  至此，大功告成。