Important packages to load

pacman::p_load(pacman, rio, sp, rgdal, rgeos, raster, tmap, ggplot2, ggmap, leaflet, spatstat, gstat, dplyr)

#sp - provide classes and methods for spatial data #rgdal & gdal - importing and exporting geospatial data formats #rgeos - handling operations on topologies #raster - working with raster data #tmap - creating thematic maps #ggplot - data visualization #ggmap - add base maps such as google maps, open street maps #leaflet - for interactive maps #spatstat - spatial point pattern analysis #gstat - Geostatistical modeling

The crime data set is loaded with ggmap package

sfdis <- import("F:/My Documents and files/Megasync/TWData/R/sfsdm.xlsx")

str(sfdis)
## 'data.frame':    30 obs. of  5 variables:
##  $ year: chr  "2017" "2017" "2017" "2017" ...
##  $ ref : chr  "TW" "TW" "TW" "TW" ...
##  $ loc : chr  "Polpithigama-Pansiyagama-Dr's land" "Maho-FA's homegarden" "Galgamuwa Termite mound" "Galgamuwa-CBNTsite" ...
##  $ lat : chr  "7°43'48.2\"N" "7°53'48.3\"N" "8°05'59.4\"N" "8°06'00.3\"N" ...
##  $ lon : chr  "80°30'59.4\"E" "80°14'48.2\"E" "80°16'51.5\"E" "80°16'50.5\"E" ...
head(sfdis, n=5)
##   year ref                                loc         lat          lon
## 1 2017  TW Polpithigama-Pansiyagama-Dr's land 7°43'48.2"N 80°30'59.4"E
## 2 2017  TW               Maho-FA's homegarden 7°53'48.3"N 80°14'48.2"E
## 3 2017  TW            Galgamuwa Termite mound 8°05'59.4"N 80°16'51.5"E
## 4 2017  TW                 Galgamuwa-CBNTsite 8°06'00.3"N 80°16'50.5"E
## 5 2017  TW               Maho - Termite mound 7°53'54.0"N 80°14'42.2"E
summary(sfdis)
##      year               ref                loc                lat           
##  Length:30          Length:30          Length:30          Length:30         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##      lon           
##  Length:30         
##  Class :character  
##  Mode  :character

Filtering data

sfdis <- filter(sfdis, lon != "NA") # Selecting only the data without NA s for lon using filter function from dplyr

Converting coordinates to decimal WGS84 system

#defining separators
dmslat <- char2dms(sfdis$lat, chd = "°", chm = "'", chs = "\"")
latnum <- as.numeric(dmslat)
 
dmslon <- char2dms(sfdis$lon, chd = "°", chm = "'", chs = "\"")
lonnum <- as.numeric(dmslon)

#Creating a final dataset
sfdis_df <- cbind(sfdis[,1:3], latnum, lonnum)

#naming column names
colnames(sfdis_df) <- c("year", "ref", "loc", "lat", "lon")

This is a local R data frame. Some packages may not understand them. Therefore, need to convert this into a spatial data frame

coords <- SpatialPoints(sfdis_df[, c("lon", "lat")]) #says these are the coordinates
sfdis_dfsp <- SpatialPointsDataFrame(coords, sfdis_df)#converting to spatial points
proj4string(sfdis_dfsp) <- CRS("+proj=longlat +no_defs +ellps=WGS84", doCheckCRSArgs = F) #This set up the projection. It gives an error if didnt add "doCheckCRSArgs = F"

class(sfdis_dfsp)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

Make sure the class of the data set is ‘SpatialPointsDataFrame’ before further analysis

SpatialPointsDataFrame has 5 slots

  1. Data - all the data we include
  2. coords.nrs - the column positions where in data the coordinates were taken from
  3. coords - coordinates generated above spatialpoints function.
  4. bbox bounding box of the coordinates
  5. proj4string Coordinate reference system

CRS() function above assign a referencing sytem

There are 2 main types of CRS’s

  1. Geographic Coordinate System (longitude, latitude) - WGS84
  2. Cartesian/Projected/Planar coordinate system (x,y)

There are some others such as NADS83, UTM (Universal Transverse Mercator)

Extracting data slot - must use @ instead of $

Saving the spatial file

saveRDS(sfdis_dfsp, "F:/My Documents and files/Megasync/TWData/R/TWshapes/sfdis_dfsp.rds")

We can save it as ESRI shapefiles to be able to use with other GIS systems.

writeOGR(sfdis_dfsp, dsn = "F:/My Documents and files/Megasync/TWData/R/TWshapes/sfdistribution.shp", driver = "ESRI Shapefile",layer = "sfdistribution", overwrite_layer = TRUE)

Exploration of data

bbox(sfdis_dfsp) #check boundaries
##           min      max
## lon 79.690333 81.12719
## lat  5.944194  9.75000
proj4string(sfdis_dfsp) #check coordinate system
## [1] "+proj=longlat +no_defs +ellps=WGS84"

Visualizing Data

by plot() function or External R libraries which are useful to create beautiful maps.

. ggplot2 . plotly . ggmap . tmap *. leaflet

sl_map <- readOGR(dsn = "F:/My Documents and files/Megasync/TWData/GIS/Shapes", layer = "lka_admbnda_adm2_slsd_20200305")
## OGR data source with driver: ESRI Shapefile 
## Source: "F:\My Documents and files\Megasync\TWData\GIS\Shapes", layer: "lka_admbnda_adm2_slsd_20200305"
## with 26 features
## It has 24 fields
class(sl_map)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
plot(sl_map, col = "lightgreen", axes = TRUE)
plot(sfdis_dfsp, pch = 21, bg = "red", cex = 1, add = TRUE)

Adding a title and a legend

plot(sl_map, col = "lightgreen", axes = TRUE)
plot(sfdis_dfsp, pch = 21, bg = "red", cex = 1, add = TRUE)
title("Phlebotomus argentipes records in Sri Lanka")
legend ("topright", title = "Legend", legend = "Recorded location", pch = 21, pt.bg = "red", pt.cex = 1, bty = "n")