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