2016年09月08日

土壌雨量指数grib2形式の読み込みプログラム

気象庁が出している土壌雨量指数を
読み込むプログラムをC#の勉強を
兼ねて作ってみた。

全然、わからなかったので、
「ロケッこが行く」のブログ
http://blog.syo-ko.com/?eid=1228
を大変参考にさせてもらった。
ありがとうございます。

grib2形式のバイナリーファイルでできている。
データ部はランレングスという形で表現されていて
この読み方に結構時間を取られた。

ちなみに、気象庁が出している技術資料は
http://www.data.jma.go.jp/add/suishin/jyouhou/pdf/302.pdf
から取れる。
でも、タンクの数が違ったり、まだまだ謎が多い。
自分の理解不足か、よくわからない。

ちょっとづつC#の勉強しながらなので、
全然クラスの使い方とか、ファイルの読み込み方とか
統一とれてない上に、動作確認のためのコンソールへの
出力とかをコメントアウトのまま残してる。

ネットで検索すると、土壌雨量指数の読み込みプログラムが、
異常な値段で売られてたりするので、
少しでも役に立てればうれしい。

ただ、実行は自己責任でお願いします。

できたテキストファイルをRで描画すると
ある日の土壌雨量指数
それっぽいので、たぶんうまくいけてる。

で、ソースは以下のとおり。

using System;
using System.IO;

