Arduinoベースの気象観測所データ分析

独自の個人気象ステーションの作成は、以前よりはるかに簡単になりました。 ニューイングランドの一貫性のない天気を考慮して、独自の気象観測所を作成し、MATLABを使用して気象データを分析することにしました。



この記事では、次の質問に答えます。

• 過去3時間で風はどの方向に吹きましたか？
• 先週、気温と露点はどのように変化しましたか？
• 雷雨が近づくと、気圧は実際に低下しますか？

議論されている問題が非常に単純であることは明らかですが、説明されている手法とコマンドは、より複雑で実際的な問題の解決に役立ちます。

ステップ1：気象観測所の場所

最初に、気象ステーションの設置場所を決定する必要がありました。 駐車場の最上階で行うことにしました。 場所は、天候の影響を受けたために選択されましたが、電子機器を取り付けるための屋根もありました。 最終的にサードパーティのデータアグリゲーターにデータを提供したかったため、場所とハードウェアを選択する際にこれを考慮する必要がありました。

ステップ2：ハードウェアの選択

選択した場所は建物のWi-Fi範囲外であったため、ステーションから隣接する建物内の受信機にデータを転送する方法を見つける必要がありました。これを行うには、Arduino UnoをXBee高出力送信機に接続しました。 次に、データは建物内に既にあるArduino MegaのXBee受信モジュールに転送されます。 このArduinoはインターネットに接続され、データは無料のデータ集約サービスThingSpeak.comに送信されます。 設置に使用した機器の完全なリストを以下に示します。

設備リスト

• 気象シールドと気象センサーを備えた Arduino Uno
• Digi Internationalの2x XBeeシールドおよび2x XBee拡張範囲モジュール
• イーサネットシールド付きArduino Mega
• 2個のArduinoマウントキット
• Arduino電源用の2x 5Vトランス
• XBeeトランスミッターをプログラミングするためのXBee Explorerドングル

ステップ3：気象ステーションの送信機を接続し、屋外のArduinoをプログラムする

外部のArduinoはガレージに設置され、測定値を収集し、屋内のArduinoにデータを送信します。 この交換を実装するために、WeatherシールドとXBeeシールドのコンタクトを最初に接続しました。 次に、シールドに付属のドキュメントに従って、XBeeシールドをArduino Unoに接続しました。 XBeeシールドには、XBee高出力トランシーバーがあります。 X-CTUソフトウェアを使用して、目的の受信者アドレスをプログラムし、目的のファームウェアをXBeeシールドにロードしました。 風速計、雨風センサーは、付属のRJ-45コネクタを使用してウェザーシールドに接続されました。

ステップ4：Weather Station Receiverと屋内Arduinoプログラミングを接続する

内部Arduinoは建物内にあり、外部Arduinoからデータを受信し、データを検証し、ThingSpeakに送信する役割を担います。ドキュメンテーションによるMega ArduinoのXBeeシールド。 再びX-CTUを使用して、目的の受信者アドレスをプログラムし、目的のファームウェアをXBeeシールドにロードしました。



次に、Arduinoをプログラムして、XBeeメッセージを受信し、パケットを1分間に1回の割合でThingSpeakに転送します。 ThingSpeakメッセージを送信する前に、アカウントを設定し、チャネルと場所の情報を構成しました。 データ転送のセットアップが終了するとすぐに、MATLABでデータの分析を開始しました。



この記事に記載されている分析を再現するために、機器を購入して構成する必要はありません、ステーションの現在のデータはチャンネル12 397で入手できます。

ステップ5：質問に答える

ThingSpeakから気象データを取得する

最初の2つの質問に答えるには、次のコマンドを使用します

thingSpeakFetch
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

利用可能なデータフィールドを表示するには、すべてのフィールドを同時にMATLABにインポートし、温度、湿度、風速、風向に関するデータを変数に保存します。 ThingSpeakサポートドキュメント 。

 [d,t,ci] = thingSpeakFetch(12397,'NumPoints',8000); % fetch last 8000 minutes of data
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

