#mosica 兩張圖像
import os, sys, gdal
from gdalconst import *
os.chdir('c:/temp/****')#改變文件夾路徑
# 注冊gdal(required)
gdal.AllRegister()

# 讀入第一幅圖像
ds1 = gdal.Open('**.img')
band1 = ds1.GetRasterBand(1)
rows1 = ds1.RasterYSize
cols1 = ds1.RasterXSize

# 獲取圖像角點(diǎn)坐標(biāo)
transform1 = ds1.GetGeoTransform()
minX1 = transform1[0]
maxY1 = transform1[3]
pixelWidth1 = transform1[1]
pixelHeight1 = transform1[5]#是負(fù)值（important）
maxX1 = minX1 + (cols1 * pixelWidth1)
minY1 = maxY1 + (rows1 * pixelHeight1)

# 讀入第二幅圖像
ds2 = gdal.Open('**.img')
band2 = ds2.GetRasterBand(1)
rows2 = ds2.RasterYSize
cols2 = ds2.RasterXSize

# 獲取圖像角點(diǎn)坐標(biāo)
transform2 = ds2.GetGeoTransform()
minX2 = transform2[0]
maxY2 = transform2[3]
pixelWidth2 = transform2[1]
pixelHeight2 = transform2[5]
maxX2 = minX2 + (cols2 * pixelWidth2)
minY2 = maxY2 + (rows2 * pixelHeight2)

# 獲取輸出圖像坐標(biāo)
minX = min(minX1, minX2)
maxX = max(maxX1, maxX2)
minY = min(minY1, minY2)
maxY = max(maxY1, maxY2)

#獲取輸出圖像的行與列
cols = int((maxX - minX) / pixelWidth1)
rows = int((maxY - minY) / abs(pixelHeight1))

# 計(jì)算圖1左上角的偏移值（在輸出圖像中）
xOffset1 = int((minX1 - minX) / pixelWidth1)
yOffset1 = int((maxY1 - maxY) / pixelHeight1)

# 計(jì)算圖2左上角的偏移值（在輸出圖像中）
xOffset2 = int((minX2 - minX) / pixelWidth1)
yOffset2 = int((maxY2 - maxY) / pixelHeight1)

# 創(chuàng)建一個(gè)輸出圖像
driver = ds1.GetDriver()
dsOut = driver.Create('mosiac.img', cols, rows, 1, band1.DataType)#1是bands，默認(rèn)
bandOut = dsOut.GetRasterBand(1)

# 讀圖1的數(shù)據(jù)并將其寫到輸出圖像中
data1 = band1.ReadAsArray(0, 0, cols1, rows1)
bandOut.WriteArray(data1, xOffset1, yOffset1)

#讀圖2的數(shù)據(jù)并將其寫到輸出圖像中
data2 = band2.ReadAsArray(0, 0, cols2, rows2)
bandOut.WriteArray(data2, xOffset2, yOffset2)
''' 寫圖像步驟'''
# 統(tǒng)計(jì)數(shù)據(jù)
bandOut.FlushCache()#刷新磁盤
stats = bandOut.GetStatistics(0, 1)#第一個(gè)參數(shù)是1的話，是基于金字塔統(tǒng)計(jì)，第二個(gè)
#第二個(gè)參數(shù)是1的話：整幅圖像重度，不需要統(tǒng)計(jì)
# 設(shè)置輸出圖像的幾何信息和投影信息
geotransform = [minX, pixelWidth1, 0, maxY, 0, pixelHeight1]
dsOut.SetGeoTransform(geotransform)
dsOut.SetProjection(ds1.GetProjection())

# 建立輸出圖像的金字塔
gdal.SetConfigOption('HFA_USE_RRD', 'YES')
dsOut.BuildOverviews(overviewlist=[2,4,8,16])#4層

 def tif_read(tifpath, bandnum):
  """
  Use GDAL to read data and transform them into arrays.
  :param tifpath:tif文件的路徑
  :param bandnum:需要讀取的波段
  :return:該波段的數(shù)據(jù)，narray格式。len(narray)是行數(shù)，len(narray[0])列數(shù)
  """
  image = gdal.Open(tifpath) # 打開該圖像
  if image == None:
   print(tifpath + "該tif不能打開!")
   return
  lie = image.RasterXSize # 柵格矩陣的列數(shù)
  hang = image.RasterYSize # 柵格矩陣的行數(shù)
  im_bands = image.RasterCount # 波段數(shù)
  im_proj = image.GetProjection() # 獲取投影信息
  im_geotrans = image.GetGeoTransform() # 仿射矩陣
  print('該tif:{}個(gè)行，{}個(gè)列，{}層波段, 取出第{}層.'.format(hang, lie, im_bands, bandnum))
  band = image.GetRasterBand(bandnum) # Get the information of band num.
  band_array = band.ReadAsArray(0,0,lie,hang) # Getting data from zeroth rows and 0 columns
  # band_df = pd.DataFrame(band_array)
  del image # 減少冗余
  return band_array, im_proj, im_geotrans

def tif_write(self, filename, im_data, im_proj, im_geotrans):
  """
  gdal數(shù)據(jù)類型包括
  gdal.GDT_Byte,
  gdal.GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
  gdal.GDT_Float32, gdal.GDT_Float64
  :param filename: 存出文件名
  :param im_data: 輸入數(shù)據(jù)
  :param im_proj: 投影信息
  :param im_geotrans: 放射變換信息
  :return: 0 
  """
  if 'int8' in im_data.dtype.name: # 判斷柵格數(shù)據(jù)的數(shù)據(jù)類型
   datatype = gdal.GDT_Byte
  elif 'int16' in im_data.dtype.name:
   datatype = gdal.GDT_UInt16
  else:
   datatype = gdal.GDT_Float32
  # 判讀數(shù)組維數(shù)
  if len(im_data.shape) == 3:
   im_bands, im_height, im_width = im_data.shape
  else:
   im_bands, (im_height, im_width) = 1,im_data.shape # 多維或1.2維
  #創(chuàng)建文件
  driver = gdal.GetDriverByName("GTiff")   #數(shù)據(jù)類型必須有，因?yàn)橐?jì)算需要多大內(nèi)存空間
  dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
  dataset.SetGeoTransform(im_geotrans)    #寫入仿射變換參數(shù)
  dataset.SetProjection(im_proj)     #寫入投影
  if im_bands == 1:
   dataset.GetRasterBand(1).WriteArray(im_data) #寫入數(shù)組數(shù)據(jù)
  else:
   for i in range(im_bands):
    dataset.GetRasterBand(i+1).WriteArray(im_data[i])
  del dataset

# 這個(gè)展示的效果并不是太好，當(dāng)做示意圖用
 def tif_display(self,im_data):
  """
  :param im_data: 影像數(shù)據(jù),narray
  :return: 展出影像
  """
  # plt.imshow(im_data,'gray') # 必須規(guī)定為顯示的為什么圖像
  plt.imshow(im_data) # 必須規(guī)定為顯示的為什么圖像
  plt.xticks([]), plt.yticks([]) # 隱藏坐標(biāo)線
  plt.show() # 顯示出來，不要也可以，但是一般都要了