[R-sig-Geo] Self-intersection error when using intersect with two shapefiles

Janka VANSCHOENWINKEL janka.vanschoenwinkel at uhasselt.be
Tue Nov 1 10:24:01 CET 2016


Dear colleagues,

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 (based on the worldwide raster data).

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

The full code can be found below. I have used it before with other
shapefiles and it worked. With the new shapefiles that I am using now,
I get the above error. I think it might be a projection problem but I
have not a lot of experience with this.

Your help is appreciated a lot! Thanks in advance!

Janka


**Below you can find the full code + data + session info**

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

grid_nuts <- gIntersects(new_spdf,nuts,byid = TRUE)
#Take all the values in a particular NUTS polygon, multiply by each
the areas of the
#corresponding new grid cells then divide by the total area of the NUTS polygon

#Added the if statement to avoid non-overlapping NUTS polygons

for(i in 1:length(nuts)){
  if(any(grid_nuts[i,])){
    nuts at data$average_spatial_value[i] <-
mean(new_spdf at data$value[grid_nuts[i,]]*

gArea(new_spdf[grid_nuts[i,],])/
                                                 gArea(nuts[i,]))
  } else {
    nuts at data$average_spatial_value[i] <- NA
  }

}



> sessionInfo() # one of the things I tried to solve the above problem is to install the latest R versions

R version 3.3.1 (2016-06-21)
Platform: x86_64-w64-mingw32/x64 (64-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] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rgeos_0.3-21 raster_2.5-8 rgdal_1.1-10 sp_1.2-3

loaded via a namespace (and not attached):
[1] tools_3.3.1     Rcpp_0.12.7     grid_3.3.1      lattice_0.20-33


P.S. I also posted the question on
http://gis.stackexchange.com/questions/216095/r-self-intersection-error-when-using-intersect-with-two-shapefiles
with no response yet.



More information about the R-sig-Geo mailing list