Температура цветения

Слишком долго глядел!

К вишневым цветам незаметно

Я прилепился душой.

Облетели... осталась одна

Печаль неизбежной разлуки

Сайгё

Обнаружил интересный подход в получении информации, можно сказать поэтический. Известно, что вишня в Японии цветет в конце марта - апреле - в начале мая, когда температура становится +7 и выше. Дни цветения вишни для японцев являются днями особенными (в США тоже отмечают. Похоже, благодаря японцам). Время Cherry blossom festival японцы и гости Японии проводят в парках, наслаждаясь тихой красотой цветущей природы (ханами - любование цветами). Этому празднику тысячи лет. Столько же лет поэтические натуры с чашкой сакэ воспевают цветение. Многие оставили свои впечатления в дневниках вместе с датами, например,

Tokistune Hiramatsu on April 14, 1644: “In Seiryoden Palace, Kyoto, we enjoyed watching cherry blossoms and took sake provided by the emperor.”

Даты цветения, оставленные аккуратными предками, стали исходными данными для косвенного определения температуры в марте на тысячу лет в прошлое. В статье "Climate Change and Cherry Tree Blossom Festivals in Japan" Richard Primack и Hiroyoshi Higuchi приводят результаты трансляции даты цветения вишни в среднюю температуру в марте в Киото. На графиках: верхний - день (в году) полного цветения, нижний - средняя температура в марте, от года 1000 до 2000-го.

В целом, корреляция явно не видна, за исключением последних 200 лет, для которых видно, что температура возрастает и вишня цветет раньше. 
Не понятно, как вычислялась температура. Не имея числовых данных для графиков, написал автору (Primack) с просьбой выслать 2 таблички, но не получил их, хотя ответил на его вопросы: “кто” и “для чего”. Пришлось оцифровать оба графика в матлабе - на каждой различимой точке жмешь кнопку, получаешь координаты и сохраняешь их, затем шкалируешь на данный диапазон по x и y.
Оцифровали, вычисляем корреляцию, но сначала надо выравнять точки по годам. Для этого создаем последовательность дат (годы, “круглые” значения в отличие от отсканированных) и отображаем на нее (аппроксимируем) точки с днем цветения и точки температуры в марте.

% read scanned data
% the blossom dates in DOY in Kyoto
gd = csvread('graphDay.csv');
% derived T Co in March in Kyoto
gt = csvread('graphTemper.csv');
% a sequence of years to interpolate to
xapr = 1000:5:1990;
gdxapr = interp1(gd(:,1),gd(:,2),xapr');
gtxapr = interp1(gt(:,1),gt(:,2),xapr');
% is there any correlation?
corrcoef(gdxapr,gtxapr)

Получаем коэффициент корреляции (-0.3939), слабая корреляция есть. Минус отражает тот факт, что чем выше температура в марте, тем раньше день цветения в году.

Как же ученые получили нижний график ("Using a complicated equation")? Что за функция, которая транслирует день цветения в среднюю температуру месяца? Можно найти аппроксимирующий полином, но: физической зависимости он не объяснит, и не подтвердит, что отображение правильное. Лучше взять официальные данные о температуре в марте в Киото (имеем, начиная с 1881 года) как калибровочные, и аппроксимировать зависимость между правдивыми данными, то есть, метеоданными и данными о дне цветения из дневников. Затем использовать полученный полином для аппроксимации температурных значений на все тысячелетие.
% read meteo data and interpolate

gt1881 = csvread('kyototemper1881_2000.csv');
gdxapr1881 = interp1q(gd(:,1),gd(:,2),gt1881(:,1));
X = gdxapr1881(1:end-3);

% input data matrix

X = [ones(length(X),1) X X.^2];

% compute polynomial koefficients

theta = inv(X'*X)*X'*gt1881(1:end-3,2);
X = gd(:,2);
X = [ones(length(X),1) X X.^2];

% compute approxomated temperature values

y = X*theta;

Картинки не сильно отличаются. Возможно, авторы применили тот же прием, но с более осмысленной зависимостью по сравнению с параболой. Осталось не ясным.

Теперь проверим опасения глобального потепления. Главный исследователь темы доктор Аоно утверждает, что наблюдаемое им увеличение средней температуры в Киото за последние 170 лет (+3.4С) определяется глобальным потеплением (+0.6С) и урбанизацией региона (+2.8С). 
Проверяем. Действительно, средняя температура за последние 100 лет в Киото значительно возросла.
% plot meteo data and trend

plot(gt1881(:,1),gt1881(:,2),'-b.');
X = gt1881(:,1); X = [ones(length(X),1) X];

% koefficients of the trend line

theta = inv(X'*X)*X'*gt1881(:,2)
hold on;
plot(gt1881(:,1),X*theta,':r');

Интерес к поэзии вечен. Ставьте даты под своими стихами.