Создание собственной личной метеостанции стало намного проще, чем раньше. С учетом непостоянной погоды в Новой Англии, мы решили, что хотим создать нашу собственную метеостанцию и использовать MATLAB для анализа метеоданных.

В статье мы ответим на следующие вопросы:

• В каком направлении дул ветер в течение последних 3-х часов?
• Как изменялись температура и точка росы в течение последней недели?
• На самом ли деле падает барометрическое давление при приближении грозы?

Понятно, что рассмотренные вопросы достаточно просты, но описанные приемы и команды помогут вам решать более сложные практические задачи.

Шаг 1: Размещение метеостанции

Во-первых, мы должны были решить, где разместить нашу метеостанцию. Мы решили сделать это на верхнем этаже автостоянки. Место было выбрано, поскольку оно было под воздействием погодных явлений, но также и имело крышу для монтажа электроники. Из-за того, что мы в конечном счете хотели отдавать данные стороннему агрегатору данных, мы должны были учитывать это при выборе места и железа.

Шаг 2: Выбор аппаратного обеспечения

Выбранное расположение было вне радиуса действия Wi-Fi нашего здания, поэтому мы должны были найти способ передачи данных от станции к приемнику, который находился внутри соседнего здания.Чтобы сделать это, мы подсоединили Arduino® Uno к передатчику высокой мощности XBee®. Затем данные передаются в принимающий модуль XBee на Arduino Mega уже внутри здания. Этот Arduino был подключен к Интернету и данные отправляются в бесплатный сервис агрегации данных, ThingSpeak.com. Полный список аппаратуры, которую мы использовали в нашей установке приведен ниже.

Список оборудования

• Arduino Uno с Weather shield и датчиками погоды
• 2x XBee Shield и 2x модулей XBee с расширенным диапазоном от Digi® International
• Arduino Mega с Ethernet Shield
• 2 комплекта для монтажа Arduino
• 2x 5В трансформатор для питания Arduino
• XBee Explorer Dongle для программирования передатчика XBee

Шаг 3: Подключение передатчика метеостанции и программирование наружной Arduino

Внешнее Arduino установлено на гараже и отвечает за сбор измерений и передачу данных на Arduino в помещении. Для реализации такого обмена, сначала мы соединили контакты от Weather shield и XBee shield. Затем мы соединили XBee shield с Arduino Uno в соответствии с документацией, предоставленной с шилдами. В XBee shield есть XBee трансивер высокой мощности. Мы использовали программное обеспечение X-CTU для программирования нужного адреса получателя и загрузки нужной прошивки в  XBee shield. Анемометр, датчики ветра дождя были подключены к weather shield с помощью разъемов RJ-45, входящими в комплект поставки.

Шаг 4: Подключение приемника погодной станции и программирование Arduino, находящегося в помещении

Внутреннее Arduino находится внутри нашего здания, и отвечает за получение данных от наружного Arduino, за проверку достоверности данных, и их отправку на ThingSpeak.Чтобы построить такую систему, сначала мы припаяли контакты на Ethernet shield и второй XBee shield.Затем мы вставили Ethernet shield и XBee shield в Мега Arduino в соответствии с документацией. Мы вновь использовали X-CTU для программирования необходимого адреса получателя и загрузки нужной прошивку на XBee shield.

Далее, мы запрограммировали Arduino получить XBee сообщения и пересылать пакеты на ThingSpeak со скоростью один раз в минуту. Перед отправкой сообщения ThingSpeak, мы настроили учетную запись и сконфигурировали информацию о канале и расположении. Как только закончилась настройка передачи данных, мы начали анализировать данные в MATLAB.

Обратите внимание, что для воспроизведения анализа, приведенного в статье, вам не нужно покупать и настраивать оборудование.Текущие данные нашей станции доступны на канале 12 397.

Шаг 5: Отвечаем на вопросы

Получение данных о погоде из ThingSpeak

Чтобы ответить на наши первые два вопроса, мы используем команду

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.См запись в блоге для дальнейшего ознакомления с этим типом данных.

Визуализация ветра за последние три часа

Так как наша погодная станция получает метеорологические сводки примерно раз в минуту, мы смотрим на последние 180 минут, что бы ответить на наш вопрос о том, какой ветер был в последние 3 часа, и мы используем Compass plot в 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})

Расчет точки росы

Теперь мы готовы ответить на наш второй вопрос о том, как температура и точка росы менялись в течении
прошлой недели.Точка росы — это температура, при которой воздух (при охлаждении) становится насыщенным водяным паром. Чем более влажность воздушной массы, тем выше точка росы. Точка росы также иногда используется как мера дискомфорта. Когда точка росы превышает 65 градусов (+18С), многие люди начинают говорить, что воздух «липкий». При точке росы за 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 градусов, что является типичной для жаркого и влажного летнего дня в Новой Англии.Если вы запустите этот код в MATLAB, вы можете получить разные ответы, по мере того, как вы будете запрашивать самые последние данные, сообщенные метеостанцией.

Получение и обработка данных о погоде за дождливый день в Новой Англии

Последний вопрос, на который мы хотели ответить — действительно ли атмосферное давление падает перед дождем? Чтобы сделать это, мы получили данные из метеостанции за известный дождливый день. На этот раз, мы заинтересованы в барометрическом давлении и количестве осадков. Мы использовали датчик самостоятельного опорожнения, который опорожняется, когда он заполнен. Наш датчик вращается и опорожняется, когда выпало 0,01 дюймов осадков.Наш код Arduino считает количество опорожнений за каждую минуту и передает соответствующее значение осадков в ThingSpeak. Мы используем MATLAB для выборки из наших данных почасовых образцов, чтобы мы могли легко увидеть накопление осадков и тенденции давления.

[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

Действительно ли шел ливень 5 июня? Это трудно сказать, если мы просто посмотрим на каждую минуту передачи данных, так как максимальное значение составляет всего 0,01 дюйма. Однако, если суммировать все осадки за 5 июня, мы видим, что мы получили 0,48 дюйма осадков, что составляет 13% от среднемесячного 3,68 дюйма, что показывает, что день был действительно очень дождливым. Чтобы получить более полное представление о том, когда выпал максимум осадков, мы преобразовали данные в часовые данные, как показано ниже.

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')

После того, как мы визуализируем эти данные и посмотрим на тренд, мы ясно видим, что атмосферное давление действительно падает ровно до начала осадков!

Очистка данных от недостоверных измерений

Наши вопросы были довольно простыми, но иногда возникают проблемы с имеющимися данными. Например, 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 обозначают количество «соседей» для сравнения. При температурах, собранных раз в минуту, мы выбираем 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

Вывод

Итак, мы проанализировали данные с нашей метеостанции и ответили на интересующие нас вопросы.
Мы также продемонстрировали некоторые способы фильтрации данных, которые не проходят первоначальную проверку.

Вопросы, поставленные в статье, элементарны, но подумайте чего можно добиться собирая данные с тысячи станций? Кроме этого можно измерять и анализировать real-time изменения промышленных и научных систем и принимать решения на основе данных.

Напоминаем, что для воспроизведения анализа, приведенного в статье, вам не нужно покупать и настраивать оборудование. Оперативные данные из нашей установки доступны на канале 12 397 на ThingSpeak.com.

MATLAB код проекта