User Tools

Site Tools



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

Link to this comparison view

wiki:rstat:r_load_rast [2017/12/05 22:53]
wiki:rstat:r_load_rast [2020/07/17 06:39] (current)
Line 1: Line 1:
 +====== Raster data processing in R - load and plot maps ======
 +[[wiki:​rstat:​r_geography|Back to spatial analyses in R]]
 +<note warning>​WORK IN PROGRESS DO NOT FOLLOW THIS PAGE</​note>​
 +We can create a temporary folder to process data frm this tutorial migrate to this folder and enter R from there so owr wR working directory will be this one
 +    mkdir visualR
 +    cd visualR
 +    R
 +===== Raster data visualization =====
 +We  start using two libraries raster and rgdal:
 +    library(rgdal) ​ # link to Gdal/Ogr functions with R to import data
 +    library(raster) # More spatial data functionality ex: create raster stack layer or compute and Moran I
 +Using the **system** R function we can access bash commands and download data using **wget** function:
 +    system("​wget http://​​ost4sem/​exercise/​r_ecoservice/​raster/​ZONES_m9.asc"​)
 +Where are the data stored now? can you see them? They are in yoyr system folder visualR not loaded in R. To load them and plot them in R we can use readGDAL and raster plot function:
 +<code R>
 +plot(MapFromRasterLib) # Visualize the map and legend
 +plot(MapFromRasterLib,​ col=c("​blue","​red","​green3","​orangeplot function of r basic package"​),​ legend=F, axes=F, main="​Map of zones",​cex.main=3,​ cex=2, cex.legend=3,​ box=F)
 +legend(131000,​110000,​c("​Coastal","​West","​Centre","​East"​),​ bty = "​n", ​ pch=c(15,​15,​15,​15),​ pt.cex=c(2.5,​2.5,​2.5), ​ cex=3, ​ col=c("​blue","​red","​green3","​orange"​))
 +{{:​wiki:​rstat:​r_map_cornwall.png?​200|Output from raster plot function}} This is the output of the plot function of the raster package which is different from the plot function from R basic package. \\
 +=== Multiple layers plotting in R  ===
 +The script below preform the following in a loop. The first task is to create a list object and a vector object of names, the n for each layer we do the following:
 +  * download layer in working directory
 +  * load layer layer from working directory to the R List 
 +Finally we are ready to store layers from the list to a RasterBrick object available from the raster library and plot the brick with default raster color scheme.
 +<code R>
 +#- Plot environmental services and zones of Cornwall
 +import12=vector("​list",​12) ​
 +for ( x in c(1:12) ){ 
 +system(paste0("​wget http://​​ost4sem/​exercise/​r_ecoservice/​raster/",​map_names[x]))
 +Here what you should see.\\
 +Now we can do some data cleaning (redefining NA values) and plotting the maps into a png file saved outside R with a specific color palette and a legend:
 +<code R>
 +length(which(complete.cases(import12[[1]]@data@values)==FALSE)) # complete.cases return a logical vector indicating which cases are complete i.e. have no missing values.
 +import12[[5]]@data@values[which(complete.cases(import12[[1]]@data@values)==FALSE)]=NA # data cleaning
 +mycolor[1]="#​a6a6a6"​ # or 808080 are grey
 +TITLE=c("​Agriculture","​Carbon above ground","​Carbon in soil","​Flood mitigation","​Plant production","​Solar energy","​Aesthetic","​Recreation","​Tourism","​Urban development","​Wind energy","​Zones of Cornwall"​)
 +png(paste("​priority_maps.png",​sep=""​),​ width = 1600, height = 1200)
 +nf <- layout(matrix(c(1:​12),​3,​4,​byrow=TRUE),​ c(1,1), c(1,1), TRUE) ; ; par(mar=c(5,​5,​3,​5))
 +for ( x in c(1,​2,​3,​4,​5,​7,​8,​9,​10,​11,​6) ){ 
 +plot(r, col=mycolor,​ legend=F, axes=F, main=TITLE[x],​cex.main=3,​ cex=2, cex.legend=3,​ box=F)
 +if (x==6){plot(r,​ legend.only=TRUE,​ col=mycolor,​ legend.width=6,​ legend.shrink=0.75,​ axis.args=list(at=c(0,​50,​89), ​ labels=c(0,​50,​100), ​ cex.axis=2.2))}
 +plot(r, col=c("​blue","​red","​green3","​orange"​),​ legend=F, axes=F, main="​Zones of Cornwall",​cex.main=3,​ cex=1.8, box=F)
 +legend(131000,​110000,​c("​Coastal","​West","​Centre","​East"​),​ bty = "​n", ​ pch=c(15,​15,​15,​15),​ pt.cex=c(2.5,​2.5,​2.5), ​ cex=3, ​ col=c("​blue","​red","​green3","​orange"​))
 +The output map:\\