ModeSは地上のATCレーダーからのコマンドに応答する形でメッセージを送信する 航空交通管制のほかローカルに設置した受信機によって収集できます。 図はこちら

		Vtas:true air speed Vgs:ground speed Vw:風速 χa:magnetic heading χg:track angle χw:風向

# (C) 2024/06/19 Dr.Hiroshi Ishikawa
import socket

def main():
	host = "192.168.1.xx"#RaspberryPi #1のアドレス
	port = 30002
	n = 0
	try:
		with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as s:
			s.connect((host, port))
			print(f"Connected to {host} on port {port}")
			while n <= 10:
				data = s.recv(1024)
				if not data:
					break
				while b'\n' in data:
					line, data = data.split(b'\n', 1)
					print(line.decode())
					n = n + 1
	except Exception as e:
		print(f"An error occurred: {e}")

if __name__ == "__main__":
	main()

+----------+----------+-------------+------------------------+-----------+
|  DF (5)  |  CA (3)  |  ICAO (24)  |         ME (56)        |  PI (24)  |
+----------+----------+-------------+------------------------+-----------+

１．ダウンロード
https://github.com/junzis/pyModeS　の　<>Code ボタンをクリックしDownloadZIP でダウンロード
展開しその中の　pyModeS ディレクトリを　RaspberryPi #3のホームディレクトリにコピー

２．テスト
RaspberryPi #3のコンソールから次のような試験を行います。

pythonを起動し、pyModeSをインポートし、いくつかfunctionを使ってみます。
上で例示したメッセージのデコードしてみます。
df機能を使ってDFが17、bds.infer機能を使ってBDSが09、すなわちADSBメッセージのAirborne velocityであることがわかります。

pi@modeS:~ $ python Python 3.11.2 (main, Mar 13 2023, 12:18:29) [GCC 12.2.0] on linux Type "help", "copyright", "credits" or "license" for more information. >>> import pyModeS as pms >>> pms.df("8D84B81C990DAD8A90083884D4F9") 17 >>> pms.bds.infer("8D84B81C990DAD8A90083884D4F9") 'BDS09

さらにtell機能により、全体の意味は次のようでGround speedが示されます。

>>> pms.tell("8D84B81C990DAD8A90083884D4F9") 　　　　Message: 8D84B81C990DAD8A90083884D4F9 　　　　ICAO address: 84B81C 　　　　Downlink Format: 17 　　　　Protocol: Mode-S Extended Squitter (ADS-B) 　　　　Type: Airborne velocity 　　　　Speed: 435 knots 　　　　Track: 259.0251267445889 degrees 　　　　Vertical rate: 64 feet/minute 　　　　Type: Ground speed

ModeSのBDS50の例を示します。Track angle、　Ground speed　True airspeed、がみえます。

>>> pms.tell("A80006ACF9363D3BBF9CE98F1E1D") 　　　　Message: A80006ACF9363D3BBF9CE98F1E1D 　　　　ICAO address: 4008B4 　　　　Downlink Format: 21 　　　　Protocol: Mode-S Comm-B identity reply 　　　　Squawk code: 6322 　　　　BDS: BDS50 (Track and turn report) 　　　　Roll angle: -9.66796875 degrees 　　　　Track angle: 140.2734375 degrees 　　　　Track rate: -0.40625 degree/second 　　　　Ground speed: 476 knots 　　　　True airspeed: 466 knots

ModeSのBDS６０の例を示します。Megnatic Heading、　Indicated airspeed、　Mach numbe　が見えます。

>>> pms.tell("A80004AAA74A072BFDEFC1D5CB4F") 　　　　Message: A80004AAA74A072BFDEFC1D5CB4F 　　　　ICAO address: 4CA53F 　　　　Downlink Format: 21 　　　　Protocol: Mode-S Comm-B identity reply 　　　　Squawk code: 4720 　　　　BDS: BDS60 (Heading and speed report) 　　　　Megnatic Heading: 110.390625 degrees 　　　　Indicated airspeed: 259 knots 　　　　Mach number: 0.7000000000000001 　　　　Vertical rate (Baro): -2144 feet/minute 　　　　Vertical rate (INS): -2016 feet/minute

