This is to demonstrate the spatial distribution of sample surveys from Yen’s ecolocal research in Sapa, Lao Cai, Vietnam
# Reading shapefile in R
library(raster)
## Loading required package: sp
## mappoint<-shapefile(file.choose())
# However,data used in this session was converted to gegraphic points and here below
df<-read.csv("https://raw.githubusercontent.com/tuyenhavan/Statistics/master/Ecologicalpointdata.csv")
# This df data is very massy and needs to be tidied up
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
df1<-df %>% dplyr::select(x=coords.x1,y=coords.x2)
df2<-df1 # Making a copy of it
# Transforming from geographic coordinate dataframe to spatial point data
coordinates(df2)<-~x+y
proj4string(df2)<-CRS("+proj=longlat +datum=WGS84") # Assining geographic coordinates to it
Getting DEM for the study area from the internet
library(ggmap)
## Warning: package 'ggmap' was built under R version 3.4.2
## Loading required package: ggplot2
library(raster)
VN<-getData("SRTM",lon=104,lat=22,download = T)
plot(VN)
Extracting elevation values from DEM using point data
df_extract<-extract(VN,df2,df=T)
# Constructing a dataframe with geographic coordinates and elevation values
df3<-data.frame(df1,df_extract$srtm_57_08)
names(df3)[3]<-"Elevation"
head(df3)
## x y Elevation
## 1 104.0281 22.21885 678
## 2 104.0403 22.20656 889
## 3 104.0406 22.20730 869
## 4 103.7695 22.34774 1915
## 5 103.7695 22.34807 1906
## 6 103.7700 22.34844 1894
Visualzing point data
# Reconstract a spatial data point with elevation data
finalmap<-df3 # Making a copy of it
coordinates(finalmap)<-~x+y
proj4string(finalmap)<-CRS("+proj=longlat +datum=WGS84")
head(finalmap@data)
## Elevation
## 1 678
## 2 889
## 3 869
## 4 1915
## 5 1906
## 6 1894
plot(finalmap,col=4,main="Spatial Point Distribution",pch=16,axes=T)
Overlaying this map into google earth
library(ggmap)
myvn<-get_map(location = c(lon=103.863603,lat=22.297664),maptype = "terrain",zoom=8)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=22.297664,103.863603&zoom=8&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
ggmap(myvn,base_layer = ggplot(data=df3)+geom_point(aes(x=x,y=y))) # I don't know why it does not work here :)
Alternative approach
library(leaflet)
leaflet(df3) %>% addTiles() %>% addProviderTiles(providers$Thunderforest.Landscape) %>% addMarkers(lng=~x,lat=~y,popup = ~as.character(Elevation),label = ~as.character(Elevation))