使用Python和GDAL给图片加坐标系的实现思路(坐标投影转换)

使用Python和GDAL给图片加坐标系可以通过以下步骤完成：

1. 安装GDAL：可以通过pip安装，命令为：pip install gdal。安装完毕后，在Python代码中用import gdal语句引入模块。

2. 读取图片：使用gdal.Open()函数打开需要添加坐标系的图片。如下所示：

```
from osgeo import gdal

filename = 'test.tif'
ds = gdal.Open(filename, gdal.GA_Update) # 打开文件，可读可写模式
```

此处使用了GDAL提供的Open()函数打开了一个名为test.tif的图片，并使用了GA_Update参数，表示可以进行读取和写入操作。

1. 添加坐标系：使用ds.SetGeoTransform()函数设置图片的投影信息。投影信息包括左上角的x坐标、每个像素的宽度（x分辨率）、旋转角度、左上角的y坐标、旋转角度、每个像素的高度（y分辨率）。示例如下：

geo_transform = (xmin, xres, x_rot, ymax, y_rot, yres)
ds.SetGeoTransform(geo_transform)

其中，xmin为图片左上角的x坐标，xres表示每个像素的宽度，ymax为图片左上角的y坐标，yres表示每个像素的高度，x_rot和y_rot为旋转角度。

1. 设置投影信息：通过ds.SetProjection()函数设置图片的投影信息。投影信息使用字符串格式表示，可以通过GDAL提供的EPSG库调用标准的投影信息。例如，设置图片的投影为WGS 84 / UTM zone 48N，则可以使用以下代码：

epsg_code = 'EPSG:32648'
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(int(epsg_code.split(':')[1]))
ds.SetProjection(target_srs.ExportToWkt())

上述代码使用epsg_code定义了投影信息代码，然后使用osr.SpatialReference()函数创建了空间参考。ImportFromEPSG()函数解析了epsg_code中的代码为整数型，然后将整数型的代码传递给空间参考中的函数，创建了“WGS 84 / UTM zone 48N”投影系的参考。

接下来，可以使用target_srs.ExportToWkt()函数将参考转换为标准的WKT（Well-Known Text）格式，并使用ds.SetProjection()函数设置图片的投影信息。

1. 关闭文件：完成以上设置后，需要使用ds.FlushCache()函数刷新数据，以确保修改生效。并且需要使用ds = None语句或del ds语句关闭文件。完整代码如下：

```
from osgeo import gdal
from osgeo import osr

filename = 'test.tif'
ds = gdal.Open(filename, gdal.GA_Update)

xmin, xres, x_rot, ymax, y_rot, yres = 100, 0.01, 0, 50, 0, -0.01
geo_transform = (xmin, xres, x_rot, ymax, y_rot, yres)

epsg_code = 'EPSG:32648'
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(int(epsg_code.split(':')[1]))
ds.SetProjection(target_srs.ExportToWkt())

ds.SetGeoTransform(geo_transform)
ds.FlushCache()
ds = None
```

这个代码将图片的左上角坐标设为(100,50)，每个像素的宽度和高度分别设为0.01和-0.01，即每个像素为正方形。这个代码还将图片的投影指定为“WGS 84 / UTM zone 48N”。

示例1：将一个无投影的tif文件添加坐标投影

```
from osgeo import gdal
from osgeo import osr

filename = 'test_no_proj.tif'
ds = gdal.Open(filename, gdal.GA_Update)

xmin, xres, x_rot, ymax, y_rot, yres = 100, 0.01, 0, 50, 0, -0.01
geo_transform = (xmin, xres, x_rot, ymax, y_rot, yres)

epsg_code = 'EPSG:32648'
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(int(epsg_code.split(':')[1]))
ds.SetProjection(target_srs.ExportToWkt())

ds.SetGeoTransform(geo_transform)
ds.FlushCache()
ds = None
```

将无痕渲染的NEM视图图像（test_no_proj.tif）从本地打开，添加了坐标投影并输出为新的文件（加坐标后得到的输出文件test.tif）。

示例2：修改一个带有投影信息的tif文件的投影

```
from osgeo import gdal
from osgeo import osr

filename = 'test.tif'
ds = gdal.Open(filename, gdal.GA_Update)

xmin, xres, x_rot, ymax, y_rot, yres = ds.GetGeoTransform()
geo_transform = (xmin, xres, x_rot, ymax, y_rot, yres)

epsg_code = 'EPSG:4326'
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(int(epsg_code.split(':')[1]))
ds.SetProjection(target_srs.ExportToWkt())

ds.SetGeoTransform(geo_transform)
ds.FlushCache()
ds = None
```

将一个已经存在“WGS 84 / UTM zone 48N”投影的tif文件（test.tif）从本地打开，将投影修改为EPSG 4326（经纬度）并输出为新的tif文件（输出文件test_4326.tif）。

以上就是使用Python和GDAL给图片加坐标系的实现思路，示例代码和相关说明。