Как мне запустить функцию Торнтвейта (которая вычисляет стандартный индекс осадков) по местоположению и широте в R?

У меня есть следующие данные:

dat <- read.table(text="
  id YEAR    MONTH   TMED    PRCP    lat
   1  1986    1      -14.5    2.3    42.4863
   1  1986    2      -13.9    5.7    42.4863
   2  1986    1      -12.9    7.2    42.46
   2  1986    2      -11.6    19.7   42.46", header=TRUE)

куда

  • id - это моя единица местоположения в диапазоне от 1 до 90
  • TMED - среднемесячная температура по местности
  • PRCP - осадки по местоположению
  • lat - широта по местоположению
  • ГОД колеблется с 1986 по 2016 год.
  • МЕСЯЦ колеблется от 1 до 12

Мне нужно запустить следующую функцию в R, чтобы рассчитать индекс SPEI (стандартный индекс осаждения-испарения) для каждого местоположения:

library(SPEI)

dat$PET <- thornthwaite(dat$TMED, dat$lat[1])

dat$BAL <- dat$PRCP-dat$PET

spei1 <- spei(dat$BAL, scale = 1)

dat$spei1=l

Этот код работает для одного места. Но мне нужно сделать петлю по широте и местоположению. Одна из проблем заключается в том, что широта должна входить в функцию thornthwaite как число (а не как список / переменную).


person Anastasia    schedule 01.10.2018    source источник
comment
Пожалуйста, не выкладывайте скриншоты, они бесполезны. Вместо этого предоставьте пример данных в виде кода, чтобы другие могли его воспроизвести и помочь.   -  person    schedule 01.10.2018
comment
Выполнено. Спасибо, @Pearly Spencer.   -  person Anastasia    schedule 01.10.2018
comment
thornthwaite() тот, что находится в SPEI библиотеке?   -  person AkselA    schedule 01.10.2018
comment
Да, это библиотека SPEI @AkselA   -  person Anastasia    schedule 01.10.2018


Ответы (1)


Я немного поигрался с уравнением Торнтвейта для своих мастеров по экологии, и эта реализация кажется немного странный. Несмотря на то, как это может выглядеть здесь, уравнение требует в качестве входных данных больше, чем просто среднюю температуру и широту. На самом деле ему нужна средняя продолжительность светового дня в данном месяце, но ее можно вычислить по широте и дате, и thornthwaite() получает дату, просто предполагая, что первая точка данных представляет январь, а остальные следуют последовательно. Уравнение Торнтуэйта также зависит от годового индекса тепла, что означает, что вам нужны среднемесячные значения температуры за весь год. thornthwaite() решает эту проблему путем агрегирования по заданному вами вектору температуры.

В итоге, для thornthwaite() работы вам потребуется последовательность значений среднемесячных температур, начиная с января и охватывая не менее одного года. Таким образом, функция не будет работать с предоставленными вами данными.

Я предлагаю вам убедиться, что ваша серия достаточно длинная, а также разбить ее на отдельные data.frames для каждого местоположения. Для этого можно использовать split() (split(dat, dat$id)).

В ?thornthwaite есть несколько примеров, в том числе один, демонстрирующий его использование для временных рядов, что полезно, если ваша серия не начинается в январе.

Я сделал макет, демонстрирующий один из возможных подходов:
(Обратите внимание, что функция будет возвращать значения, даже если данные не охватывают полный год, тогда эти значения будут весьма ненадежными.)

dat <- read.table(text="
  id YEAR    MONTH   TMED    PRCP    lat
   1  1986    1     -14.5    2.3    42.4863
   1  1986    2     -13.9    5.7    42.4863
   1  1986    3     -10.5    2.3    42.4863
   1  1986    4      -7.9    5.7    42.4863
   1  1986    5      -4.5    2.3    42.4863
   1  1986    6       0.9    5.7    42.4863
   1  1986    7      10.5    2.3    42.4863
   1  1986    8      17.9    5.7    42.4863
   2  1986    1     -12.9    7.2    42.46
   2  1986    2     -11.6   19.7    42.46
   2  1986    3      -8.9    7.2    42.46
   2  1986    4      -5.9    7.2    42.46
   2  1986    5       1.6   19.7    42.46
   2  1986    6      12.9    7.2    42.46
   2  1986    7      21.6   19.7    42.46
   2  1986    8      25.6   19.7    42.46", header=TRUE)

dat.s <- split(dat, dat$id)

lapply(dat.s, function(x) thornthwaite(x$TMED, x$lat[1]))
person AkselA    schedule 01.10.2018
comment
Спасибо за ответ @AksleA! Да, у меня есть все месяцы на каждый год - я был здесь не для экономии места. Я пробовал использовать функцию разделения, как вы предлагаете, но я не знаю, как включить в нее изменяющуюся широту. Вы бы порекомендовали, возможно, просто использовать вместо этого набор данных с координатной привязкой SPEI (уже рассчитанный авторами этого индекса)? Думаю, я могу извлечь значения на основе моей широты и долготы. Спасибо еще раз! - person Anastasia; 01.10.2018
comment
OK. Но эта команда lapply (dat.s, function (x) thornthwaite (x $ TMED, x $ lat [1])) берет первое значение широты из разделения или из исходного фрейма данных? - person Anastasia; 01.10.2018
comment
@ Анастасия: Из раскола. split() возвращает список data.frames, lapply() проходит по этому списку, применяя thornthwaite() к каждому data.frame. Уравнение Торнтуэйта даст вам лишь очень грубое приближение к потенциальной эвапотранспирации. Если у вас есть доступ к более надежным данным, вам обязательно стоит пойти на это. - person AkselA; 01.10.2018
comment
@ Анастасия: Нет проблем. Когда вы найдете ответ полезным, вы можете поставить ему зеленую галочку, щелкнув символ в левом верхнем углу под стрелками голосования. - person AkselA; 01.10.2018