PythonGDAL给图像设置投影坐标系
网站上下载的部分遥感数据只有地理坐标，没有投影坐标，本文章尝试给“只有地理坐标的图片”添加投影坐标。
原始数据：山东地区ASTER GDEM30m数据，该数据只有WGS84地理坐标
目标数据：给该数据添加UTM-50N
1.确定你的目标投影坐标信息，就可以通过下述代码实现。
（1）首先你得知道想添加的目标投影信息，可以通过网上查询的方式。
此网站查的比较费劲，我反正没在这找到：https://developers.arcgis.com/javascript/3/jshelp/pcs.htm，
该网站比较方便，可以搜索：https://spatialreference.org/ref/
我的目标投影坐标是WGS84-UTM zone 50N，https://spatialreference.org/ref/https://spatialreference.org/ref/epsg/32650/html/，见下图：
（2）上代码
from osgeo import gdal,osr,gdalnumeric
from osgeo.gdalconst import *
from pyproj import CRS
# 1.获取原数据信息
# 该数据只有地理坐标WGS84
ds = gdal.Open('D:\ProfessionalProfile\CoordinateSystemConversion\shandongDEMenvi_10.tif')
im_geotrans = ds.GetGeoTransform() #仿射矩阵信息
im_proj = ds.GetProjection() #地图投影信息
# print(im_geotrans)
# print(im_proj)
im_width = ds.RasterXSize # 栅格矩阵的列数
im_height = ds.RasterYSize # 栅格矩阵的行数
im_bands = ds.RasterCount
ds_array = ds.ReadAsArray(0, 0, im_width, im_height) # 获取原数据信息，包括数据类型int16，维度，数组等信息
# # 设置数据类型(原图像有负值)
datatype = gdal.GDT_Float32
# 2.原图像的仿射变换矩阵参数，即im_geotrans
img_transf = (114.79763889, 0.00027777777778, 0.0, 38.21347222, 0.0, -0.00027777777778)
# # 网站查询的WGS84-UTM50N坐标信息https://spatialreference.org/ref/epsg/32650/html/
img_proj = '''PROJCS["WGS 84 / UTM zone 50N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",117],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
AUTHORITY["EPSG","32650"],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]]'''
# 3.设置新文件及各项参数
filename = 'D:\ProfessionalProfile\CoordinateSystemConversion\AproUTM50fInternet.tif'
driver = gdal.GetDriverByName("GTiff") # 创建文件驱动
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(img_transf) # 写入仿射变换参数
dataset.SetProjection(img_proj) # 写入投影
# 写入影像数据
dataset.GetRasterBand(1).WriteArray(ds_array)
del dataset
（3）最后结果
2.在arcgis中定义投影
（1）Arcmap中操作
（2）结果（在ENVI中显示）
（3）对比
有一点儿差别：Python中使用的数据类型是浮点型；原数据是int16。
之前在Python中将数据类型设置为gdal.GDT_Int16时候，原数据中小于0的值都在输出影像中都成了0。没太懂这是什么原因，如果有熟悉的大佬，希望能指点一二。
3.arcgis中转换
如果有目标影像的投影信息的话：只要在ARCGIS中先加载“有目标投影信息的数据”再加载“要投影的数据”，然后再重新导出“要投影的数据”（选择当前投影）就可以了。
4.python中添加
（1）如果有目标影像的投影信息的话，也可以在Python中给“缺失投影信心的数据”添加投影信息。
from osgeo import gdal,osr,gdalnumeric
from osgeo.gdalconst import *
# 数据1是landsat8影像，投影是WGS84-UTM47
ds1 = gdal.Open('D:\ProfessionalProfile\Landsat2017\LC08_L1TP_134036_20170808_20170813_01_T1\LC08_L1TP_134036_20170808_20170813_01_T1_B1.tif')
# ds1_geotrans = ds1.GetGeoTransform()
# 数据2是山东地区的DEM数据，只有地理坐标WGS84，没有投影坐标
ds2 = gdal.Open('D:\ProfessionalProfile\CoordinateSystemConversion\shandongDEMenvi_10.tif',GA_Update)
ds2.SetProjection(ds1.GetProjection())
ds2.SetGeoTransform(ds1.GetGeoTransform())
（2）结果如下：
当然，将山东的DEM投影为UTM-47肯定不对，我这里只不过是尝试一下，看看这样做可不可行。
该山东地区DEM数据的经度范围在114.79763889E——122.70013889E之间，个人感觉比较适合的投影带是UTM-50N。
（3）对比一下UTM-47N和UTM-50N
参考：
地理坐标系WKID：https://developers.arcgis.com/javascript/3/jshelp/gcs.htm
投影坐标系WKID：https://developers.arcgis.com/javascript/3/jshelp/pcs.htm
https://spatialreference.org/ref/
<>
<>