Могу ли я извлекать частоты растровых пикселей, полигон за полигоном, по одному и сохранять, чтобы не перегружать мою оперативную память в R?

Question

Мне нужно извлечь частоты пикселей из растра с помощью SapatialPolygonsDataFrame, но мой растр представляет собой большой объем данных, и мой персональный компьютер не смог его вычислить.

Итак, если есть способ оговорить код, в котором каждый многоугольник моего SapatialPolygonsDataFrame будет рассчитываться отдельно и сохраняться по идентификатору или имени в 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)

Estructure one of my 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

Мой текущий код:

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)

}
}
}

Answer 1

Основано на примере Тима Ассала

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,]))
}

В этом случае вы можете создать таблицу с помощью

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")
}

Answer 2

Думаю, вы слишком сильно над этим задумываетесь (если я правильно понимаю, что вы пытаетесь сделать). По сути, вы можете извлечь растровые данные для каждого многоугольника и суммировать их в фрейм данных. Вот воспроизводимый пример:

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"