Friday, 22 April 2016

polygon - R: Self-intersection error when using intersect with two shapefiles


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

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)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...