python利用gdal读取、写出tif格式的遥感卫星影像，包含超大数据量的读写

gdal库特别强大，可以很方便的读写带有地理参考的影像数据

1. 使用gdal读取tif图像，读入数组data中：

from osgeo import gdal,osr
def readTif():
	#输入路径地址
	tifpath=r"D:\data.tif"
    #gdal打开影像，成为dataset数据集（这一步没有放入内存中）
    dataset=gdal.Open(tiffile,gdal.GA_ReadOnly)
    #读取仿射参数，仿射参数有详细解析
    geo_transform=dataset.GetGeoTransform()
    #获取投影
    proj = dataset.GetProjection()
    #获取空间参考
    srs = osr.SpatialReference(proj) 
    #获取坐标系的epsg编号
    epsg = srs.GetAttrValue('AUTHORITY',1)   
    
    #当图像仅有一个波段时读取会简单一些，当多波段图像时要依次读取每个波段
    if dataset.RasterCount==1:
	    #按照数组的方式读取数据（单波段）
	    data=dataset.ReadAsArray()
	else:
		data=[]
	    dataarray=dataset.ReadAsArray()
	    for i in range(1,dataset.RasterCount+1) #索引从1开始
	    	band=dataset.GstRasterBand(k)
	    	data.append(band.ReadAsArray())
	    #列表转数组
	    data=np.dstack(data)
	    
	return data

图像中除了数据矩阵（数组）外，最重要的就是仿射参数，即geo_treansform

geo_treansform(0) 左上角像素的x坐标（左上角经度）
geo_treansform(1) 东西向像素分辨率/像素宽（一个像元的经度值，通常为正值，例如哨兵2数据为：10）
geo_treansform(2) 行旋转（通常为零）
geo_treansform(3) 左上角像素的y坐标（左上角纬度）
geo_treansform(4) 列旋转（通常为零）。
geo_treansform(5) 南北向像素分辨率/像素高（一个像元的纬度值，对于上北下南图像为负值，哨兵2数据为：-10）

因为初始点默认为左上角，geo_treansform(0)和geo_treansform(3)给的也是左上角坐标，推算的时候也是以此为起点，例如：

假如一景哨兵2图像，左上角点经度10000，纬度15000，共有20行，30列
则最有一个像元（右下角）的经度为10000+30*10=10300，纬度为15000+（-20*10）=14800

2. 使用gdal的写入tif图像：

##geo_transform、porjection由读入图像确定，datadtype可以按照输出要求设定
def geotiffwrite(data,geo_transform,projection,datatype):
	outputPath=r"D:\result.tif"
	
	#定义要输出的图像命成、行列数、波段、数据类型
    driver = gdal.GetDriverByName("GTiff")
    rows,cols,bands=data.shape
    dataset=driver.Create(outputPath,cols,rows,bands,gdal.GDT_Float32)
    #定义仿射参数、投影等
    dataset.SetGeoTransform(geo_transform)
    dataset.SetProjection(projection)
    #写入数组
    if bands == 1:
        dataset.GetRasterBand(1).WriteArray(data)
    else:
        for i in range(bands):
            dataset.GetRasterBand(i+1).WriteArray(data[:,:,i])

3. 当影像超大时（超内存）：
当用到超大数据，如100GB的单景数据，计算机内存受限时，可以采用gdal提供的部分读取和部分写入的方法，每次只将少量的数据读入内存中，可以解决这个问题。
（1）读入的问题

>>> from osgeo import gdal
>>> dataset = gdal.Open("D:\data.tif")
>>> help(dataset.ReadRaster)
>Help on method ReadRaster in module osgeo.gdal:

ReadRaster(xoff=0, yoff=0, xsize=None, ysize=None, buf_xsize=None, buf_ysize=None, buf_type=None, band_list=None, buf_pixel_space=None, buf_line_space=None, buf_band_space=None, resample_alg=0, callback=None, callback_data=None) method of osgeo.gdal.Dataset instance

