[R-sig-Geo] Self-intersection error when using intersect with two shapefiles
DAnnunzio, Remi (FOA)
Remi.DAnnunzio at fao.org
Tue Nov 1 15:47:02 CET 2016
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://cdn.sstatic.net/Sites/stackoverflow/img/apple-touch-icon@2.png?v=73d79a89bded&a]<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://i0.wp.com/www.gis-blog.com/wp-content/uploads/2015/02/srt_austria_tot.jpeg?fit=641%2C656]<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 at r-project.org>> on behalf of Janka VANSCHOENWINKEL <janka.vanschoenwinkel at 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://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>
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://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>
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<http://www.uhasselt.be/musicforlife>>
P Please consider the environment before printing this e-mail
[[alternative HTML version deleted]]
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list