Зацикливание растровых функций пакета над каждым многоугольником в R - PullRequest
0 голосов
/ 23 апреля 2020

Я пытаюсь l oop несколько функций из пакета 'растр', а именно: crop (), mask (), reclassify () и unstack () / as.list (). У меня есть десять растровых слоев, которые имеют одинаковый экстент и тип данных; они соответствуют земному покрову за 10 временных периодов. Я хочу создать индивидуальные переменные списка для каждого выхода функции crop () -> mask () -> reclassify () -> as.list (). Мне удалось передать процесс для одного полигонального объекта, но мне нужно иметь возможность l oop для каждого из 10 полигональных объектов, хранящихся в шейп-файле многоугольника, так что я могу сохранить каждый выходной список в соответствии с заданным именем конвенция.

Спасибо и, пожалуйста, сообщите. Я поделился своим кодом ниже.

РЕДАКТИРОВАТЬ: Мне интересно, является ли for-l oop правильный путь к go по этому поводу, или лучше будет неудачный подход?

# Load libraries
library(raster)             # for raster processing
library(rgdal)              # for raster/vector processing
library(sf)                 # for Shapefile processing

# Stack 10 rasters together
raster.stack = stack( 
  raster("path/raster1.tif"),
  raster("path/raster2.tif"),
  raster("path/raster3.tif"),
  raster("path/raster4.tif"),
  raster("path/raster5.tif"),
  raster("path/raster6.tif"),
  raster("path/raster7.tif"),
  raster("path/raster8.tif"),
  raster("path/raster9.tif"),
  raster("path/raster10.tif")
)

# Prepare reclassification codes from 9-class raster to 3-class raster
reclasscodes = c(
  0,0, # no data 
  1,1,
  2,1,
  3,1,
  4,1,
  5,2,
  6,2,
  7,3,
  8,3,
  9,3
)

# Convert reclass codes list into n x 2 matrix
reclassmatrix = matrix(reclasscodes, ncol=2, byrow = T)

# Load multipolygon vector Shapefile
multipolygon = shapefile("path/multipolygon.shp") # Shapefile is made of n polygons

# Example subset Shapefile to polygon_1 using attribute "ID"
polygon_1 = subset(multipolygon,ID=="D-4")

# Create output for polygon_1
list_polygon_1 =
  raster.stack %>%
  crop(y = polygon_1) %>%              # crop to bounds
  mask(mask = polygon_1) %>%           # mask to polygon cutline
  reclassify(rcl = reclassmatrix) %>%  # reclassify to 3-class
  as.list() # functions the same as unstack() where raster brick is converted to list of raster layers
# I use %>% because I do not want to save any of the intermediate outputs.
# Resultant output is a variable list for polygon_1 named 'list_polygon_1' which is exactly what I want.
# Worked perfectly.

# How do I repeat this process for polygon_1 to polygon n?

# My attempt
for (i in 1:nrow(multipolygon)) {
  raster.stack %>%
  crop(y = multipolygon[i,]) %>%
  mask(mask = multipolygon[i,]) %>%
  reclassify(rcl = reclassmatrix) %>%
  as.list() %>% # up till here it is the same steps as before for polygon_1
  # now I want to save each list output as a separate variable according to i, e.g. list_polygon_2, list_polygon_3 etc.
  assign(paste(multipolygon$ID, i, sep = '_')) # assign a naming convention for each output variable 
}
# Does not work. Even without the last line of code "..assign(paste(...))" there is no output variable from the as.list() line.

1 Ответ

1 голос
/ 24 апреля 2020

Пожалуйста, всегда включайте минимальный автономный воспроизводимый пример

Пример данных

library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster")) 

xy1 <- xy2 <- xy3 <- matrix(c(10,17, 6,10,71,60,62,71), ncol=2)
xy2[,1] <- xy2[,1] + 30
xy3[,2] <- xy3[,2] - 30
p <- spPolygons(xy1, xy2, xy3)
#plot(r, 1)
#lines(p)

Что вы после

rm = matrix(c(0,100,0,100,150,2,150,255,3), ncol=3, byrow=TRUE)

out <- list()
for (i in 1:length(p)) {
    x <- crop(s, p[i,])
    x <- mask(x, p[i,])
    out[[i]] <- reclassify(x, rm)
}

То, что вы говорите о unstack, делает не имеет смысла (и unlist не работает). Я бы посоветовал против этого, но вы могли бы сделать

out2 <- lapply(out, unstack)

Я не уверен, что вы действительно после. Если вам нужны значения ячеек, вы можете сделать это намного проще (не нужно указывать al oop) и сделать

r <- reclassify(s, rm)
e <- extract(r, p)

К вашему вопросу о lapply vs al oop. С точки зрения производительности это редко имеет значение. lapply может быть кратким, но в подобных случаях лучше писать al oop, так как его легче читать и писать, особенно если вы не используете %>%.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...