2

I have two SpatialPolygonsDataFrame, one of them is the deforested area in the Amazon and the other is the city of Vila Rica, which is in the same region. I would like to calculate the total deforested area in Vila Rica, so I'm trying to delimit the deforested area inside the city.

I used the gIntersection function in rgeos package to compare these two SpatialPolygonsDataFrame, generating the following result.

enter image description here

Green is the deforested area in the Amazon and blue is the area defined by gIntersection function.

Apparently the function worked well with smaller polygons, identifying only areas inside the city, but for some reason didn't work for the largest polygon, which crosses the city. I don't know why the function didn't identify the deforested area inside the city for the largest polygon, perhaps a problem in the function?

I've tried the Over method in SP package but it doesn't seem to work like gIntersection.

The shapefiles (Rdata extension) are available in https://www.dropbox.com/sh/1jgn3b3k8hviw1i/AABGBS7zAzoI2YiYPsMisbApa?dl=0.

The "cities" data refers to all the cities in the Amazon, not only for Vila Rica. In reality my main goal is calculate the total deforested area in each city of the region. I used only Vila Rica in the example to simplify.

Here is the R code for the map:

library(rgdal);library(sp)
library(rgeos)

### Download Rdatas and save in workdirectory ########
setwd(<select wd>)

load("deforested.Rdata")
load("cities.Rdata")

proj4string(deforested) <- proj4string(cities)

### Subset Vila Rica from cities shape ########
VilaRica <- subset (cities, nome=='Vila Rica')

### Delimit the deforested areas which are included in the city ###
deforestedVR  <- over(deforested,VilaRica)
deforestedVR  <- deforested[!is.na(deforestedVR),]

### Selecting the intersections ###
intersection <- gIntersection(deforestedVR, VilaRica, byid=TRUE,drop_lower_td=TRUE)

### Plotting the map ###
plot(deforestedVR,col='green')
plot(VilaRica,lwd=2,add=TRUE)
plot(intersection,col='blue',add=TRUE)

I was wondering if the polygon class may be affecting the function, since for the simplest examples it works perfectly.

I separated the largest polygon and changed its class from SpatialPolygonsDataFrame to SpatialPolygons and did the same thing to the Vila Rica shape. Using again the gIntersection function, now works perfectly as you can see in the map:

enter image description here

Here is the R code for the second map:

### Identify and subset the biggest polygon ###
nPolys <- sapply(deforestedVR@polygons, function(x)length(x@Polygons))
bigPolygon <- deforestedVR[which(nPolys==max(nPolys)),]

### Changing class ###
VilaRica <- SpatialPolygons(VilaRica@polygons[1])
bigPolygon <- SpatialPolygons(bigPolygon@polygons[1])

### Selecting the intersections ###
intersection2 <- gIntersection(bigPolygon,VilaRica)

### Plotting the map ###
plot(bigPolygon,col="green")
plot(VilaRica,lwd=2,add=TRUE)
plot(intersection2,col='blue',add=TRUE)

Why for the same polygon the function didn't work for SpatialPolygonsDataFrame class but worked for SpatialPolygons class? Hoping you all can help. Accept other suggestions.

2

It seems like your data has some topological issues (like overlapping polygons, multi or micro polygons inside the largest one). I tried to do it with the original data in QGIS and it had the same result, so I ran a v.clean on a data subset (because original datasets are too large for testing purposes) and using this topology cleaned layer the intersection worked as it should.

  • To begin with, I just selected the polygons of deforested layer around 'Vila Rica' to get a subset of the dataset, I used the Select Features by area tool and I saved it as a new layer.

  • Then I used the v.clean Geoprocessing tool. Processing Toolbox -> Grass commands -> Vector -> v.clean, I used the rmdupl cleaning tool with a threshold of 0.0001.

  • Finally I used this topology cleaned output dataset to perform the intersection with cities layer in menu Vector->Geoprocessing->Intersection

Using the corrected dataset in R with your script, it also worked (Being deforested an object of class SpatialPolygonsDataFrame), my result in R:

enter image description here

|improve this answer|||||

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.