0

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
0

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.