[R-sig-Geo] Self-intersection error when using intersect with two shapefiles
DAnnunzio, Remi (FOA)
Remi.DAnnunzio at fao.org
Tue Nov 1 10:54:46 CET 2016
Dear Janka
1/ I would suggest to make things simpler, i.e use the original raster data directly with the shapefile you want to get your values from and use the extract() function. See post here http://stackoverflow.com/questions/22333473/how-do-i-extract-raster-values-from-polygon-data-then-join-into-spatial-data-fra
2/ it seems like you do have a geometry issue at the mentioned point in your NUTS3 shapefile. Maybe you should correct it in your favorite GIS environment, or get a better dataset with more detailed boundaries (e.g www.gadm.org<http://www.gadm.org>). You can also access the GADM data directly in R with the getData() function. See here for example http://www.gis-blog.com/r-raster-data-acquisition/.
Hope it helps
Rémi
Rémi d'Annunzio
Remote sensing and GIS - Forestry Department
UN-FAO
Phone: +44 7 94 64 35 698
Email: remi.dannunzio at fao.org
Skype: rdannunzio
________________________________
From: R-sig-Geo <r-sig-geo-bounces at r-project.org> on behalf of Janka VANSCHOENWINKEL <janka.vanschoenwinkel at uhasselt.be>
Sent: 01 November 2016 09:24
To: r-sig-geo at r-project.org
Subject: [R-sig-Geo] Self-intersection error when using intersect with two shapefiles
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
[http://cdn.sstatic.net/Sites/gis/img/apple-touch-icon@2.png?v=54e3ab1edcf3&a]<http://gis.stackexchange.com/questions/216095/r-self-intersection-error-when-using-intersect-with-two-shapefiles>
R: Self-intersection error when using intersect with two shapefiles<http://gis.stackexchange.com/questions/216095/r-self-intersection-error-when-using-intersect-with-two-shapefiles>
gis.stackexchange.com
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 inte...
with no response yet.
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
R-sig-Geo Info Page<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
stat.ethz.ch
R-sig-Geo -- R Special Interest Group on using Geographical data and Mapping About R-sig-Geo
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list