用python gdal库来读取nc文件进行波段计算，添加坐标参考系，导出成tif

有一个项目是要求对数据进行修正，我这边有原始的nc栅格数据文件以及一个修正文件，要做的事情就是两个文件一一相加，并导出结果。下边是我用python做的实力

'''
Created on 2018年1月11日

@author: cyemeng
'''
import gdal, osr, ogr, os
from os.path import split

driver_nc = gdal.GetDriverByName('netCDF')

def array2raster(newRasterfn,rasterOrigin,xstep,ystep,array):
    """
    newRasterfn：输出tif路径
    rasterOrigin:原始栅格数据路径
    xstep：x方向栅格数据的宽度
    ystep:y方向栅格数据的宽度
    array:计算后的栅格数据
    """
    cols = array.shape[1]
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn,cols, rows, 1,gdal.GDT_Float32)
    outRaster.SetGeoTransform((originX, xstep, 0, originY, 0, ystep))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

path = "f:/station/history/"
for dirpath,dirnames,filenames in os.walk(path):
    for filename in filenames:
        filepath=os.path.join(dirpath,filename) 
        origin_nc = filepath.replace('\\','/')  # "f:/station/history/1980/01/01/TEM.nc" 
        print(origin_nc) 
        
        father_path = os.path.abspath(os.path.dirname(origin_nc)+os.path.sep+"..")  #父路径
        corr_file = "f:/station/dT_correction/%d.nc"%int(split(father_path)[-1]) #修正文件路径
        ds_arr = gdal.Open( origin_nc ).ReadAsArray() #读取原始nc文件生成array
        correct_arr = gdal.Open( corr_file ).ReadAsArray() #读取修正nc文件生成array

        array = ds_arr+ correct_arr #两个arr相加
        reversed_arr = array[::-1] # reverse array so the tif looks like the array
        
        dst_filename = origin_nc+".tif" #输出路径
        xstep=0.02 #网格间距0.02
        ystep=0.02 #网格间距0.02
        
        #调用array 到raster
        array2raster(dst_filename, [74.99,14.99],xstep,ystep, reversed_arr)