R Spatial Data Analysis

This is to demonstrate the capability of R in spatial analysis

library(raster)
library(rgeos)
library(rgdal)
library(maptools)
library(sp)

Dataset will be used

data<-read.csv("https://raw.githubusercontent.com/tuyenhavan/Statistics/Dataset/viet_ppM.csv")
head(data)

Reading shape files in R

# Reading a shape file from your local computer as follow
## myshp<-shapefile(file.choose())
# Reading a shape file from web
library(ggmap)
VN<-get_map(location = "VN",zoom=5,maptype = "satellite")

Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=VN&zoom=5&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false
Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=VN&sensor=false

ggmap(VN)

Visualizing the shape file with numeric data

Vietmap<-getData(name = "GADM",country="Vietnam",level=1)

trying URL 'http://biogeo.ucdavis.edu/data/gadm2.8/rds/VNM_adm1.rds'
Content type 'application/gzip' length 2765340 bytes (2.6 MB)
downloaded 2.6 MB

spplot(Vietmap,zcol="ID_1",col="transparent",main="Vietnam Map",sub="Area of Vietnam")

Visualizing shapefile using color pallete

library(RColorBrewer)
display.brewer.all() # Display all available colors

spplot(Vietmap,"ID_1",col.regions=brewer.pal(n=8,name = "Blues"),cuts=10,col="transparent")

Interpolating point data using Thin Spline

df1<-data
library(ggplot2)
ggplot(data=data,aes(x=lon,y=lat)) + geom_point()

library(fields)

package 㤼㸱fields㤼㸲 was built under R version 3.4.2Loading required package: spam
package 㤼㸱spam㤼㸲 was built under R version 3.4.2Loading required package: dotCall64
package 㤼㸱dotCall64㤼㸲 was built under R version 3.4.2Loading required package: grid
Spam version 2.1-1 (2017-07-02) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction 
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

Attaching package: 㤼㸱spam㤼㸲

The following objects are masked from 㤼㸱package:base㤼㸲:

    backsolve, forwardsolve

Loading required package: maps

Attaching package: 㤼㸱maps㤼㸲

The following object is masked from 㤼㸱package:GISTools㤼㸲:

    map.scale

col1<-cbind(df1$lon,df1$lat)
dat<-list(X=col1,Y=df1$pptn)
thins<-Tps(dat$X,dat$Y)
surface(thins)
title("Map")

Interpolating point data using Inverse Distance Weighting (IDW)

# Getting summary
summary(df1)

      lon             lat              pptn            Month    
 Min.   :102.4   Min.   : 8.917   Min.   : 54.0   April   : 25  
 1st Qu.:105.2   1st Qu.:12.750   1st Qu.:105.2   August  : 25  
 Median :106.3   Median :17.167   Median :141.0   December: 25  
 Mean   :106.3   Mean   :16.902   Mean   :142.6   February: 25  
 3rd Qu.:107.8   3rd Qu.:21.417   3rd Qu.:176.0   January : 25  
 Max.   :109.3   Max.   :23.083   Max.   :221.0   July    : 25  
                                                  (Other) :150  

library(gstat)

package 㤼㸱gstat㤼㸲 was built under R version 3.4.2
Attaching package: 㤼㸱gstat㤼㸲

The following object is masked from 㤼㸱package:spatstat㤼㸲:

    idw

# Setting the boundary
## IDW- Inverse Distance Weighting
##1) Use a Rectangular Grid
library(gstat)
d=read.csv("https://raw.githubusercontent.com/tuyenhavan/Statistics/Dataset/viet_ppM.csv")
head(d)

names(d)[1:2]<-c("x","y")
library(ggplot2)
coordinates(d) <- ~ x + y
summary(d)

Object of class SpatialPointsDataFrame
Coordinates:
        min       max
x 102.41668 109.25002
y   8.91667  23.08334
Is projected: NA 
proj4string : [NA]
Number of points: 300
Data attributes:
      pptn            Month    
 Min.   : 54.0   April   : 25  
 1st Qu.:105.2   August  : 25  
 Median :141.0   December: 25  
 Mean   :142.6   February: 25  
 3rd Qu.:176.0   January : 25  
 Max.   :221.0   July    : 25  
                 (Other) :150  

