我可以一次提取一个多边形的光栅像素频率并保存，以避免 R 中的 RAM 过载吗？
Can I extract raster pixel frequencies, polygon by polygon, one at a time and save, to avoid overloading my RAM in R?
我需要通过SapatialPolygonsDataFrame从栅格中提取像素频率，但我的栅格数据量很大，我的个人电脑无法计算。
所以，如果有什么办法在代码中规定我的SapatialPolygonsDataFrame的每个多边形都单独计算，并通过ID或名称保存在一个DataFrame中，一个一个，我觉得会有用。
但是，我没有这样做，因为我不知道该怎么做。
我认为另一种可能的解决方案是在新的 SapatialPolygonDataFrame 中分离每个多边形。
但我认为这将是一个问题，因为我将有很多 SapatialPolygonDataFrame 并且重命名它们中的每一个都可能是一个新问题。
构造我的一个栅格（地图）：
> map
class : RasterLayer
dimensions : 86662, 111765, 9685778430 (nrow, ncol, ncell)
resolution : 0.0002689995, 0.0002690002 (x, y)
extent : -73.97832, -43.91358, -18.04061, 5.271491 (xmin, xmax, ymin, ymax)
crs : +proj=aea +lat_1=10 +lat_2=-40 +lat_0=-25 +lon_0=-50 +x_0=0 +y_0=0
+datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
source : C:/Users/kleds/OneDrive/Documentos/mestrado/PRODES_APs/pAP_PI.tif
names : pAP_PI
values : 0, 1 (min, max)
构建我的 SpatialPolygonsDataFrame (SPDF) 之一：
> summary(SPDF)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x -2425864 597166.3
y 1063407 3607311.1
Is projected: TRUE
proj4string :
[+proj=aea +lat_1=10 +lat_2=-40 +lat_0=-25 +lon_0=-50 +x_0=0 +y_0=0
+datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
ADM SIGLA ANO HECTARES
estadual :34 ESEC :18 1990 : 6 Min. : 0
federal : 2 PARNA:19 1995 : 3 1st Qu.: 77095
Federal :39 PE :23 1997 : 3 Median : 219406
municipal: 1 PM : 1 2005 : 3 Mean : 549809
REBIO:12 1998 : 2 3rd Qu.: 672743
REVIS: 1 (Other):20 Max. :4203563
NA's : 2 NA's :39
下面的每个名称都属于一个多边形，我可以为每个多边形创建 76 个新的 SPDF 吗？
>SPDF$NAME_AP
[1] PARQUE ESTADUAL SERRA DO ARAC?
[2] PARQUE ESTADUAL TUCUM?
[3] PARQUE ESTADUAL SUMA?MA
[4] PARQUE ESTADUAL DE MONTE ALEGRE
[5] ESTA??O ECOL?GICA RIO FLOR DO PRADO
[6] PARQUE ESTADUAL DO BACANGA
[7] PARQUE ESTADUAL IGARAP?S DO JURUE
[8] ESTA??O ECOL?GICA DO S?TIO RANGEDOR
[9] PARQUE ESTADUAL DE GUAJAR?-MIRIM
[10] RESERVA BIOL?GICA RIO OURO PRETO
[11] ESTA??O ECOL?GICA DO GR?O PAR?
[12] PARQUE ESTADUAL GUARIBA
[13] ESTA??O ECOL?GICA SAMUEL
[14] ESTA??O ECOL?GICA SERRA DOS TR?S IRM?OS
[15] PARQUE ESTADUAL RIO NEGRO SETOR SUL
[16] PARQUE ESTADUAL DO XINGU
[17] PARQUE ESTADUAL SERRA DOS REIS
[18] PARQUE ESTADUAL DO MATUPIRI
[19] PARQUE ESTADUAL RIO NEGRO SETOR NORTE
[20] REF?GIO DE VIDA SILVESTRE METR?POLE DA AMAZ?NIA
[21] PARQUE ESTADUAL SUCUNDURI
[22] PARQUE ESTADUAL DO UTINGA
[23] PARQUE ESTADUAL DE CORUMBIARA
[24] PARQUE ESTADUAL SERRA SANTA B?RBARA
[25] ESTA??O ECOL?GICA DO RIO ROOSEVELT
[26] PARQUE ESTADUAL CHANDLESS
[27] PARQUE ESTADUAL DO CANT?O
[28] RESERVA BIOL?GICA DE MAICURU
[29] ESTA??O ECOL?GICA DO RIO RONURO
[30] RESERVA BIOL?GICA TRA?ADAL
[31] PARQUE ESTADUAL DA SERRA DOS MART?RIOS/ANDORINHAS
[32] PARQUE ESTADUAL CRISTALINO
[33] PARQUE ESTADUAL CHARAPUCU
[34] PARQUE ESTADUAL SERRA RICARDO FRANCO
[35] PARQUE TURAL MUNICIPAL DO CANC?O
[36] Parna do Rio Novo
[37] Esec do Jari
[38] Rebio do Uatum?
[39] Parna do Viru?
[40] Esec de Marac?
[41] Parna do Pico da Neblina
[42] Parna do Ja?
[43] Parna do Monte Roraima
[44] Esec de Marac?-Jipioca
[45] Rebio do Lago Piratuba
[46] Rebio do Tapirap?
[47] Rebio do Abufari
[48] Parna de Paca?s Novos
[49] Esec Rio Acre
[50] Rebio do Guapor?
[51] Parna Serra da Cutia
[52] Esec da Terra do Meio
[53] Parna da Serra do Divisor
[54] Esec Juami-Japur?
[55] Esec de Juta?-Solim?es
[56] Esec de Caracara?
[57] Esec Niqui?
[58] Parna do Cabo Orange
[59] Rebio do Rio Trombetas
[60] Parna da Serra do Pardo
[61] Rebio Nascentes da Serra do Cachimbo
[62] Parna do Jamanxim
[63] Rebio do Jaru
[64] Parna Nascentes do Lago Jari
[65] Parna Montanhas do Tumucumaque
[66] Parna dos Campos Amaz?nicos
[67] Parna Mapinguari
[68] Parna da Amaz?nia
[69] Esec de Cuni?
[70] Parna da Serra da Mocidade
[71] Parna do Juruena
[72] Rebio do Gurupi
[73] Parna de Anavilhanas
[74] Esec Alto Mau?s
[75] RESERVA BIOLiGICA DO MANICORn
[76] PARQUE CIOL DO ACARI
76 Levels: Esec Alto Mau?s Esec da Terra do Meio ... RESERVA BIOLiGICA DO MANICORn
My current code is:
dirs<-"~/prodes/PRODES_APs"
work_dirs<-"~/prodes/PRODES_APs"
#Create a for to define the rasters directory, and to be used in the subsequent for
for (m in 1:length(dirs)) {
files<-file.path(dirs[m],list.files(path = dirs[m], pattern = ".tif"))
nomes <- list.files(path = dirs[m], pattern = ".tif")
nomes <- substr(nomes,1,nchar(nomes)-4)
}
#create a for to call simultaneously raster layer of interest, and each SPDF (initial polygons, rings and control)
#vectors to use in the for
AP<-c("PI","TI","UN","US")
AW <- c("arc","wood")
km<-c("min","1km_","2km_","3km_","4km_","5km_","6km_","7km_","8km_","9km_",10km_,"10km","20km","30km","40km","50km","60km","70km", "controle")
#empty Data Frame to save my results
results<-data.frame()
for (a in 1:(min(length(files), length(AP)))){
setwd(work_dirs)
r<-files[a]
i<- AP[a]
map<- raster(r)
for(k in AW){
for(j in km){
# deffine the directory
setwd(paste0("~/prodes/buff_",k,"/AP_rings"))
# Call each SPDF
SPDF<- readOGR(".", paste0("ring",k,j, i))
# reproject the SPDF to ALbers
SPDF <- spTransform(SPDF , CRSobj = "+proj=longlat +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +no_defs ")
#Extract the pixels values
( extrc <- extract(map, SPDF, na.rm=T) )
#proportion calculation for each class
(class.prop = lapply(extrc, function(x)
{prop.table(table(factor(x,levels=c(0,1))))}))
p.prop = setNames(
do.call(
rbind.data.frame,
class.prop),
c("Desmatado","natural"))
p.prop$ID<-seq_along(p.prop[,1])
rings$ID<- 1:length(SPDF)
freq <- merge(SPDF, p.prop) #add to polygons
frequenc<-as.data.frame(freq)
View(frequenc)
results <- rbind(results, frequenc)
setwd("~/prodes/resultados")
write.table(results, file="resultados.txt", sep="\t", row.names=F)
}
}
}
我认为您可能想多了（如果我理解您正在尝试正确地做的事情的话）。本质上，您可以为每个多边形提取栅格数据并将其汇总到一个数据框中。这是一个可重现的例子：
library(raster)
# create example raster and polygon
p <- raster(nrow=10, ncol=10)
p[] <- runif(ncell(p)) * 10
p <- rasterToPolygons(p, fun=function(x){x > 9})
r <- raster(nrow=100, ncol=100)
r[] <- runif(ncell(r))
plot(r)
plot(p, add=TRUE, lwd=4)
看起来您的所有栅格数据的值为 1。但是，如果您需要设置阈值（此处基于 0.5），则需要重新分类：
#reclassify if needed
r2<-reclassify(r, c(-Inf, 0.5, NA, 0.5, Inf, 1))
plot(r2)
plot(p, add=TRUE, lwd=4)
#extract values from each polygon and write to list
( e <- extract(r2, p) )
#summarize raster data in each polygon; again, here all values have been reclassified to 1
out.list<-( class.counts <- lapply(e, table) )
#convert to dataframe and rename column
out.df <- data.frame(matrix(unlist(out.list), nrow=length(out.list), byrow=T))
colnames(out.df) <- "Freq"
以 Tim Assal 的例子为基础
library(raster)
set.seed(91)
p <- raster(nrow=10, ncol=10)
values(p) <- runif(ncell(p)) * 10
p <- rasterToPolygons(p, fun=function(x){x > 9})
r <- raster(nrow=100, ncol=100)
values(r) <- sample(5, ncell(r), replace=TRUE)
plot(r)
plot(p, add=TRUE, lwd=4)
最简单的方法是添加减少观察数量的函数 extract
e <- extract(r, p, fun=function(i,...) table(i) )
e
# 1 2 3 4 5
# [1,] 16 22 22 20 20
# [2,] 28 22 17 21 12
# [3,] 24 15 13 24 24
# [4,] 27 20 13 18 22
# [5,] 19 20 25 21 15
# [6,] 13 18 20 22 27
# [7,] 15 28 15 22 20
# [8,] 18 22 23 25 12
# [9,] 17 22 22 18 21
#[10,] 18 20 16 23 23
如果所有情况都出现在所有多边形中，您将得到一个矩阵，否则您将得到一个列表，如下所示。
以上相当于
out <- list()
for (i in 1:length(p)) {
out[[i]] <- table(extract(r, p[i,]))
}
在这种情况下，您可以用
从外面制作一个 table
x <- do.call(rbind, out)
另一种方法，对于非常大的多边形也应该是内存安全的，可以是
out <- list()
for (i in 1:length(p)) {
x <- crop(r, p[i,])
x <- mask(r, p[i,])
out[[i]] <- freq(x, useNA="no")
}
我需要通过SapatialPolygonsDataFrame从栅格中提取像素频率，但我的栅格数据量很大，我的个人电脑无法计算。
所以，如果有什么办法在代码中规定我的SapatialPolygonsDataFrame的每个多边形都单独计算，并通过ID或名称保存在一个DataFrame中，一个一个，我觉得会有用。 但是，我没有这样做，因为我不知道该怎么做。
我认为另一种可能的解决方案是在新的 SapatialPolygonDataFrame 中分离每个多边形。 但我认为这将是一个问题，因为我将有很多 SapatialPolygonDataFrame 并且重命名它们中的每一个都可能是一个新问题。
构造我的一个栅格（地图）：
> map
class : RasterLayer
dimensions : 86662, 111765, 9685778430 (nrow, ncol, ncell)
resolution : 0.0002689995, 0.0002690002 (x, y)
extent : -73.97832, -43.91358, -18.04061, 5.271491 (xmin, xmax, ymin, ymax)
crs : +proj=aea +lat_1=10 +lat_2=-40 +lat_0=-25 +lon_0=-50 +x_0=0 +y_0=0
+datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
source : C:/Users/kleds/OneDrive/Documentos/mestrado/PRODES_APs/pAP_PI.tif
names : pAP_PI
values : 0, 1 (min, max)
构建我的 SpatialPolygonsDataFrame (SPDF) 之一：
> summary(SPDF)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x -2425864 597166.3
y 1063407 3607311.1
Is projected: TRUE
proj4string :
[+proj=aea +lat_1=10 +lat_2=-40 +lat_0=-25 +lon_0=-50 +x_0=0 +y_0=0
+datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
ADM SIGLA ANO HECTARES
estadual :34 ESEC :18 1990 : 6 Min. : 0
federal : 2 PARNA:19 1995 : 3 1st Qu.: 77095
Federal :39 PE :23 1997 : 3 Median : 219406
municipal: 1 PM : 1 2005 : 3 Mean : 549809
REBIO:12 1998 : 2 3rd Qu.: 672743
REVIS: 1 (Other):20 Max. :4203563
NA's : 2 NA's :39
下面的每个名称都属于一个多边形，我可以为每个多边形创建 76 个新的 SPDF 吗？
>SPDF$NAME_AP
[1] PARQUE ESTADUAL SERRA DO ARAC?
[2] PARQUE ESTADUAL TUCUM?
[3] PARQUE ESTADUAL SUMA?MA
[4] PARQUE ESTADUAL DE MONTE ALEGRE
[5] ESTA??O ECOL?GICA RIO FLOR DO PRADO
[6] PARQUE ESTADUAL DO BACANGA
[7] PARQUE ESTADUAL IGARAP?S DO JURUE
[8] ESTA??O ECOL?GICA DO S?TIO RANGEDOR
[9] PARQUE ESTADUAL DE GUAJAR?-MIRIM
[10] RESERVA BIOL?GICA RIO OURO PRETO
[11] ESTA??O ECOL?GICA DO GR?O PAR?
[12] PARQUE ESTADUAL GUARIBA
[13] ESTA??O ECOL?GICA SAMUEL
[14] ESTA??O ECOL?GICA SERRA DOS TR?S IRM?OS
[15] PARQUE ESTADUAL RIO NEGRO SETOR SUL
[16] PARQUE ESTADUAL DO XINGU
[17] PARQUE ESTADUAL SERRA DOS REIS
[18] PARQUE ESTADUAL DO MATUPIRI
[19] PARQUE ESTADUAL RIO NEGRO SETOR NORTE
[20] REF?GIO DE VIDA SILVESTRE METR?POLE DA AMAZ?NIA
[21] PARQUE ESTADUAL SUCUNDURI
[22] PARQUE ESTADUAL DO UTINGA
[23] PARQUE ESTADUAL DE CORUMBIARA
[24] PARQUE ESTADUAL SERRA SANTA B?RBARA
[25] ESTA??O ECOL?GICA DO RIO ROOSEVELT
[26] PARQUE ESTADUAL CHANDLESS
[27] PARQUE ESTADUAL DO CANT?O
[28] RESERVA BIOL?GICA DE MAICURU
[29] ESTA??O ECOL?GICA DO RIO RONURO
[30] RESERVA BIOL?GICA TRA?ADAL
[31] PARQUE ESTADUAL DA SERRA DOS MART?RIOS/ANDORINHAS
[32] PARQUE ESTADUAL CRISTALINO
[33] PARQUE ESTADUAL CHARAPUCU
[34] PARQUE ESTADUAL SERRA RICARDO FRANCO
[35] PARQUE TURAL MUNICIPAL DO CANC?O
[36] Parna do Rio Novo
[37] Esec do Jari
[38] Rebio do Uatum?
[39] Parna do Viru?
[40] Esec de Marac?
[41] Parna do Pico da Neblina
[42] Parna do Ja?
[43] Parna do Monte Roraima
[44] Esec de Marac?-Jipioca
[45] Rebio do Lago Piratuba
[46] Rebio do Tapirap?
[47] Rebio do Abufari
[48] Parna de Paca?s Novos
[49] Esec Rio Acre
[50] Rebio do Guapor?
[51] Parna Serra da Cutia
[52] Esec da Terra do Meio
[53] Parna da Serra do Divisor
[54] Esec Juami-Japur?
[55] Esec de Juta?-Solim?es
[56] Esec de Caracara?
[57] Esec Niqui?
[58] Parna do Cabo Orange
[59] Rebio do Rio Trombetas
[60] Parna da Serra do Pardo
[61] Rebio Nascentes da Serra do Cachimbo
[62] Parna do Jamanxim
[63] Rebio do Jaru
[64] Parna Nascentes do Lago Jari
[65] Parna Montanhas do Tumucumaque
[66] Parna dos Campos Amaz?nicos
[67] Parna Mapinguari
[68] Parna da Amaz?nia
[69] Esec de Cuni?
[70] Parna da Serra da Mocidade
[71] Parna do Juruena
[72] Rebio do Gurupi
[73] Parna de Anavilhanas
[74] Esec Alto Mau?s
[75] RESERVA BIOLiGICA DO MANICORn
[76] PARQUE CIOL DO ACARI
76 Levels: Esec Alto Mau?s Esec da Terra do Meio ... RESERVA BIOLiGICA DO MANICORn
My current code is:
dirs<-"~/prodes/PRODES_APs"
work_dirs<-"~/prodes/PRODES_APs"
#Create a for to define the rasters directory, and to be used in the subsequent for
for (m in 1:length(dirs)) {
files<-file.path(dirs[m],list.files(path = dirs[m], pattern = ".tif"))
nomes <- list.files(path = dirs[m], pattern = ".tif")
nomes <- substr(nomes,1,nchar(nomes)-4)
}
#create a for to call simultaneously raster layer of interest, and each SPDF (initial polygons, rings and control)
#vectors to use in the for
AP<-c("PI","TI","UN","US")
AW <- c("arc","wood")
km<-c("min","1km_","2km_","3km_","4km_","5km_","6km_","7km_","8km_","9km_",10km_,"10km","20km","30km","40km","50km","60km","70km", "controle")
#empty Data Frame to save my results
results<-data.frame()
for (a in 1:(min(length(files), length(AP)))){
setwd(work_dirs)
r<-files[a]
i<- AP[a]
map<- raster(r)
for(k in AW){
for(j in km){
# deffine the directory
setwd(paste0("~/prodes/buff_",k,"/AP_rings"))
# Call each SPDF
SPDF<- readOGR(".", paste0("ring",k,j, i))
# reproject the SPDF to ALbers
SPDF <- spTransform(SPDF , CRSobj = "+proj=longlat +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +no_defs ")
#Extract the pixels values
( extrc <- extract(map, SPDF, na.rm=T) )
#proportion calculation for each class
(class.prop = lapply(extrc, function(x)
{prop.table(table(factor(x,levels=c(0,1))))}))
p.prop = setNames(
do.call(
rbind.data.frame,
class.prop),
c("Desmatado","natural"))
p.prop$ID<-seq_along(p.prop[,1])
rings$ID<- 1:length(SPDF)
freq <- merge(SPDF, p.prop) #add to polygons
frequenc<-as.data.frame(freq)
View(frequenc)
results <- rbind(results, frequenc)
setwd("~/prodes/resultados")
write.table(results, file="resultados.txt", sep="\t", row.names=F)
}
}
}
我认为您可能想多了（如果我理解您正在尝试正确地做的事情的话）。本质上，您可以为每个多边形提取栅格数据并将其汇总到一个数据框中。这是一个可重现的例子：
library(raster)
# create example raster and polygon
p <- raster(nrow=10, ncol=10)
p[] <- runif(ncell(p)) * 10
p <- rasterToPolygons(p, fun=function(x){x > 9})
r <- raster(nrow=100, ncol=100)
r[] <- runif(ncell(r))
plot(r)
plot(p, add=TRUE, lwd=4)
看起来您的所有栅格数据的值为 1。但是，如果您需要设置阈值（此处基于 0.5），则需要重新分类：
#reclassify if needed
r2<-reclassify(r, c(-Inf, 0.5, NA, 0.5, Inf, 1))
plot(r2)
plot(p, add=TRUE, lwd=4)
#extract values from each polygon and write to list
( e <- extract(r2, p) )
#summarize raster data in each polygon; again, here all values have been reclassified to 1
out.list<-( class.counts <- lapply(e, table) )
#convert to dataframe and rename column
out.df <- data.frame(matrix(unlist(out.list), nrow=length(out.list), byrow=T))
colnames(out.df) <- "Freq"
以 Tim Assal 的例子为基础
library(raster)
set.seed(91)
p <- raster(nrow=10, ncol=10)
values(p) <- runif(ncell(p)) * 10
p <- rasterToPolygons(p, fun=function(x){x > 9})
r <- raster(nrow=100, ncol=100)
values(r) <- sample(5, ncell(r), replace=TRUE)
plot(r)
plot(p, add=TRUE, lwd=4)
最简单的方法是添加减少观察数量的函数 extract
e <- extract(r, p, fun=function(i,...) table(i) )
e
# 1 2 3 4 5
# [1,] 16 22 22 20 20
# [2,] 28 22 17 21 12
# [3,] 24 15 13 24 24
# [4,] 27 20 13 18 22
# [5,] 19 20 25 21 15
# [6,] 13 18 20 22 27
# [7,] 15 28 15 22 20
# [8,] 18 22 23 25 12
# [9,] 17 22 22 18 21
#[10,] 18 20 16 23 23
如果所有情况都出现在所有多边形中，您将得到一个矩阵，否则您将得到一个列表，如下所示。
以上相当于
out <- list()
for (i in 1:length(p)) {
out[[i]] <- table(extract(r, p[i,]))
}
在这种情况下，您可以用
从外面制作一个 tablex <- do.call(rbind, out)
另一种方法，对于非常大的多边形也应该是内存安全的，可以是
out <- list()
for (i in 1:length(p)) {
x <- crop(r, p[i,])
x <- mask(r, p[i,])
out[[i]] <- freq(x, useNA="no")
}