R如何在多个光栅图层中同时导出我的结果



作为第一步,我导入一个表,其中包含几个采样站的每日数据(请参阅下面的示例(。随后,我将要为每个日期(天(创建光栅图层,并在以后使用不同的名称导出这些光栅。由于我有一个包含每日数据的长时间序列(40年(,我希望制作一个循环(见下面的程序(,但我不能给输出光栅起不同的名字。所以每次他都会在前一层压垮我。我是循环编程的新手,你能帮我解决这个问题吗?

提前感谢

data=read.table(file="data.csv", sep=";", header=TRUE)
head(data)
LAMBX   LAMBY   DATE    S
600 24010   20180801    15.3
600 24010   20180802    12
600 24010   20180803    15
600 24010   20180804    14.8
600 24010   20180805    16.8
600 24010   20180806    15.1
601 24011   20180801    11
601 24011   20180802    14
601 24011   20180803    16.8
601 24011   20180804    15.1
601 24011   20180805    13.8
601 24011   20180806    12.7
602 24012   20180801    15.3
602 24012   20180802    14
602 24012   20180803    12
602 24012   20180804    16.8
602 24012   20180805    17.5
602 24012   20180806    15.1

library(dplyr)
library(sp)
library(raster)
# for each date I will do the following manipulations
for (i in data[,"DATE"]) {
# to delimit my study area
table=filter(data, DATE == i & LAMBX %in% 601:602 & LAMBY %in% 24011:24012) 
# convert geographic coordinates
table[,c("LAMBX", "LAMBY")]= 100*table[,c("LAMBX", "LAMBY")]
# spatialize the stations
xy <- table[,c("LAMBX", "LAMBY")]
sptable <- SpatialPointsDataFrame(coords = xy, data = table,
proj4string = CRS("+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs"))
# rasterize my SpatialPointsDataFrame
rsptable <- rasterFromXYZ(as.data.frame(sptable)[, c("LAMBX", "LAMBY","S")],  crs="+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs")
# export my rasters layers
writeRaster(rsptable, filename="S_date.tif", format="GTiff", overwrite=TRUE)
}

我在您的数据中没有看到T_Q变量,所以我将使用S。他们的技巧是创建一个空列表对象,在每次循环迭代时填充。然后,您可以创建结果的堆栈。堆栈的每一层都是一个单独的日期。

library(dplyr)
library(sp)
library(raster)
data <- read.table(text = "
LAMBX   LAMBY   DATE    S
600 24010   20180801    15.3
600 24010   20180802    12
600 24010   20180803    15
600 24010   20180804    14.8
600 24010   20180805    16.8
600 24010   20180806    15.1
601 24011   20180801    11
601 24011   20180802    14
601 24011   20180803    16.8
601 24011   20180804    15.1
601 24011   20180805    13.8
601 24011   20180806    12.7
602 24012   20180801    15.3
602 24012   20180802    14
602 24012   20180803    12
602 24012   20180804    16.8
602 24012   20180805    17.5
602 24012   20180806    15.1", header = TRUE)
crs_string <- "+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs"
res <- list()
for (i in seq_along(unique(data[,"DATE"]))){
table <- filter(data,
DATE == unique(data[,"DATE"])[i] &
LAMBX %in% 601:602 &
LAMBY %in% 24011:24012)
table[,c("LAMBX", "LAMBY")] <-  100 * table[,c("LAMBX", "LAMBY")]
xy       <- table[,c("LAMBX", "LAMBY")]
sptable  <- SpatialPointsDataFrame(coords = xy, data = table,
proj4string = CRS(crs_string))
rsptable <- rasterFromXYZ(
as.data.frame(sptable)[, c("LAMBX", "LAMBY","S")],  crs = crs_string)
res[[i]] <- rsptable
}
res <- stack(res)
writeRaster(res, filename="S_date.tif", format="GTiff", overwrite=TRUE)

最新更新