x.range <- as.numeric(c(102.4, 109.25))  # min/max longitude of the interpolation area
y.range <- as.numeric(c(8.9, 23.1))  # min/max latitude of the interpolation area  
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1), y = seq(from = y.range[1], to = y.range[2], by = 0.1))  # expand points to grid
coordinates(grd) <- ~x + y #assign coordinates to grid
gridded(grd) <- TRUE ## Create SpatialPixel object
plot(grd, cex = 2, col = "grey")

#idw formulae
idw <- idw(formula = pptn ~ 1, locations = d, newdata = grd)  

[inverse distance weighted interpolation]

idwO = as.data.frame(idw)  
names(idwO)

[1] "x"         "y"         "var1.pred" "var1.var" 

names(idwO)[1:3] <- c("x", "y", "var1.pred")
ggplot() + geom_tile(data = idwO, aes(x = x, y = y, fill = var1.pred)) 

ggplot() + geom_tile(data = idwO, aes(x = x, y = y, fill = var1.pred))+ scale_fill_gradient(low = "blue", high = "orange")

viet=getData("GADM",country="Vietnam",level=1)
plot(viet, axes=T)

names(viet)

 [1] "OBJECTID"  "ID_0"      "ISO"       "NAME_0"    "ID_1"     
 [6] "NAME_1"    "HASC_1"    "CCN_1"     "CCA_1"     "TYPE_1"   
[11] "ENGTYPE_1" "NL_NAME_1" "VARNAME_1"

vietC <- fortify(viet, region = "VARNAME_1")
##fortify before
##displaying shapefile in ggplot2, i.i
#convert map data to data frame
ggplot() + geom_tile(data = idwO, aes(x = x, y = y, fill = var1.pred))+ scale_fill_gradient(low = "blue", high = "orange")+ 
  geom_path(data = vietC, aes(long, lat, group = group), colour = "black")

#2)IDW Using the Vornoi method-provide user defined grid
d=read.csv("https://raw.githubusercontent.com/tuyenhavan/Statistics/Dataset/viet_ppM.csv")
names(d)[1:2]<-c("x","y")
coordinates(d) <- ~ x + y
projection(d)=CRS("+proj=longlat +ellps=WGS84")
library(dismo)

package 㤼㸱dismo㤼㸲 was built under R version 3.4.2
Attaching package: 㤼㸱dismo㤼㸲

The following object is masked from 㤼㸱package:ggmap㤼㸲:

    geocode

The following object is masked from 㤼㸱package:spatstat㤼㸲:

    domain

v <- voronoi(d)
plot(v)

#summarizes spatial variables
va <- aggregate(viet) #sp package
#set boundaries to vietnam
vca <- intersect(v, va)

non identical CRS

spplot(vca, 'pptn', col.regions=rev(get_col_regions()))

#build a raster of vietnam with stated resolution
r <- raster(va, res=0.01) #1 degree=111 sq km
projection(r)=CRS("+proj=longlat +ellps=WGS84")
#rasterize polygon
vr <- rasterize(vca, r, 'pptn')
plot(vr)

library(gstat)
gs <- gstat(formula=pptn~1, locations=d)
idw <- interpolate(r, gs)

[inverse distance weighted interpolation]

idw_disp <- mask(idw, vr)
plot(idw_disp)

Interpolating point data using Kriging

#### KRIGING
library(rgdal)
library(raster)
library(dismo)
library(rgeos)
library(maptools)
library(gstat)
library(ggplot2)
library(dplyr)

Attaching package: 㤼㸱dplyr㤼㸲

The following objects are masked from 㤼㸱package:raster㤼㸲:

    intersect, select, union

The following object is masked from 㤼㸱package:MASS㤼㸲:

    select

The following objects are masked from 㤼㸱package:rgeos㤼㸲:

    intersect, setdiff, union

The following object is masked from 㤼㸱package:nlme㤼㸲:

    collapse

The following objects are masked from 㤼㸱package:stats㤼㸲:

    filter, lag

The following objects are masked from 㤼㸱package:base㤼㸲:

    intersect, setdiff, setequal, union

d=read.csv("https://raw.githubusercontent.com/tuyenhavan/Statistics/Dataset/viet_ppM.csv")
names(d)[1:2]<-c("x","y")
coordinates(d) <- ~ x + y
summary(d)

