gIntersection은 R의 gIsValid가 TopologyException을 반환합니다.

Question

같은 질문을했습니다. elsewhere. 이제 데이터를 가져와 오류를 재현 해 볼 수 있습니다.

작은 면적의 대부분이 1997 년까지 방목 된 것을보고 싶습니다. 작은 영역은 setoresp.zip에 있습니다. 삼림 벌채는 d97.zip에 있습니다. 파일은 here입니다.

# 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 

Answer 1

PRODES 폴리곤이 MESS이므로, 래스터 파일을 사용하여 작업을 마쳤습니다. 여기 내 솔루션입니다 : 당신은 개별 개체 사이의 교차점을 할 경우이 작동 이해 할수

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