User Tools

Site Tools



This shows you the differences between two versions of the page.

Link to this comparison view

wiki:rstat:r_load_vect [2017/12/05 22:53]
wiki:rstat:r_load_vect [2020/07/17 06:39] (current)
Line 1: Line 1:
 +====== Load, plot vector layers & do spatial joins in R ======
 +[[wiki:​rstat:​r_geography|Back to spatial analyses in R]]
 +===== Vector data visualization =====
 +We have followed a [[http://​​2013/​12/​introduction-spatial-data-ggplot2/​ |]] tutorial to replicate and improve our case study. Here we work on vector layers, first download a compressed folder in linux environment and extract it and check what you have using bash functions within R.\\
 +   ​system("​wget http://​​ost4sem/​exercise/​r_ecoservice/​vector/​"​)
 +   ​system("​unzip"​)
 +   ​system("​ogrinfo -al cornwall_ward_27700.shp | more ")
 +   # ... 
 +   # Layer SRS WKT:
 +   # (unknown)
 +We downloaded a vector shape file in EPSG 27700 but this info is not available in the metadata, so we can update the metadata with projection parameters in a new layer.\\
 +   ​system("​ogr2ogr -s_srs EPSG:27700 -t_srs EPSG:27700 cornwall.shp cornwall_ward_27700.shp"​)
 +Using the raster library we can load shape files in R from linux OS folder.
 +   ​library(raster)
 +   ​ <- shapefile("​cornwall_ward_27700.shp"​)
 +   ward <- shapefile("​cornwall.shp"​)
 +   # Tell the difference between the two layers
 +   ​summary(
 +   ​summary(ward)
 +<note tip>​[[http://​​packages/​cran/​rgdal/​docs/​ogrInfo| More]] readings on using **rgdal** library to load vector maps</​note>​\\
 +The raster library has a plot function able to plot R spatial objects such as the SpatialPolygonsDataFrame which is the corresponding object to the vector shape file we just imported. raster plot is overwriting the basic plot fuction of R.\\
 +    plot(ward)
 +The ward administrative boudary in cornwall\\
 +===== Spatial Join in R =====
 +We can load  a text table in csv format, do a join by identifier and plot maps with statistics.\\
 +   stat = read.csv("​SYNT_2002.csv",​ header=T)
 +To join files we need a common identifier. Which is the table and vector file common ID? Find it and merge the files using **merge** function from raster library to then plot the map and save the join vector file in your working directory outside R.\\ 
 +    str(stat)
 +    str(ward)
 +    agri <- merge(ward, stat, by.x='​PI',​ by.y='​ID'​)
 +    writeOGR(agri,​ "​.",​ "​agristat",​ driver="​ESRI Shapefile"​)
 +Now we can plot in the ward map some features about agriculture using the statistics we have joined to our spatial database using the **ggplot2** library.\\
 +<note tip>​[[http://​​bbs/​wp-content/​uploads/​2011/​09/​handout_ggplot2.pdf|a ggplot2 tutorial]] and [[https://​​wp-content/​uploads/​2015/​03/​ggplot2-cheatsheet.pdf| ggplot2 cheatsheet ]]</​note>​
 +Plotting some features of our database (ex: ward surface vs numberof diary livestock per ward) in black and white and in colour.
 +   ​library(ggplot2)
 +   # simple visualization black and white
 +   ​p<​-ggplot(agri@data,​ aes(SUM_HA,​dairy_v))
 +   ​p+geom_point()
 +   # Add some colours
 +   ​p+geom_point(colour="​red",​ size=2)
 +   # Add a legend
 +   ​p+geom_point(aes(colour=SUM_HA,​ size=dairy_v))
 +   # Add labels ​ (these are too many for our cvisualization)
 +   ​p+geom_point(aes(colour=SUM_HA,​size=dairy_v))+geom_text(size=8,​aes(label=FIRST_NM))
 +Play with this and plot other informations...
 +Now we map the number of Diary livestock in Cornwall per each ward. We use fortify ( check ?fortify in R).
 +   ​agri_geom<​-fortify(agri,​ region="​PI"​)
 +   ​agri_geom<​-merge(agri_geom,​ agri@data, by.x="​id",​ by.y="​PI"​)
 +   ​head(agri_geom)
 +   ## The map function:
 +   ## Map<​-ggplot(agri_geom,​ aes(long,​lat,​ group=group,​ fill=dairy_v))
 +   ## + geom_polygon()+ coord_equal() + 
 +   ## labs(x="​Easting (m)", y="​Northing (m)",​fill= "# of livestocks"​)
 +   ## + ggtitle ("​Dairy cow in Cornwall"​)
 +   ​Map<​-ggplot(agri_geom,​ aes(long,​lat,​ group=group,​ fill=dairy_v))+ geom_polygon()+ coord_equal() + labs(x="​Easting (m)", y="​Northing (m)",​fill= "# of livestocks"​)+ ggtitle ("​Dairy cow in Cornwall"​)
 +   Map
 +   # Black and white map
 +   Map + scale_fill_gradient(low="​white",​ high="​black"​)