Object of class SpatialPointsDataFrame
Coordinates:
        min       max
x 102.41668 109.25002
y   8.91667  23.08334
Is projected: NA 
proj4string : [NA]
Number of points: 300
Data attributes:
      pptn            Month    
 Min.   : 54.0   April   : 25  
 1st Qu.:105.2   August  : 25  
 Median :141.0   December: 25  
 Mean   :142.6   February: 25  
 3rd Qu.:176.0   January : 25  
 Max.   :221.0   July    : 25  
                 (Other) :150  

  
x.range <- as.numeric(c(102.4, 109.25))  # min/max longitude of the interpolation area
y.range <- as.numeric(c(8.9, 23.1))  # min/max latitude of the interpolation area  
#Create a rectangular grid
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1), y = seq(from = y.range[1], 
    to = y.range[2], by = 0.1))  # expand points to grid
coordinates(grd) <- ~x + y #assign coordinates to grid
gridded(grd) <- TRUE ## Create SpatialPixel object
plot(grd, cex = 2, col = "grey")
points(d, pch = 1, col = "blue", cex = 1)

#Compute a variogram
#capture the spatial continuity in data
#spatial distribution of the response variable
v = variogram(log(pptn)~1, d) 
plot(v) #plot semi-variogram

#avlaible variogram models
vgm() 

#Fit the variogram model
v.fit=fit.variogram(v, vgm("Exp")) 
#predict unknown locations
krigeM <- krige(log(pptn) ~ 1, d, grd, model=v.fit)

[using ordinary kriging]

#display
krigeM %>% as.data.frame %>% ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
  scale_fill_gradient(low = "yellow", high="red") +theme_bw()

library(spatstat)
library(sp)
library(rgeos)
library(maptools)
library(GISTools)
library(ggmap)
library(rgeos)
library(raster)
area<-getData(name = "GADM",country="United Kingdom",level=1)

trying URL 'http://biogeo.ucdavis.edu/data/gadm2.8/rds/GBR_adm1.rds'
Content type 'application/gzip' length 1157357 bytes (1.1 MB)
downloaded 1.1 MB

plot(area)

s=read.csv("https://raw.githubusercontent.com/tuyenhavan/Statistics/master/open-plaques-United-Kingdom-2017-06-19.csv")
head(s)

s=na.omit(s) #remove NAs
#give spatial reference to the geo-location points
s$x <- s$longitude  # define x & y as longitude and latitude
s$y <- s$latitude
coordinates(s) = ~x + y
s <- remove.duplicates(s)
plot(s, col="red", pch=20)

## 1) BUILD A GENERIC DENSITY POINT
#get the spatial extent of the point data
summary(s)

Object of class SpatialPointsDataFrame
Coordinates:
       min      max
