2

I have two shapefiles. I would like to identify the part of the first shapefile that doesn't overlap with the second shapefile. Then I would like to clip that part of the first shapefile and make it its own shapefile to merge with the second shapefile in R. I have tried intersect and union in the raster package and gIntersection in the rgeos package but haven't had any luck.

After trying to use the answer below from Simbamangu I think the reason the answer is not working for my shapefiles is that the coordinates for my two shapefiles are formatted differently. The answer works with the example data.

Shapefile attributes

Shapefile 1 output from attributes Shapefile 2 output from attributes

  $bbox                                            $bbox
        min       max                                min       max
x -67.33333 -66.33739                              x  639772.8  723467.9
y  40.58634  41.50000                              y 4499609.0 4597064.1

Is there a way to format the coordinates so both shapefiles are formatted the same way and could this be the reason gIntersection is not working.

intersect1<-rgeos::gIntersection(shapefile1, shpaefile2)
NULL

I figured out to format the coordinates from the 2 shapefiles into a similar format. The gDifference function doesn't seem to be identifying all of the differences between the two shapefiles. See output below.

shapefile pictures

Output from getting answer to work. The only part of the Shapefile 1 identified as outside of Shapefile 2 is the little sliver of green, but the 2 sections in white should also be identified.

enter image description here Here are two example polygons - not the actually shapefiles

# Create SpatialPolygons objects
polygon1 <- readWKT("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20))")
polygon2 <- readWKT("POLYGON((-180 -20, -140 70, 10 0, -140 -60, -180 -20))")

R information

R version 3.3.1 (2016-06-21)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] grid      datasets  utils     stats     graphics  grDevices methods  
[8] base     

other attached packages:
 [1] taRifx.geo_1.0.7    devtools_1.12.0     rgeos_0.3-19       
 [4] raster_2.5-8        calibrate_1.7.2     spatialEco_0.1-5   
 [7] spatstat_1.46-1     rpart_4.1-10        nlme_3.1-128       
[10] rgdal_1.1-10        maptools_0.8-39     mapdata_2.2-6      
[13] maps_3.1.0          sp_1.2-3            car_2.1-2          
[16] xlsx_0.5.7          xlsxjars_0.6.1      rJava_0.9-8        
[19] plyr_1.8.4          MASS_7.3-45         RODBC_1.3-13       
[22] latticeExtra_0.6-28 RColorBrewer_1.1-2  lattice_0.20-33    
  • Your shapefiles aren't formatted differently, they have different CRSs! The second one is probably in UTM, the first in lat-lon. Use spTransform to convert one of them to the CRS of the other, then repeat the gDifference operation. Can you share your original data (zip and put in public Dropbox)? – Simbamangu Oct 13 '16 at 9:22
  • Yes - I figured that out - sorry I'm still learning geoproccessing and the correct terms. I used the spTransform to have both shapefiles in the same CRS and was then able to use your example on my shapefiles. Thanks – user41509 Oct 13 '16 at 17:59
1

You need gDifference from the rgeos package.

library(rgeos)
# Create SpatialPolygons objects
polygon1 <- readWKT("POLYGON((-180 -20, -140 55, 10 0, -140 -60, -180 -20))")
polygon2 <- readWKT("POLYGON((-180 -20, -140 70, 10 0, -140 -60, -180 -20))")

plot(polygon2, col = 'lightgray', border = 'darkgreen')
plot(polygon1, col = NA, border = 'red', add = T)

# gIntersection will find the overlapping parts:
test <- gIntersection(polygon1, polygon2)
plot(test, add = T, col = 'blue')

# gDifference will find the NON-overlapping parts:
test2 <- gDifference(polygon2, polygon1)
plot(test2, add = T, col = 'red')

You can then combine the polygons with polygon3 <- gUnion(test2, polygon1)

  • Thank you - I was able to figure things out with my shapefiles and apply for answer to get everything to work. – user41509 Oct 12 '16 at 15:10

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.