Self-intersection error when using intersect with two shapefiles

Tue, Nov 1, 2016 2:54 AM

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

Self-intersection error when using intersect with two shapefiles

Thread (6 messages)