• 参考文章
  • 文章1
  • 文章2

1. 环境条件准备

• 相关坐标为：qgis生成的网格文件shpfile文件，获取的相关坐标
  • 每个网格都有相关坐标信息
    
• 赋予的图片，使用网格文件gird_256_4490.shp切割的图片
• 赋予图片坐标系原因
  • 使用qgis切割带坐标的文件后，保存.png格式，没有了坐标信息，赋予切割图片坐标信息

2. 查找投影坐标系

• 官网链接 ：https://spatialreference.org/ref/
• 输入自己要转的坐标系名称
  • 我的是China Geodetic Coordinate System 2000
  • 点击第一个链接
    
• 点击OGC WKT
  
• 复制文本【全部复制】
  • 将文本赋予变量img_pos_proj
    

3. 相关代码

• 图片加坐标系函数

  def def_geoCoordSys(read_path, save_name, img_transf, img_proj):""":param read_path: 不带坐标的图像:param save_name: 保存带有坐标的图像:param img_transf: 投影相关坐标:param img_proj: 投影类型:return:"""array_dataset = gdal.Open(read_path)img_array = array_dataset.ReadAsArray(0, 0, array_dataset.RasterXSize, array_dataset.RasterYSize)if 'int8' in img_array.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in img_array.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(img_array.shape) == 3:img_bands, im_height, im_width = img_array.shapeelse:img_bands, (im_height, im_width) = 1, img_array.shapedriver = gdal.GetDriverByName("GTiff")  # 创建文件驱动dataset = driver.Create(save_name, im_width, im_height, img_bands, datatype)dataset.SetGeoTransform(img_transf)  # 写入仿射变换参数dataset.SetProjection(img_proj)  # 写入投影# 写入影像数据if img_bands == 1:dataset.GetRasterBand(1).WriteArray(img_array)else:for i in range(img_bands):dataset.GetRasterBand(i + 1).WriteArray(img_array[i])print(read_path, 'geoCoordSys get!')
  

• 获取shpfile文件相关坐标
  

  shpdata = gpd.read_file('data/grid_256_4490.shp',encoding='utf-8')
  for j in range(0, len(shpdata)):  # 遍历网格数据geo = shpdata.geometry[j] # 获取坐标属性# ((120.78513142722576 29.253768371972182, 120.78516670362751 29.258387207528127, 120.79043393150636 29.258356166298398, 120.79039841851554 29.25373733658219, 120.78513142722576 29.253768371972182))# 坐标位置为：0：左上，1：左下，2：右下，3：右上，4：左上minX,minY=geo.__geo_interface__['coordinates'][0][1]maxX,maxY=geo.__geo_interface__['coordinates'][0][3]numX,numY=256,256  # 图片x轴像素，y轴像素img_pos_transf=(minX,(maxX-minX)/256 , 0.0, minY, 0.0,-(minY-maxY)/256)img_pos_proj = 'GEOGCS["China Geodetic Coordinate System 2000 ",DATUM["China 2000",SPHEROID["CGCS2000",6378137,298.257222101,AUTHORITY["EPSG","1024"]],AUTHORITY["EPSG","6610"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4490"]]'def_geoCoordSys(read_path,save_name, img_pos_transf, img_pos_proj)
  

• 完整代码

  import os
  from osgeo import gdal
  import geopandas as gpddef def_geoCoordSys(read_path, save_name, img_transf, img_proj):""":param read_path: 不带坐标的图像:param save_name: 保存带有坐标的图像:param img_transf: 投影相关坐标:param img_proj: 投影类型:return:"""array_dataset = gdal.Open(read_path)img_array = array_dataset.ReadAsArray(0, 0, array_dataset.RasterXSize, array_dataset.RasterYSize)if 'int8' in img_array.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in img_array.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(img_array.shape) == 3:img_bands, im_height, im_width = img_array.shapeelse:img_bands, (im_height, im_width) = 1, img_array.shapedriver = gdal.GetDriverByName("GTiff")  # 创建文件驱动dataset = driver.Create(save_name, im_width, im_height, img_bands, datatype)dataset.SetGeoTransform(img_transf)  # 写入仿射变换参数dataset.SetProjection(img_proj)  # 写入投影# 写入影像数据if img_bands == 1:dataset.GetRasterBand(1).WriteArray(img_array)else:for i in range(img_bands):dataset.GetRasterBand(i + 1).WriteArray(img_array[i])print(read_path, 'geoCoordSys get!')label_path=''  # 无坐标图片路径
  label_tif_path=''  # 保存有坐标图片路径
  if not os.path.exists(label_tif_path):os.makedirs(label_tif_path)
  data = os.walk(label_path)
  file_list=[file_list for path, folder_list, file_list in data][0]
  shpdata = gpd.read_file('grid.shp',encoding='utf-8')
  for j in range(0, len(shpdata)):geo = shpdata.geometry[j]id_name=j+1minX,minY=geo.__geo_interface__['coordinates'][0][1]maxX,maxY=geo.__geo_interface__['coordinates'][0][3]numX,numY=256,256  # 图片x轴像素，y轴像素img_pos_transf=(minX,(maxX-minX)/256 , 0.0, minY, 0.0,-(minY-maxY)/256)read_path=os.path.join(label_path,f'{id_name}.png')  # 无坐标图片路径·save_name=os.path.join(label_tif_path,f'{id_name}.tif') # 保存有坐标图片路径save_name=save_name.replace('png','tif')img_pos_proj = 'GEOGCS["China Geodetic Coordinate System 2000 ",DATUM["China 2000",SPHEROID["CGCS2000",6378137,298.257222101,AUTHORITY["EPSG","1024"]],AUTHORITY["EPSG","6610"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4490"]]'def_geoCoordSys(read_path,save_name, img_pos_transf, img_pos_proj)