[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