class DojoUryouANAL
{


static void Main()
{

int FileYear = 2011; //入力年の指定
int[] nMonth = new int[] { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
nMonth[2 - 1] = 28 + (1 / (FileYear % 4 + 1)) * (1 - 1 / (FileYear % 100 + 1)) + (1 / (FileYear % 400 + 1)); //2月のうるう年判別

int nCount = 0;
int FileMonth = 1;
int FileDay = 1;
int FileHour = 0;
int FileMinute = 0;
int FileSecond = 0;

string InputFile1 = "InputFiles¥¥Z__C_RJTD_";
string InputFIle2 = "_SRF_GPV_Gll5km_P-swi_ANAL_grib2.bin";

string[] InputFileNameList = new string[18000];

while (FileMonth < 13)
{
while (FileDay < nMonth[FileMonth - 1] + 1)
{
while (FileHour < 24)
{
for (int minIndex = 0; minIndex < 2; minIndex++)
{
if (minIndex == 0)
{
FileMinute = 0;
InputFileNameList[nCount] = InputFile1 + FileYear.ToString("D4") + FileMonth.ToString("D2") + FileDay.ToString("D2") + FileHour.ToString("D2") + FileMinute.ToString("D2") + FileSecond.ToString("D2") + InputFIle2;
//Console.WriteLine("InputFile = " + InputFileName);
// Console.WriteLine("Year={0},Month={1},Day={2},Hour={3},Minute={4},Second={5}", Year, Month, Day, Hour, Minute, Second);
nCount++;
}
else
{
FileMinute = 30;
InputFileNameList[nCount] = InputFile1 + FileYear.ToString("D4") + FileMonth.ToString("D2") + FileDay.ToString("D2") + FileHour.ToString("D2") + FileMinute.ToString("D2") + FileSecond.ToString("D2") + InputFIle2;
//Console.WriteLine("InputFile = " + InputFileName);
//Console.WriteLine("Year={0},Month={1},Day={2},Hour={3},Minute={4},Second={5}", Year, Month, Day, Hour, Minute, Second);
nCount++;
}
}
FileHour++;
}
FileHour = 0;
FileDay++;
}

FileDay = 1;

FileMonth++;

}


for (int sCount = 0; sCount < nCount; sCount++)
{

string InputFileName = InputFileNameList[sCount];
//string InputFileName = "Z__C_RJTD_20110903170000_SRF_GPV_Gll5km_P-swi_ANAL_grib2.bin";
//string OutputFileName = "output.txt";
string FileHeadS = InputFileName.Substring(11, 10);
string FileYearS = InputFileName.Substring(21, 4);
string FileMonthS = InputFileName.Substring(25, 2);
string FileDayS = InputFileName.Substring(27, 2);
string FileHourS = InputFileName.Substring(29, 2);
string FileMinuteS = InputFileName.Substring(31, 2);
string FileSecondS = InputFileName.Substring(33, 2);


//string OutputFileName = "Dojosisu_ANAL_" + FileYearS + FileMonthS + FileDayS + FileHourS + FileMinuteS + ".txt";
//Console.WriteLine("OutputFileName = " + OutputFileName);

FileStream dataIN; //入力データファイル
StreamWriter fstr_out; //出力データファイル

//ファイルの読み込み(はじめ)
try
{
dataIN = new FileStream(InputFileName, FileMode.Open, FileAccess.Read);
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
return;
}

int fileSize = (int)dataIN.Length; //オクテット数
byte[] rByte = new byte[fileSize]; //rByte 全オクテットの保持配列

Console.WriteLine("File Size: " + fileSize);

int readSize;
int remain = fileSize;
int bufPos = 0; //読み込み開始ポイント

while (remain > 0)
{
readSize = dataIN.Read(rByte, bufPos, Math.Min(1, remain));

bufPos += readSize;
remain -= readSize;
}
//ファイルの読み込み(終わり)


//ファイルを閉じる。
dataIN.Dispose();




//////解釈の始まり//////
//1〜4オクテット
// for (int i = 0;i < 4; i++){
// Console.Write((char)rByte[i]+",");
// }

//9〜16オクテット
// long rValue = ByteToInt(rByte, 8, 8);
long rValue = ByteToSingedInt(rByte, 8, 8);
Console.Write(rValue + ",");



Console.WriteLine();
Console.WriteLine("第1節");
//第1節1〜4オクテット ※節の長さ
rValue = ByteToSingedInt(rByte, 16 + 1 - 1, 4);
Console.WriteLine(rValue);
int lenFirstchp = (int)rValue;
//第1節5オクテット
rValue = ByteToSingedInt(rByte, 16 + 5 - 1, 1);
Console.WriteLine(rValue);
//第1節6〜7オクテット
rValue = ByteToSingedInt(rByte, 16 + 6 - 1, 2);
Console.WriteLine(rValue);
//第1節13〜14オクテット
rValue = ByteToSingedInt(rByte, 16 + 13 - 1, 2);
Console.WriteLine(rValue);
//第1節15オクテット
rValue = ByteToSingedInt(rByte, 16 + 15 - 1, 1);
Console.WriteLine(rValue);
//第1節16オクテット
rValue = ByteToSingedInt(rByte, 16 + 16 - 1, 1);
Console.WriteLine(rValue);
//第1節17オクテット
rValue = ByteToSingedInt(rByte, 16 + 17 - 1, 1);
Console.WriteLine(rValue);
//第1節18オクテット
rValue = ByteToSingedInt(rByte, 16 + 18 - 1, 1);
Console.WriteLine(rValue);
//第1節19オクテット
rValue = ByteToSingedInt(rByte, 16 + 19 - 1, 1);
Console.WriteLine(rValue);
//第1節20オクテット
rValue = ByteToSingedInt(rByte, 16 + 20 - 1, 1);
Console.WriteLine(rValue);


Console.WriteLine();
Console.WriteLine("第3節");

//第3節1〜4オクテット ※節の長さ
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 1 - 1, 4);
Console.WriteLine(rValue);
int lenThirdchp = (int)rValue;
//第3節5オクテット ※節番号
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 5 - 1, 4);
Console.WriteLine(rValue);
//第3節31〜34オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 31 - 1, 4);
Console.WriteLine(rValue);
//第3節35〜38オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 35 - 1, 4);
Console.WriteLine(rValue);
//第3節47〜50オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 47 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第3節51〜54オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 51 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第3節56〜59オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 56 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第3節60〜63オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 60 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第3節64〜67オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 64 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第3節68〜71オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 68 - 1, 4);
Console.WriteLine(rValue / 1000000.0);



int lenNow = 0;

for (int NumTank = 0; NumTank < 4; NumTank++)
{



Console.WriteLine();
Console.WriteLine("第4節");



//第4節1〜4オクテット ※節の長さ
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 1 - 1, 4);
Console.WriteLine(rValue);
int lenForthchp = (int)rValue;

//第4節5オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 5 - 1, 1);
Console.WriteLine(rValue);

Console.WriteLine();
Console.WriteLine("第5節");

//第5節1〜4オクテット ※節の長さ
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 1 - 1, 4);
Console.WriteLine(rValue);
int lenFivechp = (int)rValue;

//第5節5オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 5 - 1, 1);
Console.WriteLine(rValue);
//第5節6−9オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 6 - 1, 4);
Console.WriteLine(rValue);

//第5節10−11オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 10 - 1, 2);
Console.WriteLine(rValue);

//第5節12オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 12 - 1, 1);
Console.WriteLine("nbit: " + rValue);
int nbit = (int)rValue;

//第5節13−14オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 13 - 1, 2);
Console.WriteLine("maxv: " + rValue);
int maxv = (int)rValue;

//第5節15−16オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 15 - 1, 2);
Console.WriteLine("max_M: " + rValue);
int max_M = (int)rValue;

//第5節17オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 17 - 1, 1);
Console.WriteLine(rValue);

int[] LevelValue = new int[max_M + 1];
for (int level = 1; level < max_M + 1; level++)
{
//第5節18−オクテット レベルの対応値
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 18 + (level - 1) * 2 - 1, 2);
LevelValue[level] = (int)rValue / 10;
Console.WriteLine("level: " + level + " Value: " + rValue);
Console.WriteLine("level: " + level + " LevelValue " + LevelValue[level]);
}
LevelValue[0] = 9999; ///欠測地

