Search code examples
rextractpolygonrasterspdf

Can I extract raster pixel frequencies, polygon by polygon, one at a time and save, to avoid overloading my RAM in R?

I need to extract the pixel frequencies from raster by SapatialPolygonsDataFrame, but my raster is a large volume of data and my personal computer was unable to calculate it.

So, if there is any way to stipulate in the code that each polygon of my SapatialPolygonsDataFrame will be calculated separately and saved by ID or name in a DataFrame, one by one, I think it will be useful. But, I didn't it because I don't know how can do it.

Another possible solution, which I think, is to separate each polygon, in a new SapatialPolygonDataFrame. But I think that will be a problem, because I will have a lot of SapatialPolygonDataFrame and renaming each one of them can be a new problem.

Estructure one of my rasters (map):

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

enter image description here

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

enter image description here

Each name below belongs to one polygon, could I create 76 new SPDFs one to each polygon?

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

}
}
}
Solution

  • Building on Tim Assal's example

    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)

    The simplest approach is to add the function the reduce the number of observations to 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

    If all case are present in all polygons, you get a matrix, otherwise you get a list, like below.

    The above is equivalent to

    out <- list()
for (i in 1:length(p)) {
    out[[i]] <- table(extract(r, p[i,]))
}

    In this case you can make a table from out with 

    x <- do.call(rbind, out)

    Another approach, that should also be memory safe for very large polygons, could be 

    out <- list()
for (i in 1:length(p)) {
    x <- crop(r, p[i,])
    x <- mask(r, p[i,])
    out[[i]] <- freq(x, useNA="no")
}