气象数据Grib格式解析的Python代码和Matlab代码

以.grb/.grb1/.grb2为扩展名的都是气象数据，气象数据中可以存储多个内容，如云量、雪深、气压、风速等内容，或者具有时间序列的云量等。这些文件不可以直接打开成图片，若想直观地查看grib数据，需要读取文件并将其解析出来，保存成tif或者png格式

这几天分别用matlab代码和python代码解析成tif/png格式的图片，并将其插值成任意需要的分辨率。matlab运行起来速度较慢，且容易造成数据量过大停止运行的情况，因此将python重新写了一份。例如：原数据中经纬度分辨率都是0.125°，matlab只能将其插值到0.02°，不能将其插值到0.01°；而python则可以。原数据3.95M，用python插值到0.01分辨率后文件大小为4.85G。可见数据量之大。经纬度的分辨率可以由【MatLab部分】的最后一块代码中生成的.xlsx文件中查看，也可以通过经纬度和对应的行列数计算出来。

【MatLab部分】

grib文件具有数据标签，例如该文件内部的数组格式为data（2,1441,2880），那么data(0)可能是风速，data(1)可能是云量等。下面展示框是展示了一个.grb文件的数据标签，具体代码后面会有。

可以从标签中读到该数据的维度是经度2880 x 维度1441，只有一个时间段，变量中只有Total_cloud_cover_surface一个，即地表云量信息。

netcdf E:/data/NAFP_ECMF_1_FTM-98-GLB-TCC-125X125-1-0-999998-999998-999998-2018110712-0.GRB {dimensions:lon = 2880;lat = 1441;time = 1;variables:float Total_cloud_cover_surface(time=1, lat=1441, lon=2880);:long_name = "Total cloud cover @ Ground or water surface";:units = "(0.-.1)";:missing_value = NaNf; // float:grid_mapping = "LatLon_Projection";:Grib_Variable_Id = "VAR_98-0-128-164_L1";:Grib1_Center = 98; // int:Grib1_Subcenter = 0; // int:Grib1_TableVersion = 128; // int:Grib1_Parameter = 164; // int:Grib1_Parameter_Name = "TCC";:Grib1_Level_Type = 1; // intfloat lat(lat=1441);:units = "degrees_north";:_CoordinateAxisType = "Lat";float lon(lon=2880);:units = "degrees_east";:_CoordinateAxisType = "Lon";int time(time=1);:units = "Hour since 2018-11-07T12:00:00Z";:standard_name = "time";:long_name = "Initialized analysis product for reference time";:_CoordinateAxisType = "Time";:Originating_or_generating_Center = "European Centre for Medium Range Weather Forecasts (ECMWF) (RSMC)";:Originating_or_generating_Subcenter = "0";:Conventions = "CF-1.6";:history = "Read using CDM IOSP Grib1Collection";:featureType = "GRID";:file_format = "GRIB-1";:_CoordSysBuilder = "ucar.nc2.dataset.conv.CF1Convention";
}

1.读取单个grib文件，输出数据标签、数据标签存入txt、将数据以tif图片格式保存

% Matlab codefunction readfilesetup_nctoolbox[filename, pathname] = uigetfile('*.grb', 'choose a GRB file'); %弹出对话框选择.grb文件if isequal(filename,0)    msgbox('you choose nothing');elsepathfile=fullfile(pathname, filename); %获取文件地址全名ds= ncdataset(pathfile);%读文件disp('file info:')%在命令行输出file info:字样label=char(ds.netcdf) %输出数据标签fid=fopen('E:/label.txt','w');%将数据标签存放到txt中fprintf(fid,'%s',label);fclose(fid);GPMData = ds.data(ds.variables{1});%获得数据中的云量值，{1}可能需要根据文件进行改动GPMData = squeeze(GPMData);%去掉长度为1的维度imwrite(mat2gray(GPMData),'E:/cloud.tif');%保存成tif图片end