このほかのpyModeSのすべての機能についてはこちらにあります。

本番用プログラム

RaspberryPi #3で動かす本番プログラムを作成します。Beast AVR 形式のデータから次のデータを取り出す必要があります。

adsbメッセージから ICAO,altitude,latitude,longitude
bds 50から ICAO,trackangle,groundspeed,trueair
bds　60から ICAO,heading,indicatedair,mach

adsbは1秒に2回、modeSデータは数秒に1回送信されます。1つの機体から別々のメッセージでデータが出力されるので、これらをICAOアドレスごとに一つの[icao,adsbデータ,bds50データ,bds60データ]という形に整理しないといけません。この「名寄せ」のようなもので、そのためのサブプログラムがaircraft(icao , code, data)関数です。

本番用プログラム 　modeSdata.py

# (C) 2024/06/19 Dr.Hiroshi Ishikawa import socket import pyModeS as pms import os from datetime import datetime import time #処理中の航空機のリスト　[icao,adsbデータ,bds50データ,bds60データ] aircrafts = [ [0,0,0,0], ] #受信データを1行分切り出す関数 def receive_lines(sock): buffer = b'' # 受信データを格納するバッファ while True: data = sock.recv(1024) # 1024バイトごとにデータを受信 if not data: # データが空の場合は接続が切断されたと判断してループを終了 print(timestamp() + " no data") break buffer += data # 受信したデータをバッファに追加 while b'\n' in buffer: # 改行文字が含まれるかチェック line, buffer = buffer.split(b'\n', 1) # バッファから1行分のデータを切り出す yield line.decode() # 1行分のデータをyieldする if buffer: # バッファに残っているデータがある場合はそれを処理する yield buffer.decode() #航空機リストをICAOアドレスごとに整理する関数 def aircraft(icao , code, data): global aircrafts #航空機リスト try: index = search_icao(aircrafts, icao) #aircraftsリストを検索しICAOが見つかるとリストの何番目かが返える if index is not None:#リストにICAOが見つかった場合 aircrafts[index][0] = icao #0番目にicaoを代入 if code == 1: aircrafts[index][1] = data #adsbの場合は1番目にadsbデータを代入 if code == 2: aircrafts[index][2] = data #bds50の場合は2番目にbds50データを代入 if code == 3: aircrafts[index][3] = data #bds60の場合は3番目にbds50データを代入 #航空機のadsb,bds50,bds60のデータがすべて整ったらタイムスタンプをつけてファイルに書き出す if (aircrafts[index][1] != 0) and (aircrafts[index][2] != 0) and (aircrafts[index][3]!= 0) : result = aircrafts[index][0] + "," + str(aircrafts[index][1])+ "," + str(aircrafts[index][2]) + "," + str(aircrafts[index][3]) filename = "modeSdata/" + datetime.now().strftime("%Y%m%d%H")+ "modeS_data.txt"#1時間に1ファイル with open(filename, "a" ) as f: f.write(timestamp() + "," + result + "\n") f.close() del aircrafts[index] #航空機リストにゴミがたまり50以上になったら古い10のリストを削除する if len(aircrafts) > 50: del aircrafts[0:9] else:#見つからないときは新しい航空機であるからリストに次の値を追記する if code == 1: aircrafts.append([icao,data,0,0]) if code == 2: aircrafts.append([icao,0,data,0]) if code == 3: aircrafts.append([icao,0,0,data]) """ 航空機リストの挙動を画面で確認するときはコメントアウトを解除 os.system('clear') for i in aircrafts: print(i) """ except Exception as e: print(timestamp() + f" aircraft error : {e}") #search_icao(aircrafts, icao)とすればaircraftsリストを検索しicaoが見つかるとリストの何番目かを返す def search_icao(list, icao): for index, info in enumerate(list): if info[0] == icao: return index return None def timestamp(): return datetime.now().strftime("%Y/%m/%d %H:%M:%S") def main(): # ホストとポート番号 host = "192.168.1.xx" #RaspberryPi #1のアドレス port = 30002 #ポート番号は　30002 global aircrafts #航空機リスト try: # ソケットを作成し、ポート30002に接続 with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as s: s.settimeout(20) #航空機が飛ばない時間帯があるのでタイムアウトを設定 s.connect((host, port)) print(timestamp() + f" Connected to {host} on port {port}") lat_ref = 35.0 #大まかな位置情報が必要 lon_ref = 139.0 icao = "" alt = 0.0 lat = 0.0 lon = 0.0 mach =0.0 tas = 0.0 for line in receive_lines(s):#1行すなわち1メッセージを受信 msg = line[1 : -1]#先頭の*と最後の;を取り除く if pms.df(msg) == 17:#DF=17すなわちADSBの場合 icao = str(pms.icao(msg))#icaoアドレス if pms.adsb.altitude(msg) != 0: alt = pms.adsb.altitude(msg)#高度 lat = pms.adsb.position_with_ref(msg, lat_ref, lon_ref)[0]#緯度 lon = pms.adsb.position_with_ref(msg, lat_ref, lon_ref)[1]#経度 #aircraft関数をよび航空機リストに書き込む。(alt, lat, lon)はタプル aircraft(icao ,1, (alt, lat, lon))# elif pms.bds.bds50.is50(msg):#BDS=50の場合 icao = str(pms.icao(msg))#icaoアドレス trackangle = pms.commb.trk50(msg) groundspeed = pms.commb.gs50(msg) trueair = pms.commb.tas50(msg) aircraft(icao ,2, (trackangle, groundspeed, trueair)) elif pms.bds.bds60.is60(msg):#BDS=50の場合 icao = str(pms.icao(msg))#ICAOアドレス heading = pms.commb.hdg60(msg) indicatedair = pms.commb.ias60(msg) mach = pms.commb.mach60(msg) aircraft(icao ,3, (heading, indicatedair, mach )) except Exception as e: print(timestamp() + f" error: {e}") finally: print(timestamp() + " socket closed") s.close() if __name__ == "__main__": while True: print(timestamp() + " start") main() time.sleep(30)

