0

我需要通过 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

我目前的代码是：

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)

}
}
}
4

2 回答 2

1

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

在这种情况下，您可以使用 

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")
}
于 2020-05-07T23:18:22.310 回答
0

我认为您可能想多了（如果我了解您要正确执行的操作）。本质上，您可以提取每个多边形的栅格数据并将其汇总到数据框中。这是一个可重现的示例：

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"
于 2020-05-07T19:29:18.073 回答