當應用於柵格堆棧時，柵格計算僅返回一個圖層

我有從日常氣候數據創建的柵格堆棧。可以在這裏找到：當應用於柵格堆棧時，柵格計算僅返回一個圖層

#!/bin/bash 
wget -nc -c -nd http://northwestknowledge.net/metdata/data/tmmx_1982.nc

目標是從這些日常記錄中獲得每月95％的溫度值。每當我使用raster包中的calc時，它只返回一層而不是12（例如，12個月）我錯過了什麼？

library(raster) 
library(tidyverse) 
library(lubridate) 

file = "../data/raw/climate/tmmx_1982.nc " 
rstr <- raster(file)

> rstr class : RasterBrick dimensions : 585, 1386, 810810, 366 (nrow, ncol, ncell, nlayers) resolution : 0.04166667, 0.04166667 (x, y) extent : -124.793, -67.043, 25.04186, 49.41686 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 data source : in memory names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ... min values : 1.3268673, 0.7221603, 1.8519223, 1.6214808, 0.8629752, 1.1126643, 1.8769895, 0.9587604, 1.7360761, 2.1099827, 2.1147265, 1.8696048, 1.7619936, 2.0253942, 2.6840794, ... max values : 73.20462, 60.35675, 64.68890, 53.11994, 60.15675, 55.91125, 77.29095, 64.39179, 48.26004, 64.70559, 79.85970, 62.31242, 53.89768, 52.15949, 80.23198, ...

date_seq <- date_seq[1:nlayers(rstr)] 
month_seq <- month(date_seq) 

mean_tmp <- stackApply(rstr, month_seq, fun = mean)

> mean_tmp class : RasterBrick dimensions : 585, 1386, 810810, 12 (nrow, ncol, ncell, nlayers) resolution : 0.04166667, 0.04166667 (x, y) extent : -124.793, -67.043, 25.04186, 49.41686 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 data source : /tmp/RtmpYf4pQe/raster/r_tmp_2017-09-25_182536_48012_88372.grd names : index_1, index_2, index_3, index_4, index_5, index_6, index_7, index_8, index_9, index_10, index_11, index_12 min values : 4.586111, 5.656802, 6.444234, 6.875973, 6.281896, 4.495534, 5.081545, 4.396824, 4.316368, 6.413400, 4.233641, 3.119827 max values : 49.12178, 47.61632, 44.70796, 47.57829, 46.97714, 51.61986, 37.77228, 51.30043, 42.51572, 36.86453, 37.96615, 52.15552

mean_90thtmp <- calc(mean_tmp, forceapply = TRUE, 
       fun = function(x) {quantile(x, probs = 0.90, na.rm = TRUE) })

> mean_90thtmp class : RasterLayer dimensions : 585, 1386, 810810 (nrow, ncol, ncell) resolution : 0.04166667, 0.04166667 (x, y) extent : -124.793, -67.043, 25.04186, 49.41686 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 data source : in memory names : layer values : 8.84197, 50.52144 (min, max)

建議都非常感謝！

謝謝！

來源

2017-09-25 nate-m

由於您正在計算每月平均值的第90百分位數，因此您將獲得單個圖層。我不太確定你爲什麼計算平均值。我認爲你需要在stackApply函數中運行分位數函數，或者在一個循環中運行分位數函數，在這個循環中你從堆棧中爲每個月的每日值進行分類。當我嘗試使用stackApply方法時出現錯誤，但我會嘗試更多的想法並提交答案。 –

理想情況下，這將是每天的手段，但由於我無法創建第95百分位的每月圖像，我想我會先找出解決方案。我試着在stackApply裏面的每日圖像上運行這個函數，但是一直返回一個錯誤。我真的被卡住了，任何建議都會非常有幫助。謝謝！ –

一種選擇是使用for循環：

x <- stack() # create an empty stack 
for (i in 1:nlayers(mean_tmp)){ 

mean_90thtmp <- calc(mean_tmp[[i]], forceapply = TRUE, 
       fun = function(x) {quantile(x, probs = 0.90, na.rm = TRUE) }) 

x <- stack(x , mean_90thtmp) 
}

來源

2017-09-25 22:18:07

謝謝！我假設/認爲這是我在'calc'函數內部缺少的東西，因爲我一直在'calc'文檔中讀到能夠處理在所有圖層之間應用函數的文檔。特別是當'forceapply = TRUE'時。不過謝謝 - 這是一些值得思考的食物。 –

我不能得到stackApply功能與使用quantile作爲函數工作。

下面是一個使用循環從堆棧中爲每個月選擇所有圖層的方法。

library(raster) 
rstr <- raster('tmmx_1982.nc') 
date_seq <- date_seq[1:nlayers(rstr)] 
month_seq <- month(date_seq) 

outSt <- stack() 
for (mn in 1:12){ 
    st <- subset(rstr, which(month_seq == mn)) 
    mn_90th <- calc(st, fun=function(x) raster::quantile(x, probs=0.9, na.rm=T)) 
    outSt <- addLayer(outSt, mn_90th) 
}

來源

2017-09-26 15:43:28

奇怪的是，您不能在柵格堆棧中的每個圖層上使用stackApply或calc native本地分位函數。感謝這個解決方案，並驗證stackApply問題不只是我！最受讚賞。 –

當應用於柵格堆棧時，柵格計算僅返回一個圖層

回答

相關問題