Вычислить выделенную область из сложенного растрового объекта

Есть несколько методов доступа к слоям и значениям стека растра. Я предпочитаю следующее:

#select first layer of raster stack 
stk[[1]]

#get values of second layer
stk[[2]][]

Теперь вернемся к вашему вопросу:

Когда я хочу рассчитать площадь пикселей, соответствующих определенным критериям, я делаю следующее (при работе с небольшими растрами):

numberOfPixels <- sum(stk[[1]][] <= -1000 & stk[[1]][] >= -4000 & stk[[2]][] >= 10, na.rm=T)

Это дает вам количество пикселей, соответствующих определенным критериям. Если вы будете работать в системе координат проекции (вы работаете в WGS 84 и, следовательно, вы не можете рассчитать точные площади из этого), вы просто умножите numberOfPixels на разрешение вашего растра:

area <- numberOfPixels * (res(stk)[1] * res(stk)[2])

Если вы хотите получить площадь в квадратных метрах, перепроецируйте свой растр в систему координат проекции. Например, в UTM. В вашем случае этот вариант может подойти (обратите внимание, что ваш экстент охватывает несколько зон): http://www.spatialreference.org/ref/epsg/wgs-84-utm-zone-24s/

stk <- projectRaster(from=stk, crs="+proj=utm +zone=24 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

Затем еще раз:

numberOfPixels <- sum(stk[[1]][] <= -1000 & stk[[1]][] >= -4000 & stk[[2]][] >= 10, na.rm=T)
area <- numberOfPixels * (res(stk)[1] * res(stk)[2])
area
[1] 258898897920 

Хорошо, вот возможный другой подход, основанный на преобразовании ячеек растра в полигоны и последующем вычислении площади полигонов с использованием пакета geosphere:

require(geosphere)
polys = rasterToPolygons(stk)
polys_sub = subset(polys, (layer.1 <= -1000 & layer.1 >= -4000 & layer.2 >= 10))
> sum(areaPolygon(polys_sub))/1E6  #in km2

[1] 157035.2

Это очень похоже на то, что OP получил от метода cellStats (скидка около 2%), и все же полностью отличается от метода reprojection. Насколько я понимаю из areaPolygon, это должен быть численно правильный подход. Однако я до сих пор не понимаю различных результатов в отношении ответа @maRtin.

Кстати, для OP: зачем вы настраиваете пример с такой "странной" системой отсчета (широта / долгота, другое разрешение x / y)?

Итак, сначала пара проверок на вменяемость ...

В какой области мы имеем дело? Мы основываем это на вашем стеке и будем использовать res() для получения разрешений. Достаточно просто взять их вручную, но я предпочитаю называть их.

 length(stk)*res(stk)[1]*res(stk)[2]

 70.32058

Таким образом, я получаю ~ 70,3 как фактическую площадь стека. Я получаю это только из того, что вы предоставили для данных, поэтому, если CRS или прогнозы ненадежны, это на вас. Если мы проведем работу по поиску области для нашего подмножества и обнаружим, что она сильно отличается от этой, мы сбились с пути.

Теперь я начну с вашего 4-го блока кода, в котором вы это сделали:

 a <- stk[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >= 10]

Итак, теперь у нас есть слой a, на котором мы разделили stk, так что теперь это список значений, соответствующих критериям. Нам не нужен этот список, мы просто хотим знать, сколько элементов он содержит.

Заменить:

 a<-length(stk[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >= 10])

поскольку мы только хотим знать, сколько элементов находится в списке ...

теперь, используя res() для получения разрешений x и y:

 a*res(stk)[1]*res(stk)[2]

Результаты в:

 29.80777

Судя по площади растра, которую мы вычислили ранее, это кажется разумным, ~ 42% от общей площади. Мне это кажется разумным.

Что касается проблем с проекциями и системами координат в других комментариях к ответам, ДА, это проблема. Но вы должны убедиться, что ваши прогнозы не ошибочны и имеют смысл для того, что вы пытаетесь сделать. Также помните, что простое переназначение системы проецирования - это не то же самое, что повторное проецирование. Вероятно, одна из самых распространенных ошибок, которые я вижу при работе с растрами в R, - это просто переназначение CRS вместо перепроецирования растра. Вы должны знать, в какой форме поступают ваши данные и что вам нужно с ними делать, чтобы делать с ними математические расчеты.