2.批量读取grib文件，保存成tif格式,将该文件数据标签存为txt文件。读取的是一个文件夹内的所有符合格式的数据，如果文件名的扩展名前有“SOIL”，则不对该文件进行处理。

% Matlab codefunction readfilesetup_nctoolboxfile_path='E:\data\2018110712\';%设置数据文件夹img_list=dir(strcat(file_path,'*.grib1'));%读取file_path文件夹内的所有的文件名，放入img_list中num=length(img_list);%文件的个数str1='E:\dataResult\2018110712\';str3='.tif';str4='.txt';for k=1:num%对所有文件循环img_name=img_list(k).name%获取文件名，假如该文件名为 20180101_ABC_SOIL.GRBname=strsplit(img_name,'.');%将该文件名以“.”分为纯文件名和扩展名,即20180101_ABC_SOIL和GRBa=strsplit(char(name(1)),'_');%再将该文件的纯文件名部分以“_”分隔开，即a=SOILif strcmp(char(a(length(a))),'SOIL')%如果a==SOIL即不对该文件处理continue;elsepathfile=fullfile(file_path, img_name);%获得该文件包含路径的文件名ds= ncdataset(pathfile);%读取该文件label=char(ds.netcdf);%获得数据标签GPMData = ds.data(ds.variables{1});%获取文件数据，数值1需根据文件更改GPMData = squeeze(GPMData);%去除长度1的维度str=[str1,str2,str3];%str为保存的文件名imwrite(mat2gray(GPMData),str);         %将文件保存为tif图像str2=char(name(1));%将该文件数据标签存为txttxtname=[str1,str2,str4];fid=fopen(txtname,'w');fprintf(fid,'%s',label);fclose(fid);endendend

3.批量将文件夹内的所有tif文件转换至具有透明通道的png图片

% Matlab codefunction converse%setup_nctoolboxfilepath='E:\ecGRBfiles\TCC\';img_list=dir(strcat(filepath,'*.tif'));num=length(img_list);type='.png';for k=1:numimg_name=img_list(k).name;filename=strsplit(img_name,'.');filename=char(filename(1));str=[filepath,filename,type];pathfile=fullfile(filepath, img_name);file=imread(pathfile);alpha=mapminmax(file,0,1);alpha=double(alpha);%alpha为透明值imwrite(mat2gray(file),str,'Alpha',alpha);%‘Alpha’参数为设置透明通道end
end

4.将源文件经纬度分辨率均为0.125的插值成0.01，采用三次插值，由于是图像，因此使用interp2()函数实现二维插值

% Matlab codefunction readfilesetup_nctoolbox[filename, pathname] = uigetfile('*.grb', 'choose a GRB file'); %弹出对话框选择.grb文件if isequal(filename,0)    msgbox('you choose nothing');elsepathfile=fullfile(pathname, filename); %获取文件地址全名ds= ncdataset(pathfile);%读文件disp('file info:')%在命令行输出file info:字样%将.grb输出为透明度png        0为完全透明,1为不透明GPMData = ds.data(ds.variables{1});GPMData = squeeze(GPMData);%388*190*384alpha=mapminmax(GPMData,0,1);alpha=double(alpha);imwrite(mat2gray(GPMData),'cloud_origin.png','Alpha',alpha);%插值，原数据分辨率为0.125，经度从0~359.875，纬度从90~-90x_origin=0:0.125:359.875;y_origin=90:-0.125:-90;data=GPMData;%0.05插值[x1,y1]=meshgrid(0:0.05:359.875,90:-0.05:-90);in_data1=interp2(x_origin,y_origin,data,x1,y1,'cubic');%  alpha1=mapminmax(in_data1,0,1);%  alpha1=double(alpha1);%  imwrite(mat2gray(in_data1),'cloud_005.png','Alpha',alpha1);%0.02插值[x2,y2]=meshgrid(0:0.02:359.875,90:-0.02:-90);in_data2=interp2(x1,y1,in_data1,x2,y2,'cubic');%  alpha2=mapminmax(in_data2,0,1);%  alpha2=double(alpha2);%  imwrite(mat2gray(in_data2),'cloud_002.png','Alpha',alpha2);imwrite(mat2gray(in_data2),'cloud_002.tif');% 0.01插值[x3,y3]=meshgrid(0:0.01:359.875,90:-0.01:-90);in_data3=interp2(x2,y2,in_data2,x3,y3,'cubic');alpha3=mapminmax(in_data3,0,1);alpha3=double(alpha3);imwrite(mat2gray(in_data3),'cloud_001.png','Alpha',alpha3);% write data to excelname = {'lontitude','latitude','GPM'};xlswrite(strcat(filename,'.xlsx'), name,1,'A1')xlswrite(strcat(filename,'.xlsx'), lon,1,'A2')xlswrite(strcat(filename,'.xlsx'), lat,1,'B2')xlswrite(strcat(filename,'.xlsx'), GPMData,1,'C2')msgbox('finished!');end

