R：根据使用形状文件和栅格中的变量完成的计算创建一个新的 raster/shape 文件
R: Creating a new raster/shape file from calculations done using variables from a shape file and a raster
我目前正在尝试根据条件计算创建新的光栅或形状文件，需要根据光栅文件中的值对形状值中的每个值进行计算。我通常不使用光栅和形状文件，所以我在这里很不适应。我是笼统地问这个问题，但这是我正在使用的数据，希望它能让我更好地理解我想要完成的事情：
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
主要objective是取shp_probs中每个点（经纬度）关联的概率乘以[=35中相同经纬度对应的值=]，以及之后的一些其他计算。如果我有两个 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
使用数据表简单易行。但是我在尝试使用 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 文件中。
我真的不知道该怎么做，非常感谢任何帮助。如何提取所有相应的数据点以匹配纬度和经度？我想如果我知道了，我就会知道该怎么做。
提前谢谢你。
假设您要对栅格的每个网格单元执行此计算，您可以这样做：
Download/read 数据，并添加 weight 列。请注意，这里我只使用了随机权重，因为您的示例似乎为 7 个多边形分配了 14 个权重。另外，我不确定您的 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))
对栅格像元坐标和多边形坐标进行空间叠加。 （或者，您可以使用 raster::rasterize 两次将 Prob 和 id 字段转换为栅格，然后将三个栅格相乘。）
xy <- SpatialPoints(coordinates(r), proj4string=crs(r))
o <- over(xy, s)
创建一个与原始栅格具有相同 extent/dimensions 的新栅格，并为其像元分配适当的值。
r2 <- raster(r)
r2[] <- r[] * o$Prob * o$weight
使用这些随机数据，结果如下所示：
我目前正在尝试根据条件计算创建新的光栅或形状文件，需要根据光栅文件中的值对形状值中的每个值进行计算。我通常不使用光栅和形状文件，所以我在这里很不适应。我是笼统地问这个问题，但这是我正在使用的数据，希望它能让我更好地理解我想要完成的事情：
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
主要objective是取shp_probs中每个点（经纬度）关联的概率乘以[=35中相同经纬度对应的值=]，以及之后的一些其他计算。如果我有两个 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
使用数据表简单易行。但是我在尝试使用 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 文件中。
我真的不知道该怎么做，非常感谢任何帮助。如何提取所有相应的数据点以匹配纬度和经度？我想如果我知道了，我就会知道该怎么做。
提前谢谢你。
假设您要对栅格的每个网格单元执行此计算，您可以这样做：
Download/read 数据，并添加
weight列。请注意，这里我只使用了随机权重，因为您的示例似乎为 7 个多边形分配了 14 个权重。另外，我不确定您的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))对栅格像元坐标和多边形坐标进行空间叠加。 （或者，您可以使用
raster::rasterize两次将Prob和id字段转换为栅格，然后将三个栅格相乘。）xy <- SpatialPoints(coordinates(r), proj4string=crs(r)) o <- over(xy, s)创建一个与原始栅格具有相同 extent/dimensions 的新栅格，并为其像元分配适当的值。
r2 <- raster(r) r2[] <- r[] * o$Prob * o$weight
使用这些随机数据，结果如下所示：