気象データGribフォーマット解析のPythonコードとMatlabコード
以grb/.grb1/.grb 2は拡張子であり、気象データには雲量、雪深、気圧、風速などの複数の内容、または時間系列を有する雲量などを格納することができる.これらのファイルは直接画像に開くことはできません.gribデータを直感的に表示するには、ファイルを読み出して解析し、tifまたはpng形式に保存する必要があります.
ここ数日,matlabコードとpythonコードをそれぞれtif/png形式のピクチャに解析し,任意の必要な解像度に補間した.matlabの実行速度が遅く、データ量が大きすぎて停止しやすいため、pythonを書き直しました.例えば、元のデータの緯度分解能はいずれも0.125°であり、matlabは0.02°にしか補間できず、0.01°に補間できない.pythonはできます.元のデータは3.95 Mで、pythonで0.01解像度に補間したファイルサイズは4.85 Gです.データ量の大きさがわかります.緯度の解像度は、「MatLab部分」の最後のコードから生成することができる.xlsxファイルで表示するには、緯度と対応する行列数で計算することもできます.
【MatLab部】
gribファイルには、データラベルがあります.たとえば、ファイル内部の配列フォーマットがdata(214412880)であれば、data(0)は風速、data(1)は雲量などです.次の展示枠は1つ展示されています.grbファイルのデータラベル、具体的なコードの後ろにあります.
このデータの次元は経度2880 x次元1441であり、1つの期間しかなく、変数にはTotal_cloud_cover_surfaceは、地表雲量情報である.
1.単一gribファイルを読み込み、出力データラベル、データラベルをtxtに格納し、tifピクチャ形式でデータを保存する
2.gribファイルを一括して読み取りtif形式に保存し、そのファイルデータラベルをtxtファイルとして保存する.ファイル名の拡張子に「SOIL」がある場合は、ファイルを処理しません.
3.フォルダ内のすべてのtifファイルを一括して透明なチャネルを持つpngピクチャに変換する
4.ソースファイルの緯度解像度がいずれも0.125のものを0.01に補間し、3次補間を採用し、画像であるためinterp 2()関数を用いて2次元補間を実現する
【Python部】
GDALを使用してGRIBファイルを読み込み、補間関数を使用して0.125解像度を0.05および0.01に補間します.
1.GDALで1枚のGRIBファイルを処理し、scipyモジュールのinterp 2 d関数で2次元補間を行い、0.05と0.01の解像度で補間する
2.最初と同じように、一括読取に変更するだけです
3.PILモジュールを使用してtifファイルを透明なチャネルを持つpngピクチャに変換する
GDALのCreate()関数はPNGピクチャを直接作成することはできないが、CreateCopy()関数は可能であるため、GDALでGTiffファイルを読み出し、CreateCopy()関数でPNGフォーマットに保存し、PILライブラリでPNGピクチャに透明なチャネルを与えることが主な考え方である.
ただし、この気象ファイルの数値範囲は[0-1]、フォーマットはfloat 64である.そのためGDALではGTiffからCreateCopy関数で直接PNG形式に保存することはできません.PNG形式に保存するには8 bitのByteタイプと16 bitのUIT 16タイプでなければならないので、数値data*100をUIT 16形式に再変換し、中間ファイルのtifピクチャに保存するという手段が採用されています.したがって、1つの0を読み出す.tifファイル、UInt 16タイプの0_を生成する必要がありますuint16.tif中間ファイル、最後に0を生成する.png.しかし、このときGDALで生成されたPNGピクチャには透明チャネルがなく、PILライブラリでpngピクチャを読み取り、「L」階調を「LA」モードに変換する必要があり、Aはalpha透明チャネルを表し、alphaチャネル要求範囲は[0525]、画素値を[001]範囲から[0525]範囲に拡張し、putalpha()関数で各画素に異なる透明度を与える.
1).単一のピクチャに対してtifから透明なチャネルを持つpngピクチャに移動する
2).バッチ変換
また、pythonでGRIBファイルを読み込むときにファイルデータラベルを表示する方法:まずGDALでファイルを読み込み、その後どのくらいのチャネルがあるかを確認し、各チャネルを取得した後にGetDescription()関数とGetMetadata_Dict()関数は、データラベルを取得し、GetMetadataItem()を使用して特定のフィールドのプロパティを表示することもできます.
説明:
1.透明チャンネルの設定参考:第2編Python画像処理モジュールPIL(pillow)におけるPutalpha関数
2.Python画像処理ライブラリPILにおける画像フォーマット変換(一)
3.【python画像処理】画像に透明度(alphaチャネル)を追加
ここ数日,matlabコードとpythonコードをそれぞれtif/png形式のピクチャに解析し,任意の必要な解像度に補間した.matlabの実行速度が遅く、データ量が大きすぎて停止しやすいため、pythonを書き直しました.例えば、元のデータの緯度分解能はいずれも0.125°であり、matlabは0.02°にしか補間できず、0.01°に補間できない.pythonはできます.元のデータは3.95 Mで、pythonで0.01解像度に補間したファイルサイズは4.85 Gです.データ量の大きさがわかります.緯度の解像度は、「MatLab部分」の最後のコードから生成することができる.xlsxファイルで表示するには、緯度と対応する行列数で計算することもできます.
【MatLab部】
gribファイルには、データラベルがあります.たとえば、ファイル内部の配列フォーマットがdata(214412880)であれば、data(0)は風速、data(1)は雲量などです.次の展示枠は1つ展示されています.grbファイルのデータラベル、具体的なコードの後ろにあります.
このデータの次元は経度2880 x次元1441であり、1つの期間しかなく、変数には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; // int
float 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 code
function readfile
setup_nctoolbox
[filename, pathname] = uigetfile('*.grb', 'choose a GRB file'); % .grb
if isequal(filename,0)
msgbox('you choose nothing');
else
pathfile=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
end2.gribファイルを一括して読み取りtif形式に保存し、そのファイルデータラベルをtxtファイルとして保存する.ファイル名の拡張子に「SOIL」がある場合は、ファイルを処理しません.
% Matlab code
function readfile
setup_nctoolbox
file_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.GRB
name=strsplit(img_name,'.');% “.” , 20180101_ABC_SOIL GRB
a=strsplit(char(name(1)),'_');% “_” ， a=SOIL
if strcmp(char(a(length(a))),'SOIL')% a==SOIL
continue;
else
pathfile=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));% txt
txtname=[str1,str2,str4];
fid=fopen(txtname,'w');
fprintf(fid,'%s',label);
fclose(fid);
end
end
end3.フォルダ内のすべてのtifファイルを一括して透明なチャネルを持つpngピクチャに変換する
% Matlab code
function converse
%setup_nctoolbox
filepath='E:\ecGRBfiles\TCC\';
img_list=dir(strcat(filepath,'*.tif'));
num=length(img_list);
type='.png';
for k=1:num
img_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
end4.ソースファイルの緯度解像度がいずれも0.125のものを0.01に補間し、3次補間を採用し、画像であるためinterp 2()関数を用いて2次元補間を実現する
% Matlab code
function readfile
setup_nctoolbox
[filename, pathname] = uigetfile('*.grb', 'choose a GRB file'); % .grb
if isequal(filename,0)
msgbox('you choose nothing');
else
pathfile=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*384
alpha=mapminmax(GPMData,0,1);
alpha=double(alpha);
imwrite(mat2gray(GPMData),'cloud_origin.png','Alpha',alpha);
% ， 0.125， 0~359.875， 90~-90
x_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 excel
name = {'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で1枚のGRIBファイルを処理し、scipyモジュールのinterp 2 d関数で2次元補間を行い、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 pl
os.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 ds2.最初と同じように、一括読取に変更するだけです
# 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 pl
filepath = 'E:/data'#
os.chdir(filepath)#
driver = gdal.GetDriverByName("GTiff")# driver “GTiff”
gdal.AllRegister()
pathDir = os.listdir(filepath)# ， list
for 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):# data
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])# 、 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 tif
del dataset001
del dataset005
del dataset
del ds3.PILモジュールを使用してtifファイルを透明なチャネルを持つpngピクチャに変換する
GDALのCreate()関数はPNGピクチャを直接作成することはできないが、CreateCopy()関数は可能であるため、GDALでGTiffファイルを読み出し、CreateCopy()関数でPNGフォーマットに保存し、PILライブラリでPNGピクチャに透明なチャネルを与えることが主な考え方である.
ただし、この気象ファイルの数値範囲は[0-1]、フォーマットはfloat 64である.そのためGDALではGTiffからCreateCopy関数で直接PNG形式に保存することはできません.PNG形式に保存するには8 bitのByteタイプと16 bitのUIT 16タイプでなければならないので、数値data*100をUIT 16形式に再変換し、中間ファイルのtifピクチャに保存するという手段が採用されています.したがって、1つの0を読み出す.tifファイル、UInt 16タイプの0_を生成する必要がありますuint16.tif中間ファイル、最後に0を生成する.png.しかし、このときGDALで生成されたPNGピクチャには透明チャネルがなく、PILライブラリでpngピクチャを読み取り、「L」階調を「LA」モードに変換する必要があり、Aはalpha透明チャネルを表し、alphaチャネル要求範囲は[0525]、画素値を[001]範囲から[0525]範囲に拡張し、putalpha()関数で各画素に異なる透明度を与える.
1).単一のピクチャに対してtifから透明なチャネルを持つpngピクチャに移動する
from osgeo import gdal
from PIL import Image
from gdalconst import *
import numpy as np
inputFile = 'E:/0_005.tif' # tif
midputFile = 'E:/0_005_tif.tif' # Uint16 tif
outputFile = 'E:/0_005_png.png' # png
verbose = True
gdal.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 uint16
tifds = 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 tifds
mid_ds = gdal.Open(midputFile) # ， mid_us uint16
dst_ds = driverpng.CreateCopy(outputFile, mid_ds) # png
del mid_ds
del dst_ds
img = 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_ds2).バッチ変換
# read the folder all files and convert them to png with transparent channel
# completed at 2018/11/24
# writen by Zhou Dengji
from osgeo import gdal
from PIL import Image
from gdalconst import *
import numpy as np
import os
import sys
gdal.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 tifds
mid_ds = gdal.Open(midputFile) #
dst_ds = driverpng.CreateCopy(outputFile, mid_ds) # png
del mid_ds
del dst_ds
img = 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.透明チャンネルの設定参考:第2編Python画像処理モジュールPIL(pillow)におけるPutalpha関数
2.Python画像処理ライブラリPILにおける画像フォーマット変換(一)
3.【python画像処理】画像に透明度(alphaチャネル)を追加