【Python部分】

主要思路是采用GDAL读取GRIB文件，然后用插值函数将0.125分辨率插值成0.05和0.01的。

1.采用GDAL处理单张的GRIB文件，并用scipy模块的interp2d函数进行二维插值，插值为0.05和0.01分辨率

# read .grib1 and .grib2 file and output a .tif image
# completed at 2018/11/20
# writen by Zhou Dengji
import os
import sys
from osgeo import gdal
from gdalconst import *
from scipy import interpolate
import numpy as np
import pylab as plos.chdir('E:/ecGRBfiles/TCC')#设置当前工作目录
driver = gdal.GetDriverByName("GTiff")#设置driver驱动为GTiff格式
gdal.AllRegister()filename = '0.GRB'#文件名
ds = gdal.Open(filename, GA_ReadOnly)#读文件
if ds is None:print('Cannot open this .grb file.')sys.exit(1)cols = ds.RasterXSize#列数
rows = ds.RasterYSize#行数
bands = ds.RasterCount#波段数
data = np.zeros((bands, rows, cols), dtype=np.float64)#存放数据的数组
for i in range(bands):data[i] = ds.GetRasterBand(i + 1).ReadAsArray(0, 0, cols, rows)#逐波段读数据datax = np.arange(0, 360, 0.125)  # 原图的x方向尺寸，从0~360，间隔0.125
datay = np.arange(90, -90.125, -0.125)  # 原图的y方向尺寸，从90~-90，间隔-0.125
dataxx, datayy = np.meshgrid(datax, datay)  # 以datax/datay建立格网############## 0.05间隔的插值
dataxnew005 = np.arange(0, 360, 0.05)  # 插值后的x方向尺寸，从0~360，间隔0.05
dataynew005 = np.arange(90, -90.125, -0.05)  # 插值后的y方向尺寸，从90~-90，间隔-0.05
rows005 = len(dataynew005)  # y方向个数，即行数
cols005 = len(dataxnew005)  # x方向个数，即列数
datanew005 = np.zeros((bands, rows005, cols005), dtype=np.float64)  # 建立新data存放插值后的数据############## 0.01间隔的插值
dataxnew001 = np.arange(0, 360, 0.01)  # 插值后的x方向尺寸，从0~360，间隔0.01
dataynew001 = np.arange(90, -90.125, -0.01)  # 插值后的y方向尺寸，从90~-90，间隔-0.01
rows001 = len(dataynew001)  # y方向个数，即行数
cols001 = len(dataxnew001)  # x方向个数，即列数
datanew001 = np.zeros((bands, rows001, cols001), dtype=np.float64)  # 建立新data存放插值后的数据############# 创建新文件
datatype = gdal.GDT_Float64
dataset = driver.Create(filename.split('.')[0] + '.tif', cols, rows, bands, datatype)
dataset005 = driver.Create(filename.split('.')[0] + '_005.tif', cols005, rows005, bands, datatype)  # 0.05插值文件
dataset001 = driver.Create(filename.split('.')[0] + '_001.tif', cols001, rows001, bands, datatype)  # 0.03插值文件for i in range(bands):dataset.GetRasterBand(i + 1).WriteArray(data[i])# 0.05插值f1 = interpolate.interp2d(datax, datay, data[i], kind='cubic')datanew005[i] = f1(dataxnew005, dataynew005)dataset005.GetRasterBand(i + 1).WriteArray(datanew005[i])# 0.01插值f2 = interpolate.interp2d(dataxnew005, dataynew005, datanew005[i], kind='cubic')datanew001[i] = f2(dataxnew001, dataynew001)dataset001.GetRasterBand(i + 1).WriteArray(datanew001[i])del dataset001
del dataset005
del dataset
del ds