x -7.65396  1.74302
y  0.00000 57.59584
Is projected: NA 
proj4string : [NA]
Number of points: 2265
Data attributes:
       id                                     title     
 Min.   :    1   Charles Dickens blue plaque     :   7  
 1st Qu.:  589   Alan Mathison Turing blue plaque:   6  
 Median : 3181   Charles Darwin blue plaque      :   4  
 Mean   : 8675   John Wesley blue plaque         :   4  
 3rd Qu.: 9010   Dorothy L. Sayers blue plaque   :   3  
 Max.   :42885   Emily Davies blue plaque        :   3  
                 (Other)                         :2238  
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    inscription  
 'Estcourt' - 195 Holdenhurst Road. The final home of Inspector Frederick George Abberline 1843-1929. During his 29 years with the Metropolitan Police Abberline gained 84 commendations and awards and became well-known for his work on the case of Jack-The-Ripper                                                                                                                                                                                                                     :   1  
 'Father' Henry Willis 1821-1901 organ builder lived here                                                                                                                                                                                                                                                                                                                                                                                                                                 :   1  
  Don Redfern Memorial Bandstand  funded by  High Peak Borough Council    The Don Redfern Memorial Committee                                                                                                                                                                                                                                                                                                                                                                              :   1  
 "A very gallant gentleman". To commemorate Captain Lawrence E. G. Oates a member of Capt. Scott's expedition to the South Pole 1910-1912 a frequent visitor to Meanwoodside, the Oates family home. Died 17th March 1912                                                                                                                                                                                                                                                                 :   1  
 "Chance Meeting"  These sculptures were inspired by two legendary Liverpudlians.    Ken Dodd O.B.E.    One of Liverpool's greatest entertainers, bringing laughter and joy to millions for more than 50 years.    Bessie Braddock MP. 1899 - 1970.    Labour MP for Liverpool Exchange for over 24 years. She campaigned tirelessly to improve the conditions for her constituents. Awarded the Freedom of the City in 1970.    Unveiled by Ken Dodd in June 2009.    Sculptor Tom Murphy:   1  
 "Epworth Villa" 14 Gloucester Street Reverend John Wesley, AM preached here on 6 September 1776 thereby making it the first Methodist meeting house in Weymouth (Melcombe Regis)                                                                                                                                                                                                                                                                                                         :   1  
 (Other)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  :2259  
    latitude       longitude                 country    
 Min.   : 0.00   Min.   :-7.6540   United Kingdom:2265  
 1st Qu.:51.50   1st Qu.:-1.9625                        
 Median :51.53   Median :-0.2304                        
 Mean   :51.98   Mean   :-1.1251                        
 3rd Qu.:52.59   3rd Qu.:-0.1449                        
 Max.   :57.60   Max.   : 1.7430                        
                                                        
            area                  address        erected    
 London       :1055   High Street     :   9   Min.   :   1  
 Manchester   :  79   ?               :   6   1st Qu.:1984  
 Birmingham   :  78   Cathedral Street:   5   Median :1999  
 Oxford       :  49   Church Street   :   4   Mean   :1990  
 Wolverhampton:  47   Granada Studios :   4   3rd Qu.:2009  
 Belfast      :  40   Hyde Park       :   4   Max.   :2017  
 (Other)      : 917   (Other)         :2233                 
                                                           main_photo  
                                                                : 363  
 http://farm4.staticflickr.com/3603/3412004633_ffa34dd0c6_z.jpg :   2  
 http://farm1.static.flickr.com/179/368146620_c9acec875d_b.jpg  :   1  
 http://farm1.static.flickr.com/6/9146592_21f13e419a_b.jpg      :   1  
 http://farm1.staticflickr.com/17/22241637_349d4fff4d_z.jpg     :   1  
 http://farm1.staticflickr.com/22/34622156_3782c8e823_z.jpg?zz=1:   1  
 (Other)                                                        :1896  
       colour                         organisations 
 blue     :1663   ["English Heritage"]       : 388  
 green    : 158   ["London County Council"]  : 216  
 black    :  91   ["Greater London Council"] : 190  
 bronze   :  69   []                         : 172  
 film cell:  55   ["Ulster History Circle"]  :  82  
 brown    :  52   ["Manchester City Council"]:  64  
 (Other)  : 177   (Other)                    :1153  
            language   
 English        :2248  
 Welsh & English:  13  
 Welsh          :   2  
 Irish          :   1  
 Irish & English:   1  
                :   0  
 (Other)        :   0  
                                                                                                   series    
                                                                                                      :2146  
 Centenary Of Cinema 1996                                                                             :  56  
 English Heritage outside London                                                                      :  32  
 Islington People's Plaques                                                                           :   8  
 The Great Exhibition of the works of industry of all nations & Crystal Palace building Hyde Park 1851:   4  
 Transport Heritage Site 'Red Wheel'                                                                  :   4  
 (Other)                                                                                              :  15  
   series_ref   geolocated.  photographed. number_of_subjects
        :2205   false:   0   false: 359    Min.   : 1.000    
 003    :   1   true :2265   true :1906    1st Qu.: 1.000    
 005    :   1                              Median : 1.000    
 011    :   1                              Mean   : 1.255    
 013    :   1                              3rd Qu.: 1.000    
 021    :   1                              Max.   :76.000    
 (Other):  55                                                
               lead_subject_name lead_subject_born_in
 Charles Dickens        :   9    Min.   : 849        
 Alan Mathison Turing   :   8    1st Qu.:1809        
 Isambard Kingdom Brunel:   7    Median :1859        
 Charles Darwin         :   6    Mean   :1839        
 John Lennon            :   6    3rd Qu.:1893        
 John Wesley            :   6    Max.   :2011        
 (Other)                :2223                        
 lead_subject_died_in lead_subject_type     lead_subject_roles
 Min.   : 899               :   0       []           :  87    
 1st Qu.:1878         animal:   0       ["architect"]:  34    
 Median :1929         group :  14       ["artist"]   :  21    
 Mean   :1909         man   :1967       ["poet"]     :  21    
 3rd Qu.:1966         place :  66       ["composer"] :  18    
 Max.   :2016         thing :  21       ["writer"]   :  15    
                      woman : 197       (Other)      :2069    
                                           lead_subject_wikipedia
                                                      : 170      
 https://en.wikipedia.org/wiki/Charles_Dickens        :   9      
 https://en.wikipedia.org/wiki/Alan_Mathison_Turing   :   8      
 https://en.wikipedia.org/wiki/Isambard_Kingdom_Brunel:   7      
 https://en.wikipedia.org/wiki/Charles_Darwin         :   6      
 https://en.wikipedia.org/wiki/John_Lennon            :   6      
 (Other)                                              :2059      
                                          lead_subject_dbpedia
                                                    : 170     
 http://dbpedia.org/resource/Charles_Dickens        :   9     
 http://dbpedia.org/resource/Alan_Mathison_Turing   :   8     
 http://dbpedia.org/resource/Isambard_Kingdom_Brunel:   7     
 http://dbpedia.org/resource/Charles_Darwin         :   6     
 http://dbpedia.org/resource/John_Lennon            :   6     
 (Other)                                            :2059     
                                                                                      lead_subject_image
                                                                                               :1064    
 https://commons.wikimedia.org/wiki/Special:FilePath/Dickens_Gurney_head.jpg?width=640         :   9    
 https://commons.wikimedia.org/wiki/Special:FilePath/Alan_Turing_Aged_16.jpg?width=640         :   8    
 https://commons.wikimedia.org/wiki/Special:FilePath/IKBrunelChains.jpg?width=640              :   7    
 https://commons.wikimedia.org/wiki/Special:FilePath/Charles_Darwin_01.jpg?width=640           :   6    
 https://commons.wikimedia.org/wiki/Special:FilePath/John_Wesley_by_George_Romney.jpg?width=640:   6    
 (Other)                                                                                       :1165    
                                                                                                                                                                                                                                                        subjects   
 ["Charles Dickens|(1812-1870)|man|novelist, journalist, policeman, son of John Dickens, son of Elizabeth Dickens, and friend of Thomas Latimer"]                                                                                                           :   8  
 ["Alan Mathison Turing|(1912-1954)|man|code-breaker, pioneer of computer science, founder of computer science, cryptographer, creator of computer science, mathematician, Reader in Mathematics, and designer of The Bombe"]                               :   7  
 ["Charles Darwin|(1809-1882)|man|naturalist, father of Anne Elizabeth Darwin, author of The Origin of Species, traveller, adventurer, writer, friend of Gideon Lincecum, son of Robert Darwin, father of George Howard Darwin, and husband of Emma Darwin"]:   5  
 ["flying bomb (V1/V2)|(1944-1945)|thing|bomb"]                                                                                                                                                                                                             :   5  
 ["John Wesley|(1703-1791)|man|evangelist, founder of Methodism, Master of Arts, son of Susanna Annesley, brother of Charles Wesley, and son of Samuel Wesley"]                                                                                             :   5  
 ["Edward Elgar|(1857-1934)|man|composer, Knight Bachelor, friend of Alfred Rodewald, friend of George Robertson Sinclair, organist, violinist, conductor, teacher, and husband of Caroline Alice Elgar"]                                                   :   4  
 (Other)                                                                                                                                                                                                                                                    :2231  

#Convert spatial data to point pattern object
#general form: ppp(x.coordinates, y.coordinates, x.range, y.range)
mypattern <- ppp(s$longitude, s$latitude, c(-7 ,1), c(0,57))

51 points were rejected as lying outside the specified window

plot(mypattern)

plot(density(mypattern, sigma = 500)) # sigma sets the diamater of the kernel in map units

## try changing the kernel diameter to 700
##2) BUILD A DENSITY PLOT WHICH MAPS SPATIAL CONCENTRATION AT UK-WIDE SCALE
window <- as.owin(area)
plot(window)

mypattern <- ppp(s$longitude, s$latitude, window=window)

11 points were rejected as lying outside the specified window

plot(density(mypattern, sigma = 500))