python批量赋予图片坐标系及根据坐标合并图片

• 参考文章
  • gdal给栅格图像加地理坐标
  • 使用Python和GDAL给图片加其原有坐标系，进行坐标投影
  • Python遥感影像拼接

1. python批量赋予图片坐标系

1.1 环境条件准备

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

1.2. 查找投影坐标系

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

1.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_Byte
      elif 'int16' in img_array.dtype.name:
          datatype = gdal.GDT_UInt16
      else:
          datatype = gdal.GDT_Float32
  
      if len(img_array.shape) == 3:
          img_bands, im_height, im_width = img_array.shape
      else:
          img_bands, (im_height, im_width) = 1, img_array.shape
      driver = 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.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 gpd
  
  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_Byte
      elif 'int16' in img_array.dtype.name:
          datatype = gdal.GDT_UInt16
      else:
          datatype = gdal.GDT_Float32
  
      if len(img_array.shape) == 3:
          img_bands, im_height, im_width = img_array.shape
      else:
          img_bands, (im_height, im_width) = 1, img_array.shape
      driver = 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+1
      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)
      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)
  
  

2.python批量根据坐标合并图片并重新赋予坐标

• 读取文件名

  dirpath = " "   # 图片文件绝对路径
  save_path= os.path.join(dirpath, 'merge_tif')  # 合并图片文件保存路径
  if not os.path.exists(save_path):
          os.makedirs(save_path)
  out_path = os.path.join(save_path, "merge.tif")
  coordinate_path=os.path.join(save_path, "coordinate_merge.tif") # 赋予合并坐标
  tif_file = glob.glob(os.path.join(dirpath, "*.tif"))  # 读取所有文件名
  

• 将文件加入列表

  src_files_to_mosaic = []
  for tif_f in tif_file:
      src = rasterio.open(tif_f)
      src_files_to_mosaic.append(src)
  

• 合并及文件保存
  • tif文件为灰度图片，所以我需要mosaic.reshape(mosaic.shape[1],mosaic.shape[2])

  mosaic, out_trans = merge(src_files_to_mosaic)
  Image.fromarray(mosaic.reshape(mosaic.shape[1],mosaic.shape[2]),'L').save(out_path) 
  

• 重新赋予图片坐标

  def def_geoCoordSys(read_path, save_name,img_transf, img_proj):
          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_Byte
          elif 'int16' in img_array.dtype.name:
              datatype = gdal.GDT_UInt16
          else:
              datatype = gdal.GDT_Float32
  
          if len(img_array.shape) == 3:
              img_bands, im_height, im_width = img_array.shape
          else:
              img_bands, (im_height, im_width) = 1, img_array.shape
          driver = 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!')
  img_pos_transf=((out_trans[2],out_trans[0],0.0,out_trans[5],0.0,out_trans[4]))
  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(out_path,coordinate_path, img_pos_transf, img_pos_proj)
  

• 完整代码

  import os
  import glob
  import rasterio
  from osgeo import gdal
  import PIL.Image as Image
  from rasterio.merge import merge
  def def_geoCoordSys(read_path, save_name,img_transf, img_proj):
          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_Byte
          elif 'int16' in img_array.dtype.name:
              datatype = gdal.GDT_UInt16
          else:
              datatype = gdal.GDT_Float32
  
          if len(img_array.shape) == 3:
              img_bands, im_height, im_width = img_array.shape
          else:
              img_bands, (im_height, im_width) = 1, img_array.shape
          driver = 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!')
  dirpath = ""
  save_path= os.path.join(dirpath, 'merge_tif')
  if not os.path.exists(save_path):
          os.makedirs(save_path)
  out_path = os.path.join(save_path, "merge.tif")
  coordinate_path=os.path.join(save_path, "coordinate_merge.tif")
  tif_file = glob.glob(os.path.join(dirpath, "*.tif"))
  src_files_to_mosaic = []
  for tif_f in tif_file:
  	src = rasterio.open(tif_f)
  	src_files_to_mosaic.append(src)
  mosaic, out_trans = merge(src_files_to_mosaic)
  Image.fromarray(mosaic.reshape(mosaic.shape[1],mosaic.shape[2]),'L').save(out_path) 
  img_pos_transf=((out_trans[2],out_trans[0],0.0,out_trans[5],0.0,out_trans[4]))
  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(out_path,coordinate_path, img_pos_transf, img_pos_proj)