R：通過使用形狀文件和柵格中的變量完成的計算創建新的光柵/形狀文件

Question

我目前正在嘗試基於有條件的計算創建一個新的光柵或形狀文件，需要對每個值基於光柵文件中的值的形狀值。我通常不會使用光柵和形狀文件，所以我在這裏完全沒有我的元素。我問這籠統，但這裏是我用的，所以希望它會給一個更好的理解什麼，我試圖完成的數據：

rast_norm <- ftp://prism.nacse.org/normals_4km/tmean/PRISM_tmean_30yr_normal_4kmM2_04_bil.zip 
shp_probs <- ftp://ftp.cpc.ncep.noaa.gov/GIS/us_tempprcpfcst/seastemp_201603.zip 

主要目標，是要與相關的概率shp_probs中的每個點（經度和緯度）並乘以與rast_norm中相同經度和緯度對應的值以及其後的一些其他計算。如果我有兩個data.tables，我可以做類似如下：使用數據表

dt1 <- data.table(col1 = c(0:3), col2 = c(1:4)*11, factor1 = sqrt(c(285:288)) 

# # Output # # 
# col1 col2 factor1 
# 0 11 16.88194 
# 1 22 16.91153 
# 2 33 16.94107 
# 3 44 16.97056 

dt2 <- data.table(col1 = c(0:3), col2 = c(1:4)*11, factor2 = abs(sin(c(1:4)))) 

# # Output # # 
# col1 col2 factor1 
# 0 11 0.8414710 
# 1 22 0.9092974 
# 2 33 0.1411200 
# 3 44 0.7568025 

dt3 <- merge(dt1, dt2, by = c("col1", "col2")) 
dt3$factor1 <- dt3$factor1 * dt3$factor2 
dt3$factor2 <- NULL 

# # Output # # 
# col1 col2 factor1 
# 0 11 14.205665 
# 1 22 15.377615 
# 2 33 2.390725 
# 3 44 12.843364  

易peasy。但是我試圖用Raster和SpatialPolygonsDataFrame來做到這一點。這裏是我迄今爲止在閱讀和清理文件：

# Importing the "rast_norm" file, the first listed above with a link 
rast_norm <- "/my/file/path/PRISM_tmean_30yr_normal_4kmM2_04_bil.zip" 
zipdirec <- "/my/zip/directory" 
unzip(rast_norm, exdir = zipdirec) 

# Get the correct file from the file list 
rast_norm <- list.files(zipdirec, full.names = TRUE, pattern = ".bil") 
rast_norm <- rast_norm[!grepl("\\.xml", rast_norm)] 

# Convert to raster 
rast_norm <- raster(rast_norm) 

Plotting rast_norm on its own gives this map.

# Importing the "shp_probs" file, the second listed above with a link 
shp_probs <- "/my/file/path/seastemp_201603.zip" 
zipdirec <- "/my/zip/directory" 
unzip(shp_probs, exdir = zipdirec, overwrite = TRUE) 

# Get the correct file from the list of file names and find the layer name 
layer_name <- list.files(zipdirec, pattern = "lead14") 
layer_name <- layer_name[grepl(".shp", layer_name)] 
layer_name <- layer_name[!grepl("\\.xml", layer_name)] 
layer_name <- do.call("rbind", strsplit(layer_name, "\\.shp"))[,1] 
layer_name <- unique(layer_name) 

# Use the layer name to read in the shape file 
shp_probs <- readOGR(shp_probs, layer = layer_name) 
names_levels <- paste0(shp_probs$Cat, shp_probs$Prob) 
names_levels <- gsub("Below", "-", names_levels) 
names_levels <- gsub("Above", "+", names_levels) 
names_levels <- as.integer(names_levels) 
shp_probs@data$id <- names_levels 
shp_probs <- as(shp_probs, "SpatialPolygons") 

# Create a data frame of values to use in conjunction with the existing id's 
weights <- data.table(id = shp_probs$id, weight = shp_probs$id) 
weights$weight <- c(.80, .80, .10, .10, .10, .10, .10, .10, .80, .10, .10, .10, .10, .10) 
shp_probs <- SpatialPolygonsDataFrame(otlk_sp, weights, match.ID = FALSE) 

Plotting shp_probs on its own gives this map.

我現在要採取與該shp_probs文件相關聯的概率和將其乘以與rast_norm文件相關的降雨量，再乘以與shp_probs文件中概率相關的權重。

我真的不知道該怎麼做，任何幫助將非常感激。如何提取匹配緯度和經度的所有相應數據點？我想如果我知道，我會知道該怎麼做。

謝謝，提前。

Answer 1

假設要執行此計算爲你的柵格每個網格單元，你可以做這樣的事情：

1. 下載/讀取數據，並添加weight列。請注意，在這裏我使用了隨機權重，因爲您的示例似乎將14個權重分配給7個多邊形。另外，我不確定您的id欄的用途是什麼，所以我跳過了這一部分。

   library(raster) 
   library(rgdal) 
   
   download.file('ftp://prism.nacse.org/normals_4km/tmean/PRISM_tmean_30yr_normal_4kmM2_04_bil.zip', 
          fr <- tempfile(), mode='wb') 
   download.file('ftp://ftp.cpc.ncep.noaa.gov/GIS/us_tempprcpfcst/seastemp_201603.zip', 
          fs <- tempfile(), mode='wb') 
   
   unzip(fr, exdir=tempdir()) 
   unzip(fs, exdir=tempdir()) 
   
   r <- raster(file.path(tempdir(), 'PRISM_tmean_30yr_normal_4kmM2_04_bil.bil')) 
   s <- readOGR(tempdir(), 'lead14_Apr_temp') 
   s$weight <- runif(length(s)) 
   

2. 執行柵格單元和所述多邊形的座標的空間疊置。 （或者，您也可以使用raster::rasterize兩次到Prob和id領域轉換成柵格，然後再乘以三個柵格。）

   xy <- SpatialPoints(coordinates(r), proj4string=crs(r)) 
   o <- over(xy, s) 
   

3. 具有相同程度的創建一個新的光柵/尺寸與原來的柵格和將適當的值分配給它的單元格。

   r2 <- raster(r) 
   r2[] <- r[] * o$Prob * o$weight 
   

有了這些隨機數據，結果看起來是這樣的：