2.与第一个一样，只是改为了批量读取

# read .grib1 and .grib2 file and output a .tif image
# completed at 2018/11/23
# writen by Zhou Dengji
import os
import sys
from osgeo import gdal
from gdalconst import *
from scipy import interpolate
import numpy as np
import pylab as plfilepath = 'E:/data'# 数据文件夹os.chdir(filepath)# 当前工作目录driver = gdal.GetDriverByName("GTiff")#设置driver为“GTiff”格式
gdal.AllRegister()pathDir = os.listdir(filepath)#读取文件夹内的所有文件，形成listfor filename in pathDir:#对每个文件处理ds = gdal.Open(filename, GA_ReadOnly)#打开if ds is None:print('Cannot open this .grb file.')sys.exit(1)cols = ds.RasterXSize#列数rows = ds.RasterYSize#行数bands = ds.RasterCount#波段数data = np.zeros((bands, rows, cols), dtype=np.float64)#用来存放结果的数组，与源文件的size具有相同的bands/rows/cols，float64类型for i in range(bands):#获取每个波段的datadata[i] = ds.GetRasterBand(i + 1).ReadAsArray(0, 0, cols, rows)datax = np.arange(0, 360, 0.125)  # 原图的x方向尺寸，从0~360，间隔0.125datay = np.arange(90, -90.125, -0.125)  # 原图的y方向尺寸，从90~-90，间隔-0.125dataxx, datayy = np.meshgrid(datax, datay)  # 以datax/datay建立格网############## 0.05间隔的插值dataxnew005 = np.arange(0, 360, 0.05)  # 插值后的x方向尺寸，从0~360，间隔0.05dataynew005 = np.arange(90, -90.125, -0.05)  # 插值后的y方向尺寸，从90~-90，间隔-0.05rows005 = len(dataynew005)  # y方向个数，即行数cols005 = len(dataxnew005)  # x方向个数，即列数datanew005 = np.zeros((bands, rows005, cols005), dtype=np.float64)  # 建立新data存放插值后的数据############## 0.01间隔的插值dataxnew001 = np.arange(0, 360, 0.01)  # 插值后的x方向尺寸，从0~360，间隔0.01dataynew001 = np.arange(90, -90.125, -0.01)  # 插值后的y方向尺寸，从90~-90，间隔-0.01rows001 = len(dataynew001)  # y方向个数，即行数cols001 = len(dataxnew001)  # x方向个数，即列数datanew001 = np.zeros((bands, rows001, cols001), dtype=np.float64)  # 建立新data存放插值后的数据############# 创建新文件datatype = gdal.GDT_Float64dataset = driver.Create(filename.split('.')[0] + '.tif', cols, rows, bands, datatype)dataset005 = driver.Create(filename.split('.')[0] + '_005.tif', cols005, rows005, bands, datatype)  # 0.05插值文件dataset001 = driver.Create(filename.split('.')[0] + '_001.tif', cols001, rows001, bands, datatype)  # 0.03插值文件for i in range(bands):dataset.GetRasterBand(i + 1).WriteArray(data[i])#将没有插值的数据、直接转换成额数据转换成tif图片# 0.05插值f1 = interpolate.interp2d(datax, datay, data[i], kind='cubic')datanew005[i] = f1(dataxnew005, dataynew005)dataset005.GetRasterBand(i + 1).WriteArray(datanew005[i])#保存0.05插值文件成tif# 0.01插值f2 = interpolate.interp2d(dataxnew005, dataynew005, datanew005[i], kind='cubic')datanew001[i] = f2(dataxnew001, dataynew001)dataset001.GetRasterBand(i + 1).WriteArray(datanew001[i])#保存0.01插值文件成tifdel dataset001
del dataset005
del dataset
del ds

