library(sp)
## Warning: package 'sp' was built under R version 3.2.5
# get spatial data classes
getClass("Spatial")
## Class "Spatial" [package "sp"]
##
## Slots:
##
## Name: bbox proj4string
## Class: matrix CRS
##
## Known Subclasses:
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2
library(maps)
## Warning: package 'maps' was built under R version 3.2.5
##
## # maps v3.1: updated 'world': all lakes moved to separate new #
## # 'lakes' database. Type '?world' or 'news(package="maps")'. #
data("canada.cities")
canada<-canada.cities
head(canada)
## name country.etc pop lat long capital
## 1 Abbotsford BC BC 157795 49.06 -122.30 0
## 2 Acton ON ON 8308 43.63 -80.03 0
## 3 Acton Vale QC QC 5153 45.63 -72.57 0
## 4 Airdrie AB AB 25863 51.30 -114.02 0
## 5 Aklavik NT NT 643 68.22 -135.00 0
## 6 Albanel QC QC 1090 48.87 -72.42 0
# Assign coordinate systems
coordinates(canada)<-~long +lat
# assign geographic coordinate system
proj4string(canada)<-CRS("+proj=longlat +datum=WGS84")
# plot the data point
plot(canada)
# Get coordinate data points
head(coordinates(canada))
## long lat
## 1 -122.30 49.06
## 2 -80.03 43.63
## 3 -72.57 45.63
## 4 -114.02 51.30
## 5 -135.00 68.22
## 6 -72.42 48.87
# adding a new column to spatial dataframe
# creating a column with population>20 thousand
canada$pop_25<-canada@data$pop>25000
head(canada@data)
## name country.etc pop capital pop_25
## 1 Abbotsford BC BC 157795 0 TRUE
## 2 Acton ON ON 8308 0 FALSE
## 3 Acton Vale QC QC 5153 0 FALSE
## 4 Airdrie AB AB 25863 0 TRUE
## 5 Aklavik NT NT 643 0 FALSE
## 6 Albanel QC QC 1090 0 FALSE
# Summarize the total number of counties with pop >25 thoudsand
sum(canada@data$pop_25)
## [1] 88
# Import and export filters
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.2.5
## rgdal: version: 1.2-5, (SVN revision 648)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: C:/Users/Tuyen/Documents/R/win-library/3.2/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: C:/Users/Tuyen/Documents/R/win-library/3.2/rgdal/proj
## Linking to sp version: 1.2-4
ogrDrivers()
## name long_name
## 1 AeronavFAA Aeronav FAA
## 2 ARCGEN Arc/Info Generate
## 3 AVCBin Arc/Info Binary Coverage
## 4 AVCE00 Arc/Info E00 (ASCII) Coverage
## 5 BNA Atlas BNA
## 6 CSV Comma Separated Value (.csv)
## 7 DGN Microstation DGN
## 8 DXF AutoCAD DXF
## 9 EDIGEO French EDIGEO exchange format
## 10 ESRI Shapefile ESRI Shapefile
## 11 Geoconcept Geoconcept
## 12 GeoJSON GeoJSON
## 13 Geomedia Geomedia .mdb
## 14 GeoRSS GeoRSS
## 15 GML Geography Markup Language (GML)
## 16 GPKG GeoPackage
## 17 GPSBabel GPSBabel
## 18 GPSTrackMaker GPSTrackMaker
## 19 GPX GPX
## 20 HTF Hydrographic Transfer Vector
## 21 Idrisi Idrisi Vector (.vct)
## 22 JML OpenJUMP JML
## 23 KML Keyhole Markup Language (KML)
## 24 MapInfo File MapInfo File
## 25 Memory Memory
## 26 MSSQLSpatial Microsoft SQL Server Spatial Database
## 27 ODBC ODBC
## 28 ODS Open Document/ LibreOffice / OpenOffice Spreadsheet
## 29 OGR_GMT GMT ASCII Vectors (.gmt)
## 30 OGR_PDS Planetary Data Systems TABLE
## 31 OGR_SDTS SDTS
## 32 OGR_VRT VRT - Virtual Datasource
## 33 OpenAir OpenAir
## 34 OpenFileGDB ESRI FileGDB
## 35 OSM OpenStreetMap XML and PBF
## 36 PCIDSK PCIDSK Database File
## 37 PDF Geospatial PDF
## 38 PGDUMP PostgreSQL SQL dump
## 39 PGeo ESRI Personal GeoDatabase
## 40 REC EPIInfo .REC
## 41 S57 IHO S-57 (ENC)
## 42 SEGUKOOA SEG-P1 / UKOOA P1/90
## 43 SEGY SEG-Y
## 44 Selafin Selafin
## 45 SQLite SQLite / Spatialite
## 46 SUA Tim Newport-Peace's Special Use Airspace Format
## 47 SVG Scalable Vector Graphics
## 48 SXF Storage and eXchange Format
## 49 TIGER U.S. Census TIGER/Line
## 50 UK .NTF UK .NTF
## 51 VFK Czech Cadastral Exchange Data Format
## 52 Walk Walk
## 53 WAsP WAsP .map format
## 54 XLSX MS Office Open XML spreadsheet
## 55 XPlane X-Plane/Flightgear aeronautical data
## write copy isVector
## 1 FALSE FALSE TRUE
## 2 FALSE FALSE TRUE
## 3 FALSE FALSE TRUE
## 4 FALSE FALSE TRUE
## 5 TRUE FALSE TRUE
## 6 TRUE FALSE TRUE
## 7 TRUE FALSE TRUE
## 8 TRUE FALSE TRUE
## 9 FALSE FALSE TRUE
## 10 TRUE FALSE TRUE
## 11 TRUE FALSE TRUE
## 12 TRUE FALSE TRUE
## 13 FALSE FALSE TRUE
## 14 TRUE FALSE TRUE
## 15 TRUE FALSE TRUE
## 16 TRUE TRUE TRUE
## 17 TRUE FALSE TRUE
## 18 TRUE FALSE TRUE
## 19 TRUE FALSE TRUE
## 20 FALSE FALSE TRUE
## 21 FALSE FALSE TRUE
## 22 TRUE FALSE TRUE
## 23 TRUE FALSE TRUE
## 24 TRUE FALSE TRUE
## 25 TRUE FALSE TRUE
## 26 TRUE FALSE TRUE
## 27 TRUE FALSE TRUE
## 28 TRUE FALSE TRUE
## 29 TRUE FALSE TRUE
## 30 FALSE FALSE TRUE
## 31 FALSE FALSE TRUE
## 32 FALSE FALSE TRUE
## 33 FALSE FALSE TRUE
## 34 FALSE FALSE TRUE
## 35 FALSE FALSE TRUE
## 36 TRUE FALSE TRUE
## 37 TRUE TRUE TRUE
## 38 TRUE FALSE TRUE
## 39 FALSE FALSE TRUE
## 40 FALSE FALSE TRUE
## 41 TRUE FALSE TRUE
## 42 FALSE FALSE TRUE
## 43 FALSE FALSE TRUE
## 44 TRUE FALSE TRUE
## 45 TRUE FALSE TRUE
## 46 FALSE FALSE TRUE
## 47 FALSE FALSE TRUE
## 48 FALSE FALSE TRUE
## 49 TRUE FALSE TRUE
## 50 FALSE FALSE TRUE
## 51 FALSE FALSE TRUE
## 52 FALSE FALSE TRUE
## 53 TRUE FALSE TRUE
## 54 TRUE FALSE TRUE
## 55 FALSE FALSE TRUE
# read shapefile from raster package
library(raster)
## Warning: package 'raster' was built under R version 3.2.5
mymap<-shapefile(system.file("external/lux.shp",package="raster"))
plot(mymap,col=palette(),axes=T,asp=1) # asp=1 means that setting up world coordinates for graphic window
# get online database
VN<-getData(name="GADM",country="Vietnam",level=1)
plot(VN,axes=T,col=rainbow(nrow(VN)))
# Functions from rgeos package
# union, distance, intersection, convex hull, envelope, buffer, simplify, polygon assembly, valid, area, length, intersects, touches, disjoint, crosses, within, contains, overlaps, equals, covers, etc.
# Subset one province of Vietnam, let's say "Vinh Phuc"
vp<-VN[VN@data$VARNAME_1=="Vinh Phuc" ,]
par(mar=rep(2,4))
plot(vp,axes=T)
# checking coordinates
crs(vp)
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
# transforming longlat to utm
vp_utm<-spTransform(vp,CRS("+proj=utm +zone=48 +datum=WGS84"))
#making buffer 1 km
library(rgeos)
## Warning: package 'rgeos' was built under R version 3.2.5
## rgeos version: 0.3-22, (SVN revision 544)
## GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084
## Linking to sp version: 1.2-4
## Polygon checking: TRUE
vp_buffer<-gBuffer(vp_utm,width = 1000)
plot(vp_utm,border=3)
plot(vp_buffer,border=4,add=T,lty=2)
# writable ogrDrivers
x <- ogrDrivers()
as.character(x[x$write==TRUE,"name"])
## [1] "BNA" "CSV" "DGN" "DXF"
## [5] "ESRI Shapefile" "Geoconcept" "GeoJSON" "GeoRSS"
## [9] "GML" "GPKG" "GPSBabel" "GPSTrackMaker"
## [13] "GPX" "JML" "KML" "MapInfo File"
## [17] "Memory" "MSSQLSpatial" "ODBC" "ODS"
## [21] "OGR_GMT" "PCIDSK" "PDF" "PGDUMP"
## [25] "S57" "Selafin" "SQLite" "TIGER"
## [29] "WAsP" "XLSX"
# Raster data: get online data
vn_raster <- getData('SRTM', lon=105,lat=21, download=TRUE)
# clip vinh phuc
vp_extent<-c(105.3217, 105.7858, 21.12178, 21.57398)
vp_crop<-crop(vn_raster,vp_extent,snap="out")
# mask and get Vinh Phuc
vp_mask<-mask(vp_crop,vp)
plot(vp_mask)
# Derive slope, aspect and hillshade from DEM
vp_slope<-terrain(vp_mask,opt="slope")
plot(vp_slope)
# aspect
vp_aspect<-terrain(vp_mask,opt="aspect")
plot(vp_aspect)
# hillshade
vp_tpi<-terrain(vp_mask,40,270,opt="TPI")
plot(vp_tpi)
# histogram
hist(getValues(vp_slope),col=4)
# get slope >0.5
plot(vp_slope>0.5)
# transforming raster coordinates
library(raster)
vp_raster_utm<-projectRaster(vp_mask,crs="+proj=utm +zone=48 +datum=WGS84")
# resample raster
resam<-aggregate(vp_mask,fact=11,fun=mean)
plot(resam)