xoff,yoff ：指定想要读取的部分原点位置在整张图像中距离全图原点的位置(以像元为单位）。
xsize,ysize ： 指定要读取部分图像的矩形的长和宽（以像元为单位）。
buf_xsize，buf_ysize ：可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高，gdal将帮你缩放到这个大小。
buf_type ：可以对读出的数据的类型进行转换（比如原图数据类型是short，你要把它们缩小成byte）。
band_list ：适应多波段的情况。可以指定要读取的波段。

也就是说对于一个超大数据，可以只读其中一部分，甚至只读一个点，再结合循环和部分写出的方法就能实现：将大图像分批次读入，计算处理，并分批次写出成一个完整数据。

（2）写出的问题

WriteArray()除了可以写出一整景数据，还可以只写一小部分。
假设最终要写出的影像是一个10000列、20000行的超大影像，如果按照整景写出的方法，

dataset.GetRasterBand(1).WriteArray(data)

data过大，超出内存，无法写入成功！
可以采用依次写入的方法：

##geo_transform、porjection由读入图像确定，datadtype可以按照输出要求设定
def geotiffwrite(data,geo_transform,projection,datatype):
	driver = gdal.GetDriverByName("GTiff")
	##定义一个最终要写出的影像是一个10000列、20000行的超大影像
	newdataset = driver.Create(outputFile,10000,20000,1,gdal.GDT_Float32)
	newdataset.SetProjection(projection())
	newdataset.SetGeoTransform(c_geotrans)
	out_band = newdataset.GetRasterBand(1)
	block = np.random.rand(100, 100)
	#这里对于这个大影像，只将左上角的100*100的区域写入了数据，其他地方则没有写。
	out_band.WriteArray(block,xoff=0,yoff=0)

xoff和yoff规定这个block写在大图像的什么位置（0，0）就代表在左上角写入。

这样结合读入、写出、循环的方式，就可以实现一个超大影像的完整处理过程，完整代码如下：

以超大影像的ndvi计算为例：

def NDVI(inputfile,outputfile):
	#gdal打开数据
    dataset=gdal.Open(inputfile)
    rows=dataset.RasterYSize
    cols=dataset.RasterXSize
    bands=dataset.RasterCount
  	
  	#新建影像，用来保存ndvi
	driver = gdal.GetDriverByName("GTiff")
    newdataset = driver.Create(outputFile,cols,rows,1,gdal.GDT_Float32)
    newdataset.SetProjection(dataset.GetProjection())
    newdataset.SetGeoTransform(dataset.GetGeoTransform())
    out_band = newdataset.GetRasterBand(1)
	
	#确定一个切片的大小，可以根据计算机的可用内存确定
    pixelnum=2048
    
    #计算需要切成多少块
    xnum = ceil(cols/pixelnum)
    ynum = ceil(rows/pixelnum)
    
    #切块开始
    for i in range(ynum):
        if ynum == 1:               #垂直方向不需要切的时候   
            yoff = 0
            ysize = rows
        else:                          
            if i == 0:              #垂直方向需要切，最上面一条
                yoff = 0
                ysize = pixelnum       
            elif i == ynum-1:       ##垂直方向需要切，最下面一条
                yoff = i*pixelnum
                ysize = rows - i*pixelnum
            else:					##垂直方向需要切，中间的条
                yoff = i*pixelnum
                ysize = pixelnum
        for j in range(xnum):      ##水平方向不需要切
            if xnum == 1:                     
                xoff = 0
                xsize = cols
            else:                         
                if j == 0:         ##水平方向需要切，最左边一条
                    xoff=0
                    xsize = pixelnum       
                elif j == xnum-1:  ##水平方向需要切，最右边一条
                    xoff = j*pixelnum
                    xsize = cols - j*pixelnum
                else:				 ##水平方向需要切，中间的条
                    xoff = j*pixelnum
                    xsize = pixelnum 
            
            #读取当前block的数据
            data=[]
            for k in range(1,bands+1):
                band=dataset.GetRasterBand(k)
                data.append(band.ReadAsArray(xoff,yoff,xsize,ysize))
            data=np.dstack(data)
            #计算当前block的数据对应的ndvi
            ndvi=(data[:,:,3]-data[:,:,2])/(data[:,:,3]+data[:,:,2])
            
            #写出当前block的数据计算结果
           	out_band.WriteArray(ndvi,xoff,yoff) 
 

另外切片的大小该怎样确定呢？我是按照自己的可用内存实验出来的，有没有一种方法可以确定：内存和其能处理的数组大小之间的关系，如果有懂的，希望可以不吝赐教！
附：简单的可用内存大致确定切块大小

##获取计算机可用内存，并按照可用，确定切片大小
 import psutil

 available_memory = psutil.virtual_memory().available
 available_memory_gb = available_memory / (1024 ** 3)  #可用内存
 gbFor2=int(log(available_memory_gb,2))  ##是2的多少次方  
 pixelnum=gbFor2*2048
if gbFor2==0:
   pixelnum=512