Console.WriteLine();
Console.WriteLine("第6節");


//第6節1−4オクテット ※節の長さ
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenNow + 1 - 1, 4);
Console.WriteLine(rValue);
int lenSixchp = (int)rValue;

//第6節5オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenNow + 5 - 1, 1);
Console.WriteLine(rValue);
//第6節6オクテット
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenNow + 6 - 1, 1);
Console.WriteLine(rValue);


Console.WriteLine();
Console.WriteLine("第7節");


//第7節1−4オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenSixchp + lenNow + 1 - 1, 4);
Console.WriteLine("length of this chapter: " + rValue);
int n_data = (int)rValue - 6 + 1;
int lenSevenchp = (int)rValue;

//第7節5オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenSixchp + lenNow + 5 - 1, 1);
Console.WriteLine("number of chapter: " + rValue);



//出力ファイル名の操作
string OutputFileName = "OutputFIles¥¥Dojosisu_ANAL_" + "TANK" + NumTank.ToString() + "_" + FileYearS + FileMonthS + FileDayS + FileHourS + FileMinuteS + "UTC.txt";
Console.WriteLine("OutputFileName = " + OutputFileName);


int j;
int p = -1, m = 1, n = 0, k = 0;
int v = 0;
int l = ipow(2, nbit) - 1 - maxv;

// n_data = 30;

try
{
fstr_out = new StreamWriter(OutputFileName); //土壌雨量指数の出力ファイル
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
return;
}

for (int iRun = 1; iRun < n_data + 1; iRun++)
{
// Console.WriteLine("iRun: " + iRun);

//第7節6−オクテット
// rValue = ByteToSingedInt(rByte, 16 + 21 + 72 + 58 + lenFivechp + 6 + iRun + 5 - 1, 1);
// Console.WriteLine("iRun: " + iRun + " Value: " + rValue);
// Console.WriteLine("iRun: " + iRun + " Value: " + rByte[16 + 21 + 72 + 58 + lenFivechp + 6 + iRun + 5 - 1]);
// Console.WriteLine("iRun: " + iRun + " Value: " + rByte[16 + 21 + 72 + 82 + lenFivechp + 6 + iRun + 5 - 1].ToString("X"));

v = (int)ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenSixchp + lenNow + iRun + 5 - 1, 1);
// Console.Write(v + ",");
if (v <= maxv)
{
if (p >= 0)
{
// Console.WriteLine(" ,m: " + m);

for (j = 0; j < m; j++)
{
k++;
// Console.Write(p+" , ");
try
{
fstr_out.Write(LevelValue[p] + " , ");
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
break;
}
}
// Console.WriteLine("p: " + p + " ,m: " + m );
p = v;
m = 1;
n = 0;
}
p = v;
m = 1;
n = 0;

}
else
{
m += ipow(l, n) * (v - maxv - 1);
// Console.WriteLine(ipow(l, n) * (v - maxv - 1)+","+m);
n++;
}

}

