R Как экспортировать мои результаты в несколько растровых слоев одновременно? - PullRequest
0 голосов
/ 28 ноября 2018

В качестве первого шага я импортирую таблицу, которая содержит ежедневные данные для нескольких станций отбора проб (см. Образец ниже).Впоследствии я хочу создать растровые слои для каждой даты (дня) и экспортировать эти растры позже с разными именами.Поскольку у меня есть длинный временной ряд (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)
}

1 Ответ

0 голосов
/ 29 ноября 2018

Я не вижу переменную 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)
...