ホームディレクトリに　modeSdata というフォルダーを作っておきます。
次のように起動します。

pi@modeS:~ $ python modeSdata.py 2024/06/19 10:32:19 start 2024/06/19 10:32:19 Connected to 192.168.1.xx on port 30002

1時間に1本
modeSdata/2024061909modeS_data.txt
というファイルができます。
出力ファイルは1航空機に対し次のような形式です。()はpythonのタプル
　　　yyyy/MM/dd HH:mm:ss,ICAO,(altitude,latitude,longitude),(trackangle,groundspeed,trueair),(heading,indicatedair,mach)

気温、風向・風速を算出するプログラム（主要部分）

　気温、風向・風速を算出するプログラムは、Windowsの親機で実行することとしました。Javaで作りました。上で述べた、重要な気温算出、風の計算はプログラムのおしまいにある5つのサブプログラム（メソッド）です。
trueAirSpeedと double machから気温を算出
　　　double temperature(double trueAirSpeed, double mach)
trackangle とgroundspeedと headingと trueairspeedから風のx成分を算出
　　　double wind_x(double trackangle ,double groundspeed, double heading, double trueairspeed)
trackangle とgroundspeedと headingと trueairspeedから風のy成分を算出
　　　double wind_y(double trackangle ,double groundspeed, double heading, double trueairspeed)
風のx成分と風のy成分から風向を求める
　　　double windangle(double speedx ,double speedy)
風のx成分と風のy成分から風速を求める
　　　 double windspeed(double speedx ,double speedy)

数分に１回動かします。標高100mごとに集計し平均し、図に描くようにしていますが、この部分は省略します。

気温、風向・風速を算出するプログラム 　modeS_weatherdata.java

