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 ポリゴンは 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)
于 2013-08-08T16:50:58.583 に答える