for (j = 0; j < m; j++)
{
k++;
if (k != 286720)
{
// Console.Write(p);
try
{
fstr_out.Write(LevelValue[p] + " , ");
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
break;
}
}
else
{
// Console.Write(p+" , ");
try
{
fstr_out.Write(LevelValue[p] + " ");
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
break;
}
}
}
Console.WriteLine("lenNow:" + lenNow);
lenNow += lenForthchp + lenFivechp + lenSixchp + lenSevenchp;
Console.WriteLine("lenNow:" + lenNow);
fstr_out.Close();
Console.WriteLine("k: " + k); //286,720=512x560

}

//第8節1−4オクテット ※-1
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 0 - 1, 1);
Console.WriteLine("¥n -1");
Console.WriteLine("Message: " + rValue);


//第8節1−4オクテット ※1
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 1 - 1, 1);
Console.WriteLine("¥nEighth Chapter 1");
Console.WriteLine("Message: " + rValue);

//第8節1−4オクテット ※2
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 2 - 1, 1);
Console.WriteLine("¥nEighth Chapter 2");
Console.WriteLine("Message: " + rValue);

//第8節1−4オクテット ※3
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 3 - 1, 1);
Console.WriteLine("¥nEighth Chapter 3");
Console.WriteLine("Message: " + rValue);

//第8節1−4オクテット ※4
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 4 - 1, 1);
Console.WriteLine("¥nEighth Chapter 4");
Console.WriteLine("Message: " + rValue);

}
}
//Mainメソッド 終了


//整数の累乗を計算するメソッド
static int ipow(int i, int j)
{
int k, l;
for (k = 0, l = 1; k < j; k++) l *= i;
return (l);
}


//バイトを、整数に変換するメソッド
static long ByteToInt(byte[] source, int Start, int Length)
{
long rVal = 0;
int intB;
int intE = Start + Length;

for (intB = Start; intB < intE; intB++)
{
rVal *= 256;
rVal += Convert.ToInt64(source[intB]);

//Console.WriteLine("intB: " + intB + " ,souce: " + Convert.ToInt64(source[intB]) + " ,rval: " + rVal);
}

return rVal;

}


//バイトを、符号付き整数に変換するメソッド
static long ByteToSingedInt(byte[] source, int Start, int Length)
{
long rVal = 0;
int intB;
int intE = Start + Length;
bool bSng = true;

//先頭オクテットの確認
if ((source[Start] & 0x80) == 0x80)
{
// Console.WriteLine("負数である");
bSng = true;
source[Start] = (byte)(source[Start] & 0x7F);
}
else
{
bSng = false;
// Console.WriteLine("正数である");
}


for (intB = Start; intB < intE; intB++)
{
rVal *= 256;
rVal += Convert.ToInt64(source[intB]);

//Console.WriteLine("intB: " + intB + " ,souce: " + Convert.ToInt64(source[intB]) + " ,rval: " + rVal);

}

if (bSng == true)
{
rVal *= -1;
}

return rVal;

}


}




同じカテゴリー(スタッフ日記)の記事画像
GMT5.3.1で陰影図や鳥瞰図を作る
GMT5.3.1で矢印とか地図をかく
QGISで矢印を表現させる
Mac mini Late2012のHDDをSSDに換装
zenfone3 の購入
河川シミュレーションソフトiRIC
同じカテゴリー(スタッフ日記)の記事
 2017初頭のウルトラブック比較 (2017-02-26 19:08)
 Latexをwindowsにインストール (2017-02-18 21:28)
 Javaを使った有限体積法による数値計算(1) (2017-02-16 21:07)
 GMT5.3.1で陰影図や鳥瞰図を作る (2016-12-12 15:55)
 GMT5.3.1で矢印とか地図をかく (2016-12-09 00:08)
 QGISで矢印を表現させる (2016-12-07 20:38)

Posted by 和歌山サイエンスカフェインフィニティ at 20:31│Comments(0)スタッフ日記
上の画像に書かれている文字を入力して下さい
 
<ご注意>
書き込まれた内容は公開され、ブログの持ち主だけが削除できます。