[R-sig-Geo] Self-intersection error when using intersect with two shapefiles
Janka VANSCHOENWINKEL
janka.vanschoenwinkel at uhasselt.be
Tue Nov 1 22:47:58 CET 2016
Dear Remi,
Awesome!
I first got the following error: "identicalCRS(x,y) is not TRUE"
I therefore added a projection change with the following code: grid <-
spTransform(grid, CRSobj = CRS(proj4string(nuts)))
That did the job!
Thank you very much for the effort to write a small example code. That made
a huge difference!
Janka
2016-11-01 15:47 GMT+01:00 DAnnunzio, Remi (FOA) <Remi.DAnnunzio at fao.org>:
> This below may do what you are looking for, to be adapted
>
> # 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 at 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 at 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)
>
>
> *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
> DAnnunzio, Remi (FOA) <Remi.DAnnunzio at fao.org>
> *Sent:* 01 November 2016 14:16
> *To:* Janka VANSCHOENWINKEL
> *Cc:* r-sig-geo at r-project.org
> *Subject:* Re: [R-sig-Geo] Self-intersection error when using intersect
> with two shapefiles
>
> Dear Janka,
>
> Then if you have two vector based layers, use the over() function
>
> Run vignette("over") for detailed explanations
>
>
> Best
>
> Remi
>
> PS : for info, the http://biogeo.ucdavis.edu/ that hosts the datasets
> from getData() is often down these days. expect issues...
> Ecology, Geography, and Agriculture - Biogeographic Research
> <http://biogeo.ucdavis.edu/>
> biogeo.ucdavis.edu
> Ecology, Geography, and Agriculture¶ We use quantitative methods to study
> crop ecology, agricultural biodiversity, and agricultural development and
> their effects on ...
>
>
>
> 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: Janka VANSCHOENWINKEL <janka.vanschoenwinkel at uhasselt.be>
> Sent: 01 November 2016 14:07
> To: DAnnunzio, Remi (FOA)
> Cc: r-sig-geo at r-project.org
> Subject: Re: [R-sig-Geo] Self-intersection error when using intersect with
> two shapefiles
>
> 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<
> mailto: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
>
> <http://stackoverflow.com/questions/22333473/how-do-i-extract-raster-values-from-polygon-data-then-join-into-spatial-data-fra>
> How do I extract raster values from polygon data then join ...
> <http://stackoverflow.com/questions/22333473/how-do-i-extract-raster-values-from-polygon-data-then-join-into-spatial-data-fra>
> stackoverflow.com
> I would like to merge polygon data and raster data into one dataframe for
> purposes of then using randomForests package in R. This involves first
> extracting the mean ...
>
>
>
> 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/.
> Global Administrative Areas | Boundaries without limits
> <http://www.gadm.org/>
> www.gadm.org
> GADM is a spatial database of the location of the world's administrative
> areas (or adminstrative boundaries) for use in GIS and similar software.
> <http://www.gis-blog.com/r-raster-data-acquisition/>
> R raster: getData SRTM, WorldClim, GADM | GIS-Blog.com
> <http://www.gis-blog.com/r-raster-data-acquisition/>
> www.gis-blog.com
> Hey! Today I will show how powerful the R {raster} package is on another
> example. The raster package is not only a great tool for raster processing
> and calculation ...
>
>
>
> 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<mailto:remi.dannunzio at fao.org>
> Skype: rdannunzio
>
>
>
> ________________________________
> From: R-sig-Geo <r-sig-geo-bounces at r-project.org<mailto:r-sig-geo-bounces@
> r-project.org>> on behalf of Janka VANSCHOENWINKEL <janka.vanschoenwinkel@
> uhasselt.be<mailto:janka.vanschoenwinkel at uhasselt.be>>
> Sent: 01 November 2016 09:24
> To: r-sig-geo at r-project.org<mailto: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...
>
> [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>
>
> <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...
>
> gis.stackexchange.com<http://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<mailto: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
>
> R-sig-Geo Info Page<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
>
> stat.ethz.ch<http://stat.ethz.ch>
> R-sig-Geo -- R Special Interest Group on using Geographical data and
> Mapping About R-sig-Geo
>
>
>
>
>
> --
>
> [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<http://www.uhasselt.be/eec>
>
> Universiteit Hasselt | Campus Diepenbeek
> Agoralaan Gebouw D | B-3590 Diepenbeek
> Kantoor F11
>
> Postadres: Universiteit Hasselt | Martelarenlaan 42 | B-3500 Hasselt
>
>
> [Music For Life] Maak van UHasselt de #warmsteunief |
> www.uhasselt.be/musicforlife<http://www.uhasselt.be/musicforlife
> <http://www.uhasselt.be/musicforlife%3Chttp://www.uhasselt.be/musicforlife>
> >
>
>
>
> P Please consider the environment before printing this e-mail
>
>
> [[alternative HTML version deleted]]
>
>
--
[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