[R-sig-Geo] Self-intersection error when using intersect with two shapefiles
Janka VANSCHOENWINKEL
janka.vanschoenwinkel at uhasselt.be
Tue Nov 1 15:07:29 CET 2016
Dear Remi,
Thank you very much for your reply.
Unfortunately, the first solution does not do it for me. I have two
shapefiles. One of them contains raster-style data but they are saved as a
shapefile. So I want to overlap two shapefiles. The extract function only
works with a real raster and a shapefile if I understand it correctly. I
tried it, and it didn't work.
For the second solution: I have never heard of this and I will have a
thorough look at it.
Janka
2016-11-01 10:54 GMT+01:00 DAnnunzio, Remi (FOA) <Remi.DAnnunzio at fao.org>:
> 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). 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://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
>
>
--
[image: Logo UHasselt] Mevrouw Janka Vanschoenwinkel
*Doctoraatsbursaal - PhD *
Milieueconomie - Environmental economics
T +32(0)11 26 86 96 | GSM +32(0)476 28 21 40
www.uhasselt.be/eec
Universiteit Hasselt | Campus Diepenbeek
Agoralaan Gebouw D | B-3590 Diepenbeek
Kantoor F11
Postadres: Universiteit Hasselt | Martelarenlaan 42 | B-3500 Hasselt
[image: Music For Life] Maak van UHasselt de #warmsteunief |
www.uhasselt.be/musicforlife
P Please consider the environment before printing this e-mail
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list