3.使用PIL模块将tif文件转成带有透明通道的png图片

由于GDAL中的Create()函数不能直接创建PNG图片，但是CreateCopy()函数可以，因此主要思路就是用GDAL读取GTiff文件，然后用CreateCopy()函数保存成PNG格式，再用PIL库赋予PNG图片一个透明通道。

但是该气象文件的数值范围是[0-1]，格式是float64.因此在GDAL中不能直接从GTiff用CreateCopy函数直接保存成PNG格式，因为保存成PNG格式必须是8bit的Byte类型和16bit的UInt16类型，因此采用一个手段是将数值data*100再转换成UInt16格式，保存成中间文件的tif图片。因此，读取一个0.tif文件，需要生成UInt16类型的0_uint16.tif中间文件，最后生成0.png。但是此时GDAL生成的PNG图片没有透明通道，需要用PIL库读取该png图片，将其有“L”灰度转换至“LA”模式，A就代表alpha透明通道,alpha通道要求范围为[0,255]，把像素值由[0,100]范围扩展到[0,255]范围，用putalpha()函数对每个像素赋予不同的透明度。

1).对单张图片由tif转至带有透明通道的png图片

from osgeo import gdal
from PIL import Image
from gdalconst import *
import numpy as npinputFile = 'E:/0_005.tif'  # 输入的tif文件
midputFile = 'E:/0_005_tif.tif'  # 中间生成的Uint16格式的tif文件
outputFile = 'E:/0_005_png.png'  # 最后生成的png图片
verbose = Truegdal.AllRegister()
driverpng = gdal.GetDriverByName('PNG')  # PNG驱动
drivertif = gdal.GetDriverByName('GTiff')  # GTiff驱动src_ds = gdal.Open(inputFile)  # 读取tif文件cols = src_ds.RasterXSize  # 列数
rows = src_ds.RasterYSize  # 行数
band1 = src_ds.GetRasterBand(1)  # 读取波段，该文件只有一个波段data = band1.ReadAsArray()  # 读取数据
data = (data * 100).astype(np.uint16)  # 将数据*100并设置格式为uint16tifds = drivertif.Create(midputFile, cols, rows, 1, GDT_UInt16)  # 保存中间tif文件,uint16格式
band = tifds.GetRasterBand(1)  # 读取中间文件的波段
band.WriteArray(data, 0, 0)  # 将[0,100]范围的uint16类型data写入中间文件，否则不能保存为png图片
del tifdsmid_ds = gdal.Open(midputFile)  # 读取中间文件，这时的mid_us是uint16类型
dst_ds = driverpng.CreateCopy(outputFile, mid_ds)  # 保存成png图片
del mid_ds
del dst_dsimg = Image.open(outputFile)  # 用PIL模块打开
img = img.convert('LA')  # 转换至含有alpha格式的LA模式
r, a = img.split()  # 分离每个波段
ra = r.point(lambda i: i * 255 / 100)  # 将[0,100]像素值扩展到[0,255]当成透明度值
img.putalpha(ra)  # 将透明度值设到透明通道中
img.save(outputFile)  # 保存del src_ds

