In this example I will use QGIS geoprocessing scripts through the RQGIS library in R. I will use data from the 2005 Peru Demographic and Health Survey, where I construct Voronoi polygons from the primary sampling unit locations and generate estimates of educational attainment for women.
library(sf)## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3library(RQGIS) #you need Qgis 2.18 installed, not 3!## Loading required package: reticulatelibrary(mapview)eth_dhs_points<-st_read("C:/Users/HP/Google Drive/fikre/ethiopia_dhs/ETGE71FL",layer = "ETGE71FL")## Reading layer `ETGE71FL' from data source `C:\Users\HP\Google Drive\fikre\ethiopia_dhs\ETGE71FL' using driver `ESRI Shapefile'
## Simple feature collection with 643 features and 20 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 5.684342e-14 ymin: 5.684342e-14 xmax: 47.00792 ymax: 14.50213
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defseths_dhs_points<-st_transform(eth_dhs_points, crs=20136 )
#project the data into a meter based system
mapview(eth_dhs_points["DHSCLUST"])This lets R find where your QGIS binaries are located. set_env() should work without you specifying the path
my_env<-set_env(root = "c:/Program Files/QGIS 2.18/")First we find the name of the algorithm for point buffering
find_algorithms(search_term = "buffer", qgis_env = my_env)## [1] "Fixed distance buffer-------------------------------->qgis:fixeddistancebuffer"
## [2] "Variable distance buffer----------------------------->qgis:variabledistancebuffer"
## [3] "r.buffer - Creates a raster map layer showing buffer zones surrounding cells that contain non-NULL category values.--->grass7:r.buffer"
## [4] "r.buffer.lowmem - Creates a raster map layer showing buffer zones surrounding cells that contain non-NULL category values (low-memory alternative).--->grass7:r.buffer.lowmem"
## [5] "v.buffer.column - Creates a buffer around features of given type.--->grass7:v.buffer.column"
## [6] "v.buffer.distance - Creates a buffer around features of given type.--->grass7:v.buffer.distance"
## [7] "Fixed distance buffer-------------------------------->saga:shapesbufferfixeddistance"
## [8] "Raster buffer---------------------------------------->saga:gridbuffer"
## [9] "Raster proximity buffer------------------------------>saga:gridproximitybuffer"
## [10] "Threshold raster buffer------------------------------>saga:thresholdbuffer"
## [11] "Variable distance buffer----------------------------->saga:shapesbufferattributedistance"
## [12] "Buffer vectors--------------------------------------->gdalogr:buffervectors"
## [13] "Single sided buffers (and offset lines) for lines---->gdalogr:singlesidedbuffersandoffsetlinesforlines"
## [14] "r.buffer - Creates a raster map layer showing buffer zones surrounding cells that contain non-NULL category values.--->grass:r.buffer"
## [15] "v.buffer.column - Creates a buffer around features of given type.--->grass:v.buffer.column"
## [16] "v.buffer.distance - Creates a buffer around features of given type.--->grass:v.buffer.distance"For this, we’ll use the qgis:fixeddistancebuffer function, but we need to see what the arguments to the function are:
get_usage(alg="qgis:fixeddistancebuffer", qgis_env = my_env, intern = F)## ALGORITHM: Fixed distance buffer
## INPUT <ParameterVector>
## DISTANCE <ParameterNumber>
## SEGMENTS <ParameterNumber>
## DISSOLVE <ParameterBoolean>
## OUTPUT <OutputVector>so we have 5 arguments, when we use the function, we need to specify all of these:
params <- get_args_man(alg = "qgis:fixeddistancebuffer", qgis_env = my_env)
params## $INPUT
## [1] "None"
##
## $DISTANCE
## [1] "10.0"
##
## $SEGMENTS
## [1] "5"
##
## $DISSOLVE
## [1] "False"
##
## $OUTPUT
## [1] "None"Here I do a 5km buffer around each PSU location.
wd<-"C:/Users/HP/Google Drive/fikre"
params$INPUT <-eths_dhs_points
params$DISTANCE<-5000 #5km around each point
params$OUTPUT<-file.path(wd, "eths_psu_buffer_5k.shp") # path to the output shapefilenow we have our parameters defined, we run the script:
eths_buff <- run_qgis(alg = "qgis:fixeddistancebuffer", params = params, load_output = TRUE, qgis_env = my_env)## $OUTPUT
## [1] "C:/Users/HP/Google Drive/fikre/eths_psu_buffer_5k.shp"mapview(eths_buff["DHSCLUST"], legend=T,map.types="OpenStreetMap")5km buffers done! Now, If I only had some other point data, I could do some point in polygon operations!
find_algorithms(search_term = "Voronoi")## [1] "Voronoi polygons------------------------------------->qgis:voronoipolygons"
## [2] "v.voronoi - Creates a Voronoi diagram from an input vector layer containing points.--->grass7:v.voronoi"
## [3] "v.voronoi - Creates a Voronoi diagram from an input vector layer containing points.--->grass:v.voronoi"get_usage(alg="qgis:voronoipolygons", qgis_env = my_env, intern = F)## ALGORITHM: Voronoi polygons
## INPUT <ParameterVector>
## BUFFER <ParameterNumber>
## OUTPUT <OutputVector>so we have 3 arguments, when we use the function,:
params <- get_args_man(alg = "qgis:voronoipolygons", qgis_env = my_env)
params## $INPUT
## [1] "None"
##
## $BUFFER
## [1] "0.0"
##
## $OUTPUT
## [1] "None"#wd<-"/Users/fikrewoldbitew/Google Drive/fikre"
params$INPUT <-eths_dhs_points
params$OUTPUT<-file.path(wd, "eths_psu_von_poly.shp") # path to the output shapefilenow we have our parameters defined, we run the script:
eths_von <- run_qgis(alg = "qgis:voronoipolygons",
params = params,
load_output = TRUE,
qgis_env = my_env)## $OUTPUT
## [1] "C:/Users/HP/Google Drive/fikre/eths_psu_von_poly.shp"mapview(eths_von["DHSCLUST"],map.types="OpenStreetMap")The polygons obviously have issues, so we can clip them to the Peruvian border:
find_algorithms(search_term = "clip")## [1] "Clip------------------------------------------------->qgis:clip"
## [2] "Clip Data-------------------------------------------->lidartools:clipdata"
## [3] "Poly Clip Data--------------------------------------->lidartools:polyclipdata"
## [4] "lasclip---------------------------------------------->lidartools:lasclip"
## [5] "Clip points with polygons---------------------------->saga:clippointswithpolygons"
## [6] "Clip raster with polygon----------------------------->saga:clipgridwithpolygon"
## [7] "Polygon clipping------------------------------------->saga:polygonclipping"
## [8] "Clip raster by extent-------------------------------->gdalogr:cliprasterbyextent"
## [9] "Clip raster by mask layer---------------------------->gdalogr:cliprasterbymasklayer"
## [10] "Clip vectors by extent------------------------------->gdalogr:clipvectorsbyextent"
## [11] "Clip vectors by polygon------------------------------>gdalogr:clipvectorsbypolygon"get_usage(alg="qgis:clip", qgis_env = my_env, intern = F)## ALGORITHM: Clip
## INPUT <ParameterVector>
## OVERLAY <ParameterVector>
## OUTPUT <OutputVector>so we have 3 arguments, when we use the function,:
params <- get_args_man(alg = "qgis:clip", qgis_env = my_env)
params## $INPUT
## [1] "None"
##
## $OVERLAY
## [1] "None"
##
## $OUTPUT
## [1] "None"library(maptools)## Loading required package: sp## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()eths_border=sf::st_read("C:/Users/HP/Google Drive/fikre/ETH_adm_DIVA_GIS/ETH_adm0.shp")## Reading layer `ETH_adm0' from data source `C:\Users\HP\Google Drive\fikre\ETH_adm_DIVA_GIS\ETH_adm0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 70 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 32.99994 ymin: 3.322099 xmax: 47.98618 ymax: 14.89996
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defseths_border<-st_transform(eths_border,crs=20136) #your CRS didn't match the one you used above
wd= "C:/Users/HP/Google Drive/fikre"
params$INPUT <-eths_von
params$OVERLAY<-eths_border
params$OUTPUT<-file.path(wd, "eth_psu_von_poly_clip.shp") # path to the output shapefilenow we have our parameters defined, we run the script:
eths_von_clip <- run_qgis(alg = "qgis:clip",
params = params,
load_output = TRUE,
qgis_env = my_env)## $OUTPUT
## [1] "C:/Users/HP/Google Drive/fikre/eth_psu_von_poly_clip.shp"mapview(eths_von_clip["DHSCLUST"])