gdal合成葵花八号卫星nc格式数据为envi格式、tiff格式图像

背景

在看文献时发现作者使用Satpy合成葵花八号数据为三通道图像，但是经过自己尝试后发现Satpy只能支持葵花八号标准数据HSD

该类型数据只有最近三十天，而我需要使用一整年的数据，故选用的是/jma/netcdf文件夹下的nc格式数据，然而Satpy并不能识别这种格式，所以经过查找资料，选用gdal可以成功融合。

1.数据下载

可见上一文章葵花八号卫星多任务数据下载

2.python包准备

需要的包和版本如下

1.gdal
gdal包下载需要使用oseg

2.netcdf

3.numpy

要维护好环境，不然随时会崩，缺少关键包导致跑不起来，特别是numpy和matplotlib经常导致版本冲突。

3.gdal函数分析

driver = gdal.GetDriverByName('ENVI')

创建一个gdal驱动，名字可选ENVI、GTiff、JPEG、PNG、Shp等

out_tif = driver.Create(fileName,width,height,bands,eType)

fileName:文件地址
width:图像宽度
height:图像高度
bands:通道数
eType:图像格式（float32、uint32、int16等）

out_tif.SetProjection(srs.ExportToWkt())

设置投影

out_tif.GetRasterBand(1).WriteArray(band1)

通道读取并写入

4.合成关键代码

'''
1-图像融合
'''
import calendar
from osgeo  import gdal,osr
import numpy as np
import netCDF4 as nc
import os

def layerStack(file_path):
    file_path = file_path
    f_list=os.listdir(file_path)
    images=[]
    for i in f_list:
        # os.path.splitext():分离文件名与扩展名
        if os.path.splitext(i)[1] == '.nc':
            images.append(i)
    count = 1
   
    for dataName in images:
        # 通过netcdf读取数据
        d = nc.Dataset(file_path + '/'+dataName)
        # 选取通道 band1、band2、band3分别是R、G、B
        band1 = np.asarray(d['albedo_01'][:])
        band2 = np.asarray(d['albedo_02'][:])
        band3 = np.asarray(d['albedo_03'][:])
        # 获取地理信息
        lat = d.variables['latitude'][:]
        lon = d.variables['longitude'][:]
        latMin, latMax, lonMin, lonMax = lat.min(), lat.max(), lon.min(), lon.max()
        # 分辨率
        lat_Res = (latMax - latMin) / (lat.shape[0] - 1)
        lon_Res = (lonMax - lonMin) / (lon.shape[0] - 1)

        # 保存的路径和名称，可选GTiff、ENVI
        # driver = gdal.GetDriverByName('GTiff')
        driver = gdal.GetDriverByName('ENVI')
        split=dataName.split('_')
        #合成图像保存路径
        out_path=file_path+'/dat/'
        out_tif_name = str(count-1)+'.dat'
        out_tif = driver.Create(out_path+out_tif_name, lon.shape[0], lat.shape[0], bands=3, eType=gdal.GDT_Float32)

        # 输出结果
        geotransform = (lonMin, lon_Res, 0, latMax, 0, -lat_Res)
        out_tif.SetGeoTransform(geotransform)
        srs = osr.SpatialReference()
        srs.SetWellKnownGeogCS('WGS84')
        out_tif.SetProjection(srs.ExportToWkt())

        out_tif.GetRasterBand(1).WriteArray(band1)

        out_tif.GetRasterBand(2).WriteArray(band2)

        out_tif.GetRasterBand(3).WriteArray(band3)

        print('融合完成第',count)
        count=count+1
        out_tif.FlushCache()
        out_tif = None