8000ポイントは、ThingSpeakが一度にリクエストできる最大ポイント数です。 測定頻度では、これは約6日間の測定に相当します。

 tempF = d(:,4); % field 4 is temperature in deg F baro = d(:,6); % pressure in inches Hg humidity = d(:,3); % field 3 is relative humidity in percent windDir = d(:,1); windSpeed = d(:,2); tempC = (5/9)*(tempF-32); % convert to Celsius availableFields = ci.FieldDescriptions'
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

 availableFields = 'Wind Direction (North = 0 degrees)' 'Wind Speed (mph)' '% Humidity' 'Temperature (F)' 'Rain (Inches/minute)' 'Pressure ("Hg)' 'Power Level (V)' 'Light Intensity'
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

観測期間の基本的な統計の計算とその視覚化

データをよりよく理解するために、インポートしたデータの最小値、最大値、平均値を見つけ、最大値と最小値の時間を見つけます。これにより、気象観測所からのデータを迅速に検証できます。

 [maxData,index_max] = max(d); maxData = maxData'; times_max = datestr(t(index_max)); times_max = cellstr(times_max); [minData,index_min] = min(d); minData = minData'; times_min = datestr(t(index_min)); times_min = cellstr(times_min); meanData = mean(d); meanData = meanData'; % make column vector summary = table(availableFields,maxData,times_max,meanData,minData,times_min) % display
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

 summary = availableFields maxData times_max ____________________________________ _______ ______________________ 'Wind Direction (North = 0 degrees)' 338 '10-Jul-2014 05:01:32' 'Wind Speed (mph)' 6.3 '10-Jul-2014 12:47:14' '% Humidity' 86.5 '15-Jul-2014 04:51:24' 'Temperature (F)' 96.7 '12-Jul-2014 16:28:55' 'Rain (Inches/minute)' 0.04 '15-Jul-2014 13:47:13' 'Pressure ("Hg)' 30.23 '11-Jul-2014 09:25:07' 'Power Level (V)' 4.44 '10-Jul-2014 10:25:01' 'Light Intensity' 0.06 '12-Jul-2014 13:23:38' meanData minData times_min _________ _______ ______________________ NaN 0 '10-Jul-2014 04:54:32' 3.0272 0 '10-Jul-2014 01:33:14' 57.386 25.9 '12-Jul-2014 13:39:39' 80.339 69.6 '11-Jul-2014 06:59:54' 5.625e-05 0 '10-Jul-2014 01:02:11' 30.04 29.78 '15-Jul-2014 13:04:08' 4.4149 4.38 '11-Jul-2014 09:22:06' 0.0092475 0 '10-Jul-2014 01:02:11'
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

最大気圧40気圧や最大温度1700度などの予期しない値が得られた場合は、データが間違っていると想定できます。 このようなエラーは、伝送エラー、電圧サージなどの理由で発生する可能性があります。 記事の最後に、このような「外れ値」を処理する方法をいくつか示しますが、アップロードされたデータについては、このレポートが公開されたときにすべてが正常に見えます。



上記のコードでは、MATLAB R2013bで最初に導入された表形式のデータ型を使用しました。このデータ型の詳細については、 ブログエントリを参照してください。

過去3時間の風の可視化

気象ステーションは1分に1回程度の気象レポートを受信するため、過去180分を見て過去3時間の風についての質問に答え、MATLABのコンパスプロットを使用して速度と方向を確認します。興味深い時に風。 これは、北（0度）が右側にあり、度が反時計回りに増加する数学コンパスです。 90度は東（上）、180度は南（左）、270度は西（下）を表します。 雷雨または正面通過中にウィンドシアがない場合、コンパスは、原則として、グラフ上の高密度の矢印で示される、風の優先方向を示します。 この例で要求されたデータの場合、主にニューイングランドの夏時間に典型的な南および南西から風の方向を観察します。

 figure(1) windDir = windDir((end-180):end); % last 3 hours windSpeed = windSpeed((end-180):end); rad = windDir*2*pi/360; u = cos(rad) .* windSpeed; % x coordinate of wind speed on circular plot v = sin(rad) .* windSpeed; % y coordinate of wind speed on circular plot compass(u,v) titletext = strcat(' for last 3 hours ending on: ',datestr(now)); title({'Wind Speed, MathWorks Weather Station'; titletext})
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

露点計算