2).批量转换

# read the folder all files and convert them to png with transparent channel
# completed at 2018/11/24
# writen by Zhou Dengjifrom osgeo import gdal
from PIL import Image
from gdalconst import *
import numpy as np
import os
import sysgdal.AllRegister()
driverpng = gdal.GetDriverByName('PNG')
drivertif = gdal.GetDriverByName('GTiff')filepath = 'E:/data/ecGRBfiles/TCC_TIF'  # 文件夹
os.chdir(filepath)  # 设置当前工作目录
pathDir = os.listdir(filepath)  # 读取文件夹下的所有文件for filename in pathDir:  # 对每个文件ds = gdal.Open(filename, GA_ReadOnly)if ds is None:print('Cannot open this file:' + filename)sys.exit(1)inputFile = filename  # 输入文件midputFile = inputFile.split('.')[0] + '_tif.tif'  # 生成的中间文件outputFile = inputFile.split('.')[0] + '_png.png'  # 最终生成的png图片src_ds = gdal.Open(inputFile)  # GDAL打开源文件cols = src_ds.RasterXSize  # 列数rows = src_ds.RasterYSize  # 行数band1 = src_ds.GetRasterBand(1)  # 读波段，只有一个波段data = band1.ReadAsArray()  # 读数值data = (data * 100).astype(np.uint16)  # 将[0,1]转换至[0,100]并设置为uint16格式，否则不能保存为png图片tifds = drivertif.Create(midputFile, cols, rows, 1, GDT_UInt16)  # 创建uint16类型的中间文件band = tifds.GetRasterBand(1)  # 获取中间文件的波段band.WriteArray(data, 0, 0)  # 向中间文件写数值del tifdsmid_ds = gdal.Open(midputFile)  # 打开中间文件dst_ds = driverpng.CreateCopy(outputFile, mid_ds)  # 保存成png图片del mid_dsdel dst_dsimg = Image.open(outputFile)  # 用Image读取文件img = img.convert('LA')  # 转换至含有alpha格式的LA模式r, a = img.split()  # 分离每个波段ra = r.point(lambda x: x * 255 / 100)  # 将[0,100]像素值扩展到[0,255]当成透明度值img.putalpha(ra)  # 将透明度值设到透明通道中img.save(outputFile)  # 保存del src_ds

另：用python读取GRIB文件时查看文件数据标签的方法：先用GDAL读取文件，然后查看有多少的通道，获取每个通道后用GetDescription()函数和GetMetadata_Dict()函数获取数据标签，还可以针对性地用GetMetadataItem()查看特定字段的属性。

>>>from osgeo import gdal
>>>from gdalconst import *
>>>filename = 'E:/0.GRB'
>>>ds = gdal.Open(filename, GA_ReadOnly)
>>>ds.RasterXSize
Out[6]: 2880
>>>ds.RasterYSize
Out[7]: 1441
>>>ds.RasterCount
Out[8]: 1
>>>band=ds.GetRasterBand(1)
>>>band.GetDescription()
Out[10]: '0[-] SFC (Ground or water surface)'
>>>band.GetMetadata_Dict()
Out[11]:
{'GRIB_COMMENT': 'Total cloud cover (0 - 1) [-]','GRIB_ELEMENT': 'TCC','GRIB_FORECAST_SECONDS': '0 sec','GRIB_REF_TIME': '  1541592000 sec UTC','GRIB_SHORT_NAME': '0-SFC','GRIB_UNIT': '[-]','GRIB_VALID_TIME': '  1541592000 sec UTC'}
>>>band.GetMetadataItem('GRIB_COMMENT')
Out[12]: 'Total cloud cover (0 - 1) [-]'

关于：

1.透明通道的设置参考：第二篇 Python图片处理模块PIL（pillow）中Putalpha函数

2.Python图像处理库PIL中图像格式转换（一）

3.【python图像处理】给图像添加透明度（alpha通道）