R Polygon Self-Intersection – Fix Self-Intersection Error When Using Intersect with Two Shapefiles in R

intersectpolygonrself-intersectionshapefile

I am trying to "translate" worldwide raster data to European nuts3 polygons. That means I try to overlap two shapefiles to find the average value on nuts3 level.

Unfortunately, when I use the intersect function, I get the following error:

new_spdf <- intersect(grid,nuts)
    Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  : 
      TopologyException: Input geom 1 is invalid: Self-intersection at or near point 13.3081265 48.608032999999999 at 13.3081265 48.608032999999999

Below you can find the full code. 

### download the files available through this link:
# https://drive.google.com/open?id=0By9u5m3kxn9yUy1xVDF2NV9TajA

> library(rgdal)
> nuts <- readOGR(".", layer = "NUTS_RG_60M_2010")
OGR data source with driver: ESRI Shapefile 
Source: ".", layer: "NUTS_RG_60M_2010"
with 1920 features
It has 4 fields
> grid <- readOGR(".", layer = "a3_mean_annual_runoff_1950_2000")
OGR data source with driver: ESRI Shapefile 
Source: ".", layer: "a3_mean_annual_runoff_1950_2000"
with 41553 features
It has 2 fields
> summary(nuts)
Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x -61.74369 55.83663
y -21.37656 71.17308
Is projected: FALSE 
proj4string : [+proj=longlat +ellps=GRS80 +no_defs]
Data attributes:
    NUTS_ID       STAT_LEVL_     SHAPE_Leng         SHAPE_Area      
 AT     :   1   Min.   :0.00   Min.   :  0.1326   Min.   : 0.00057  
 AT1    :   1   1st Qu.:3.00   1st Qu.:  1.6480   1st Qu.: 0.11795  
 AT11   :   1   Median :3.00   Median :  2.9772   Median : 0.35707  
 AT111  :   1   Mean   :2.66   Mean   :  5.1573   Mean   : 1.56752  
 AT112  :   1   3rd Qu.:3.00   3rd Qu.:  5.3525   3rd Qu.: 0.97955  
 AT113  :   1   Max.   :3.00   Max.   :132.3810   Max.   :81.09899  
 (Other):1914                                                       
> summary(grid)
Object of class SpatialPolygonsDataFrame
Coordinates:
   min   max
x -180 180.0
y  -55  83.5
Is projected: FALSE 
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
       ID             CODE       
 Min.   :    1   Min.   :   0.0  
 1st Qu.:10389   1st Qu.: 118.6  
 Median :20777   Median : 260.0  
 Mean   :20777   Mean   : 432.1  
 3rd Qu.:31165   3rd Qu.: 528.8  
 Max.   :41553   Max.   :6854.2
grid <- spTransform(grid, CRSobj = CRS(proj4string(nuts)))
> new_spdf <- intersect(grid,nuts)
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  : 
  TopologyException: Input geom 1 is invalid: Self-intersection at or near point 13.3081265 48.608032999999999 at 13.3081265 48.608032999999999

Best Answer

I found the answer through the R-sig-geo forum (I'll post the link when it is available).

#projection
grid <- spTransform(grid, CRSobj = CRS(proj4string(nuts)))    

# Select only the elements of grid that intersect nuts

grid_nuts <- grid[nuts,]

# Count the number of elements from grid_nuts that fall within each zone of nuts
nuts@data$count_grid <- unlist(over(x =nuts ,y=grid_nuts[,"ID"],fn="length"))

# Compute the average of the value CODE of elements from grid_nuts that fall within each zone of nuts
nuts@data$mean_grid <-  unlist(over(x =nuts ,y=grid_nuts[,"CODE"],fn="mean"))

# Export
writeOGR(obj=nuts,dsn="nuts_grid.shp",layer="nuts_grid",driver="ESRI Shapefile",overwrite_layer = T)
Related Question