これで、温度と露点がどのように変化したかについての2番目の質問に答える準備ができました。

先週。露点は、空気（冷却した場合）が水蒸気で飽和する温度です。 気団の湿度が高いほど、露点が高くなります。 露点は、不快感の尺度として使用されることもあります。 露点が65度（+ 18C）を超えると、多くの人が空気が「べたつく」と言い始めます。 70の露点では、多くの人が不快に感じます。 以下に示すように、方程式と定数を使用して、露点の一般的な推定値TDPを見つけることができます。



どこで



と

さて、上記の方程式をMATLABコードとして書いたので、次のことができます。

 b = 17.62; c = 243.5; gamma = log(humidity/100) + b*tempC ./ (c+tempC); tdp = c*gamma ./ (b-gamma); tdpf = (tdp*1.8) + 32; % convert back to Fahrenheit
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

温度、空気湿度、露点対時間の視覚化

露点を計算したので、データを視覚化し、過去5日間または6日間の挙動を観察する準備ができました。

 figure(2) [ax, h1, h2] = plotyy(t,[tempF tdpf],t,humidity); set(ax(2),'XTick',[]) set(ax(2),'YColor','k') set(ax(1),'YTick',[0,20,40,60,80,100]) set(ax(2),'YTick',[0,20,40,60,80,100]) datetick(ax(1),'x','keeplimits','keepticks') set(get(ax(2),'YLabel'),'String',availableFields(3)) set(get(ax(1),'YLabel'),'String',availableFields(4)) grid on legend('Location','South','Temperature','Dew point', 'Humidity')
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

典型的な週の間に、温度と湿度の毎日の変化をはっきりと見ることができます。 予想通り、相対湿度は夜（温度が露点まで下がる）に上昇する傾向があり、通常は午後に最高温度に達します露点温度は気団の湿度を示します。 この例の公開用にデータを収集すると、露点が70度を超えていることがわかります。これは、ニューイングランドの夏の暑くて蒸し暑い日の典型的な例です。気象観測所によって報告された最新のデータをどのように要求するかについて。

ニューイングランドの雨の日の気象データの受信と処理

私たちが答えたかった最後の質問は-雨の前に気圧が本当に下がるのか？ これを行うために、有名な雨の日の気象観測所からデータを受け取りました。 今回は、気圧と降水量に興味があります。 満杯になると空になる自己排出センサーを使用しました。 0.01インチの降水量が降ると、センサーが回転して空になり、Arduinoコードが1分あたりの空の数をカウントし、適切な降雨値をThingSpeakに送信します。 MATLABを使用してデータから1時間ごとのサンプルをサンプリングし、降水量と圧力の傾向の蓄積を簡単に確認できるようにします。

 [d,t,ci] = thingSpeakFetch(12397,'DateRange',{'6/4/2014','6/6/2014'}); % get data baro = d(:,6); % pressure extraData = rem(length(baro),60); % computes excess points beyond the hour baro(1:extraData) = []; % removes excess points so we have even number of hours rain = d(:,5); % rainfall from sensor in inches per minute
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

6月5日に本当に雨が降っていましたか？ 最大値は0.01インチしかないため、データ転送の1分ごとを見るだけでは言いにくいです。 しかし、6月5日のすべての降雨量を要約すると、0.48インチの雨が降ったことがわかります。これは、平均月間3.68インチの13％で、この日は本当に雨が多いことを示しています。 最大降雨量がいつ低下したかをよりよく理解するために、以下に示すように、データを時間ごとのデータに変換しました。

 rain(1:extraData) = []; t(1:extraData) = []; rainHourly = sum(reshape(rain,60,[]))'; % convert to hourly summed samples maxRainPerMinute = max(rain) june5rainfall = sum(rainHourly(25:end)) % 24 hours of measurements from June 5 baroHourly = downsample(baro,60); % hourly samples timestamps = downsample(t,60); % hourly samples
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

 maxRainPerMinute = 0.0100 june5rainfall = 0.4800
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

時間ごとのクラウドおよび気圧データの視覚化と気圧傾向の検出

