遥感影像的镶嵌（初试）
在做PM2.5的机器学习的时候，不同的aod数据的利用率问题就显得十分重要，能多要一个数据都显得弥足珍贵，而几个类型的数据中，modis3km的数据较为杂乱，因为他是5分钟的采样时间，也就是一天会产生大于24*60/5的数据量，具体怎么解释我现在也迷迷糊糊，过境数据每次位置都有不同，一个矩形的范围内，最多可以有好几个影像的范围。
好比如此。
所以数据的镶嵌，拼接就比较重要，gdal有一个镶嵌的merged,py文件，调用起来暂时没想通怎么用python调用，应该比较好用才是。参照网上的一些方法，利用行列号把数据写入新的栅格，但是出现的偏差很大，因为缺少了投影计算的步骤，所以目前的代码应该还不能用。后期新增这部分应该也算是创新吧。
上代码：
import os,sys,gdal
from gdalconst import *
import glob
def get_extent(fn):
ds=gdal.Open(fn)
height=ds.RasterYSize
weight=ds.RasterXSize
gt=ds.GetGeoTransform()
minx=gt[0]
maxy=gt[3]
maxx=gt[0]+gt[1]*weight
miny=gt[3]+gt[5]*height
res=gt[1]
return(minx,maxy,maxx,miny,res)
in_files=os.listdir('C:/pytemp/modismosictemp/modis3kre')
minx,maxy,maxx,miny,resb=get_extent('C:/pytemp/modismosictemp/modis3kre/'+in_files[0])
for i in in_files[1:]:
minXa,maxYa,maxXa,minYa,resa=get_extent('C:/pytemp/modismosictemp/modis3kre/'+i)
minX=min(minXa,minx)
maxY=max(maxYa,maxy)
maxX=max(maxXa,maxx)
minY=min(minYa,miny)
resmin=min(resa,resb)
resmax=max(resa,resb)
in_ds=gdal.Open('C:/pytemp/modismosictemp/modis3kre/'+in_files[0])
gt=in_ds.GetGeoTransform()
weight=int(abs(int(maxX-minX))/resmin)
height=int(abs(int(maxY-minY))/resmin)
format=("GTiff")
driver=gdal.GetDriverByName(format)
outDataRaster=driver.Create('D:/b/mo.tif',weight,height,1,gdal.GDT_Int16)####cols=0的话就没法继续看了
gt = [minX, resmax, 0, minY, 0, resmax]
outDataRaster.SetGeoTransform(gt)
outDataRaster.SetProjection(in_ds.GetProjection())
out_band = outDataRaster.GetRasterBand(1)
import numpy as np
raster=np.zeros((height,weight))
in_ds=gdal.Open('C:/pytemp/modismosictemp/modis3kre/'+in_files[0])
a=in_ds.ReadAsArray()
yk,xk=a.shape
x1,y1=int(abs(minx-minX)/resmax),int(abs(maxy-maxY)/resmax)
x2,y2=int(abs(maxx-minX)/resmax),int(abs(miny-maxY)/resmax)
for x in range(xk):
for y in range(yk):
raster[y+min(y1,y2),x+min(x1,x2)]=a[y,x]
in_ds1=gdal.Open('C:/pytemp/modismosictemp/modis3kre/'+in_files[1])
a1=in_ds1.ReadAsArray()
yk1,xk1=a1.shape
x3,y3=int(abs(minXa-minX)/resmax),int(abs(maxYa-maxY)/resmax)
x4,y4=x3+xk1,y3+yk1
for xx in range(xk1):
for yx in range(yk1):
raster[yx+min(y3,y4),xx+min(x3,x4)]=a1[yx,xx]
raster1=np.transpose(np.transpose(raster)[::-1])[::-1][:,::-1]##其中这个部分是数组转置翻转等操作
out_band.WriteArray(raster1)
outDataRaster.FlushCache() ## remove from memory...important
del outDataRaster ## delete the data (not the actual geotiff)
上图的应该镶嵌的形状，下图是代码出来的形状。。。前后有差异，等我更新吧。。。