0

我在其他地方也提出了同样的问题。现在我带来数据，以便您尝试重现错误。

我只是想看看到 1997 年有多少小区域被砍伐。小区域在 setoresp.zip 中。砍伐森林在 d97.zip 中。文件在这里。将它们与代码一起解压缩到工作目录：

# This script reads sectors of investigated cities and calculates the
#  deforestated area (up to 1997) for the first sector

library(rgdal)
library(rgeos)
library(raster)

setp <- readOGR(".","setoresp")
# to make it faster, I'll try to intersect with the bounding box,
# instead of the actual polygon
rect <- extent(setp[1,]@bbox)
rect <- as(rect, 'SpatialPolygons')

# Deforestation of brazilian Amazon, clipped to the interest area
#  and projected to UTM (so I'll get the intersection area in meters)
d97 <- readOGR(".","d97")
rect@proj4string <- d97@proj4string

if (gIntersects(rect, d97)) {
    print("Intersects!")
    flush.console() # so I'll know the error is below, not above
    rect97 <- gIntersection(rect, d97)
}
#Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : 
#  TopologyException: found non-noded intersection between LINESTRING (533036 -314770, 533036 -314770) and LINESTRING (533036 -314780, 533036 -314770) at 533035.88110651996 -314769.97350772272
4

1 回答 1

0

由于 PRODES 多边形是一团糟，我最终使用了他们的光栅文件。这是我的解决方案：

# This script reads sectors of investigated cities and calculates the
#  deforestated area (in different periods) for each sector    

library(rgdal)
library(rgeos)
library(raster)

desmat <- raster('Prodes2012_AM_90m.tif')

set <- readOGR(".","setoresp",encoding="UTF-8",stringsAsFactors=F)
set@proj4string

#classes
d97 <- 17 #cinza
d06 <- c(2,3,5,6,8,9,11:13,15,16,18,19,24,29,36:39,41,42,44,47,50,54,55,60,62,66,68:71,74,79,81) #rosa
d11 <- c(4,7,14,20:23,25:28,30,33:35,40,43,45,46,48,49,51:53,56:59,61,63,65,72,73,75:78,82,84) #vermelho
d12 <- 80 #roxo
dagua <- 67 #azul
dnuv <- c(32,64) #azul claro (nuvem e resíduo)
dflor <- 31 #verde escuro
dnfl <- 10 #laranja
detc <- c(1,83) #marrom (DV? e DSF_ANT?)

par(mar=c(0,0,0,0))

d <- data.frame(codsetor=character(),
    codmun=character(),
    nomemun=character(),
    areaset=numeric(),
    a97=numeric(), #areas
    a06=numeric(),
    a11=numeric(),
    a12=numeric(),
    afl=numeric(),
    anfl=numeric(),
    anuv=numeric(),
    aagua=numeric(),
    p97=numeric(), #percentages
    p06=numeric(),
    p11=numeric(),
    p12=numeric(),
    pfl=numeric(),
    pnfl=numeric(),
    pnuv=numeric(),
    pagua=numeric(),
    stringsAsFactors=FALSE)

for (i in 1:length(set)) {

    r <- raster(ext=extent(set[i,]))
    setR <- rasterize(set[i,],r)
    desmatP <- projectRaster(from=desmat,to=setR,method='ngb')
    setI <- mask(desmatP,set[i,])
    plot(setI)
    freqI <- freq(setI,useNA='no')
    areas <- freqI[,2] * prod(res(setI))
    freqI <- cbind(freqI,areas)

    print(set@data[i,c(2,19,20)],row.names=F)
    flush.console()

    gArea(set[i,])
    d[nrow(d)+1,] <- c(as.character(set@data[i,2]), #codsetor
        as.character(set@data[i,19]), #codmun
        as.character(set@data[i,20]), #nomemun
        gArea(set[i,])/1e6, #areaset (all areas in Km2)
        sum(freqI[freqI[,1] %in% d97,3])/1e6, #a97
        sum(freqI[freqI[,1] %in% d06,3])/1e6, #a06
        sum(freqI[freqI[,1] %in% d11,3])/1e6, #a11
        sum(freqI[freqI[,1] %in% d12,3])/1e6, #a12
        sum(freqI[freqI[,1] %in% dflor,3])/1e6, #afl
        sum(freqI[freqI[,1] %in% dnfl,3])/1e6, #anfl
        sum(freqI[freqI[,1] %in% dnuv,3])/1e6, #anuv
        sum(freqI[freqI[,1] %in% dagua,3])/1e6, #aagua
        sum(freqI[freqI[,1] %in% d97,3])/gArea(set[i,]), #p97
        sum(freqI[freqI[,1] %in% d06,3])/gArea(set[i,]), #p06
        sum(freqI[freqI[,1] %in% d11,3])/gArea(set[i,]), #p11
        sum(freqI[freqI[,1] %in% d12,3])/gArea(set[i,]), #p12
        sum(freqI[freqI[,1] %in% dflor,3])/gArea(set[i,]), #pfl
        sum(freqI[freqI[,1] %in% dnfl,3])/gArea(set[i,]), #pnfl
        sum(freqI[freqI[,1] %in% dnuv,3])/gArea(set[i,]), #pnuv
        sum(freqI[freqI[,1] %in% dagua,3])/gArea(set[i,])) #pagua
}
write.table(d,"mapsdesm/tabela.txt",quote=FALSE,sep="\t",row.names=FALSE,col.names=T)
于 2013-08-08T16:50:58.583 回答