データを前処理した後、それらを視覚化する準備が整いました。 ここでは、MATLABを使用して、気圧データに対する傾向線を見つけます。 プロットすると、豪雨の前に圧力が低下しますか？

 figure(3) subplot(2,1,1) bar(timestamps,rainHourly) % plot rain xlabel('Date and Time') ylabel('Rainfall (inches /per hour)') grid on datetick('x','dd-mmm HH:MM','keeplimits','keepticks') title('Rainfall on June 4 and 5') subplot(2,1,2) hold on plot(timestamps,baroHourly) % plot barometer xlabel('Date and Time') ylabel(availableFields(6)) grid on datetick('x','dd-mmm HH:MM','keeplimits','keepticks') detrended_Baro = detrend(baroHourly); baroTrend = baroHourly - detrended_Baro; plot(timestamps,baroTrend,'r') % plot trend hold off legend('Barometric Pressure','Pressure Trend') title('Barometric Pressure on June 4 and 5')
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

このデータを視覚化し、傾向を確認すると、降水前に気圧が本当に低下していることがはっきりとわかります。

無効な測定値からデータを消去する

質問は非常に簡単でしたが、利用可能なデータに問題がある場合があります。 たとえば、5月30日に温度に関するいくつかの誤ったデータを記録しました。 MATLABを使用してフィルター処理する方法を見てみましょう。

 [d,t] = thingSpeakFetch(12397,'DateRange',{'5/30/2014','5/31/2014'}); rawTemperatureData = d(:,4); newTemperatureData = rawTemperatureData; minTemp = min(rawTemperatureData) % wow that is cold!
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

 minTemp = -1.7662e+03
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

しきい値フィルターを使用してデータスパイクを削除する

しきい値テストに合格しなかったアイテムを削除します。 この場合、華氏-1766度など、明らかに信頼できない値がいくつかあります。 したがって、0〜120の温度値のみを含むデータを使用できます。これは、ニューイングランドの春の日に許容される値です。

 tnew = t'; outlierIndices = [find(rawTemperatureData < 0); find(rawTemperatureData > 120)]; tnew(outlierIndices) = []; newTemperatureData(outlierIndices) = [];
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

クリーンアップされたソースデータをプロットします。

 figure(4) subplot(2,1,2) plot(tnew,newTemperatureData,'-o') datetick xlabel('Time of Day') ylabel(availableFields(4)) title('Filtered Data - outliers deleted') grid on subplot(2,1,1) plot(t,rawTemperatureData,'-o') datetick xlabel('Time of Day') ylabel(availableFields(4)) title('Original Data') grid on
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

メディアンフィルターを使用して無効なデータを削除する

不良データを削除する別の方法は、メディアンフィルターを適用することです。 メディアンフィルターは、データセットに関する知識をそれほど必要としません。 最近傍の中央の外側にある値を単純に削除します。 結果のフィルタリングは、データポイントの削除とは対照的に、元と同じ長さのベクトルであり、データギャップとレコードの短縮につながります。 このタイプのフィルターは、信号からノイズを除去するためにも使用できます。

 n = 5; % this value determines the number of total points used in the filter
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

nの大きな値は、比較のための「隣接」の数を示します。 1分に1回収集された温度では、n = 5を選択します。これは、通常、温度が5分以内に大幅に変化しないためです。

 f = medfilt1(rawTemperatureData,n); figure(5) subplot(2,1,2) plot(t,f,'-o') datetick xlabel('Time of Day') ylabel(availableFields(4)) title('Filtered Data - using median filter') grid on subplot(2,1,1) plot(t,rawTemperatureData,'-o') datetick xlabel('Time of Day') ylabel(availableFields(4)) title('Original Data') grid on
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

おわりに

そこで、気象観測所からのデータを分析し、質問に答えました。

また、初期検証に合格しないデータのフィルタリング方法もいくつか示しました。



この記事で提起された質問は初歩的ですが、何千ものステーションからデータを収集することで何が達成できるかを考えてください。 さらに、産業システムおよび科学システムのリアルタイムの変化を測定および分析し、データに基づいて意思決定を行うことができます。



この記事に記載されている分析を再現するために、機器を購入して構成する必要はありません。 インストールのライブデータは、ThingSpeak.comのチャンネル12 397で入手できます。



MATLABプロジェクトコード