/*　(C) 2024/06/19 Dr.Hiroshi Ishikawa * 1時間分のyyyyMMddHHmodeSdata.txtを集計 *出力は yyyyMMddHHmodeS_weather.csv　 */ package adsb; import java.io.*; import java.text.SimpleDateFormat; import java.util.Calendar; class modeS_weatherdata{ public static void main(String args[]) throws IOException { SimpleDateFormat DF = new SimpleDateFormat("yyyyMMddHH"); Calendar cal = Calendar.getInstance(); String yyyyMMddHH = DF.format(cal.getTime());//今 String filename0 = "D:\\MyDocuments\\adsb\\modeSweather\\" + yyyyMMddHH + "modeS_weather.csv";//出力ファイル String filename1 = "\\\\192.168.1.zz\\pi\\modeSdata\\" + yyyyMMddHH + "modeS_data.txt";//入力ファイル zzは RaspberryPi#3のアドレス String ss1 ="" ; String[] strAry = new String[40]; String Time = ""; String icao = ""; double altitude = 0.0; double latitude = 0.0; double longitude = 0.0; double trackangle = 135.0; double groundspeed = 14.0; double heading = 180; double trueairspeed = 10;//m/s double mach = 0.0; double speedx = 0.0; double speedy = 0.0; double temperature = 0.0; try { File f0 = new File(filename0); PrintWriter output0 = new PrintWriter(new BufferedWriter(new FileWriter(f0))); //2024/05/31 00:00:02,AC2D32,(33000, 35.56581115722656, 137.96710968017578),(57.3046875, 524, 494),(59.765625, 293, 0.8240000000000001) //yyyy/MM/dd HH:mm:ss,ICAO,(altitude,latitude,longitude),(trackangle,groundspeed,trueair),(heading,indicatedair,mach) File f1 = new File(filename1); BufferedReader input1 = new BufferedReader(new FileReader(f1)); while ((ss1 = input1.readLine()) != null){ ss1 = ss1.replace("(", "");//タプルのカッコを取り除く ss1 = ss1.replace(")", "");//タプルのカッコを取り除く ss1 = ss1.replace("None", "0"); strAry = ss1.split(",");//カンマ区切りをstrAryにいれる Time = strAry[0]; icao = strAry[1]; if (!(strAry[2].equals("None"))){ altitude = Double.parseDouble(strAry[2]) * 0.3048; }//feetをmに換算 if (!(strAry[3].equals("None"))){latitude = Double.parseDouble(strAry[3]);} if (!(strAry[4].equals("None"))){longitude = Double.parseDouble(strAry[4]);} if (!(strAry[5].equals("None"))){trackangle = Double.parseDouble(strAry[5]);} if (!(strAry[6].equals("None"))){groundspeed = Double.parseDouble(strAry[6]) * 0.514 ;}//ノットをm/sに換算 if (!(strAry[7].equals("None"))){ trueairspeed = Double.parseDouble(strAry[7]) * 0.514 ;}//ノットをm/sに換算 if (!(strAry[8].equals("None"))){heading = Double.parseDouble(strAry[8]);} if (!(strAry[10].equals("None"))){mach = Double.parseDouble(strAry[10]);} temperature = temperature(trueairspeed, mach) ; speedx = wind_x(trackangle , groundspeed, heading, trueairspeed) ; speedy = wind_y(trackangle , groundspeed, heading, trueairspeed); if (((latitude <= 36.5) && (latitude >= 34.75)) && ((longitude <= 140.5) && (longitude >= 137.5) )){//空域を限定 if (heading > 0 && trueairspeed > 0 && groundspeed > 0 && mach > 0 ){//ノイズ除去 output0.println(Time+ "," +icao +"," +altitude + "," + trackangle +"," +groundspeed +"," +heading + "," + trueairspeed + "," + temperature + "," + windspeed(speedx , speedy)+ "," + windangle( speedx , speedy)); } } } input1.close(); output0 .close(); System.out.println(filename0 + "の書き込み終了"); } catch(Exception e){ System.out.println("エラー:　" + e + " " + ss1); } } public static double temperature(double trueAirSpeed, double mach){//trueAirSpeedとmachから気温を算出 try{ double k = 1.403 ; //比熱比（空気） double M = 28.966E-3; //分子量 kg/mol　（空気） double R = 8.314472;// 気体定数 double K = M / k / R; return (trueAirSpeed /mach ) * (trueAirSpeed /mach ) * K - 273.15; }catch(Exception e){ System.out.println( e ); return 0; } } public static double wind_x(double trackangle ,double groundspeed, double heading, double trueairspeed){//風のx成分 double MAGNETIC_DECLINATION=7.64;//日本における磁気偏角 double ground_dir = 0.5 * Math.PI - Math.toRadians(trackangle) ; double ground_x = groundspeed * Math.cos(ground_dir); double air_dir = 0.5 * Math.PI - Math.toRadians(heading)+ Math.toRadians( MAGNETIC_DECLINATION ); double air_x = trueairspeed * Math.cos(air_dir); return ground_x - air_x; } public static double wind_y(double trackangle ,double groundspeed, double heading, double trueairspeed){//風のy成分 double MAGNETIC_DECLINATION =7.64; double ground_dir = 0.5 * Math.PI - Math.toRadians(trackangle) ; double ground_y = groundspeed * Math.sin(ground_dir); double air_dir = 0.5 * Math.PI - Math.toRadians(heading)+ Math.toRadians( MAGNETIC_DECLINATION); double air_y = trueairspeed* Math.sin(air_dir); return ground_y - air_y; } public static double windangle(double speedx ,double speedy){//風のx成分と風のy成分から風向を求める return Math.toDegrees((0.5 * Math.PI - Math.atan2(speedy , speedx) + 2.0 * Math.PI + Math.PI ) % (2.0 * Math.PI)) ; } public static double windspeed(double speedx ,double speedy){//風のx成分と風のy成分から風速を求める return Math.sqrt(speedx * speedx + speedy *speedy) ; } }

高層気象台との比較

　筑波館野高層気象台の計測データは気象庁ホームページで公開されています（ただし1日後に）。それと今回のデータを比較します。筑波館野高層気象台と日野は約80kmはなれていることを考慮すると、下のように温度、風速、風向ともに、よく一致しています。
ここで、筑波のデータは2024年6月11日の9時のラジオゾンデによる計測データ（最高高度に達するまで約1時間かかるので時間差がある）。当方日野のデータは、９時から９時59分まで間に受信した気象データを、標高100mごとに集計し平均化したものです。
この1時間の間に、90万メッセージを受信し、その中から5500の有意な気象データを抽出しています。一方その間のACRSからの気象データはわずか49個でした。

次のようにコメントできるでしょう。

１．ラジオゾンデは2万メートル以上の高度まで（風船が破裂するまで）計測されますが、最高高度に達するまで約1時間かかるそうです。航空機の巡航高度は最高1万2千メートル前後であり、限定されます。しかし、地表から高 さ約1万メートルまでの範囲を対流圏（troposphere）と言い、気象現象はほとんどこの中で起こっていますので充分です。
２．当地から東京国際空港（羽田）まで、３５ｋｍですが、低高度の航空機からの電波が受信できていません。そのため高度1000m以下のデータがあまり取れません。航空機の速度が低いと、得られる気温は誤差が多いようです。
３．航空機はラジオゾンデのように垂直に上昇するわけではなく、高高度では遠隔のデータを受信することがあります。広範囲の受信は必要ではなく、空域を限定すべきです。（緯度36.5〜 34.75、経度140.5〜 137.5にしています）
４．ラジオゾンデのように、1日に2回だけではなく、リアルタイムで有効な情報が得られます。
５．ただし、航空機の飛ばない深夜時間帯にはデーターがえられません。また、台風や雷雲を避けて飛行しますので、その中のデータは得られません。
６．データは、気温、風向、風速のみで、湿度はありません。

　以上のことを配慮しても、多数のデータが得られ、コストがかからないので、航空機の気象データは大いに活用されるべきと考えます。

リアルタイム表示

ナチュラル研究所のホームページではリアルタイムにModeSデータにもとずく高層気象を表示しています。高度1000m以下のデータを補完するため、ACARSデータを合わせて表示しています。さらに充実していくつもりです。
ModeSのメッセージは混雑時は衝突したり、伝送エラーのため、誤ったデータが受信される時があります。そのため現実的でない異常値が混ざる場合があり、適切なフィルター処理により、排除する必要があります。
また使用したpyModeSライブラリーはそのようなデータを受信したとき、ストール（プログラムが停止）する場合があります。回避する改良を加えたpyModeSを用意しています。問い合わせください。