library(rgdal)
## Loading required package: sp
## rgdal: version: 1.2-8, (SVN revision 663)
## 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/Jorn/Documents/R/win-library/3.4/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: C:/Users/Jorn/Documents/R/win-library/3.4/rgdal/proj
## Linking to sp version: 1.2-5
library(maptools)
## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()
library(sp)
#specify the working directory dsn and the layer (without extension .shp)
wd <- "C:\\Users\\Jorn\\Documents\\R_training" # wd <- "F:/R_Training" setting the work space
setwd(wd)
dsn<-getwd()
#check the projection system - will look for the prj file to retrieve projection information
ogrInfo(dsn=dsn,layer="TLA")
## Source: "C:/Users/Jorn/Documents/R_training", layer: "TLA"
## Driver: ESRI Shapefile; number of rows: 9
## Feature type: wkbPolygon with 2 dimensions
## Extent: (1691730 5480255) - (1920888 5757382)
## CRS: +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs
## LDID: 0
## Number of fields: 4
## name type length typeName
## 1 TA2014 4 254 String
## 2 TA2014_NAM 4 254 String
## 3 SHAPE_Leng 2 19 Real
## 4 SHAPE_Area 2 19 Real
#read the actual shapefile, do not add the .shp
prov.ogr<-readOGR(dsn=dsn,layer="TLA")
## OGR data source with driver: ESRI Shapefile
## Source: "C:/Users/Jorn/Documents/R_training", layer: "TLA"
## with 9 features
## It has 4 fields
summary(prov.ogr)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 1691730 1920888
## y 5480255 5757382
## Is projected: TRUE
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000
## +y_0=10000000 +ellps=GRS80 +units=m +no_defs]
## Data attributes:
## TA2014 TA2014_NAM SHAPE_Leng
## 021 :1 Horowhenua District :1 Min. :108406
## 034 :1 Manawatu District :1 1st Qu.:317559
## 036 :1 Palmerston North City:1 Median :422018
## 037 :1 Rangitikei District :1 Mean :400731
## 038 :1 Ruapehu District :1 3rd Qu.:492256
## 039 :1 Stratford District :1 Max. :646203
## (Other):3 (Other) :3
## SHAPE_Area
## Min. :3.946e+08
## 1st Qu.:2.163e+09
## Median :2.567e+09
## Mean :3.457e+09
## 3rd Qu.:4.484e+09
## Max. :6.964e+09
##
class(prov.ogr)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
# with maptools - specify the projection system yourself, do not add the .shp extension (source: https://gis.stackexchange.com/questions/20389/converting-nzmg-or-nztm-to-latitude-longitude-for-use-with-r-map-library/116682)
prov<-readShapePoly("TLA",proj4string=CRS("+proj=tmerc +lat_0=0.0 +lon_0=173.0 +k=0.9996 +x_0=1600000.0 +y_0=10000000.0 +datum=WGS84 +units=m"))
## Warning: use rgdal::readOGR or sf::st_read
summary(prov)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 1691730 1920888
## y 5480255 5757382
## Is projected: TRUE
## proj4string :
## [+proj=tmerc +lat_0=0.0 +lon_0=173.0 +k=0.9996 +x_0=1600000.0
## +y_0=10000000.0 +datum=WGS84 +units=m +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## TA2014 TA2014_NAM SHAPE_Leng
## 021 :1 Horowhenua District :1 Min. :108406
## 034 :1 Manawatu District :1 1st Qu.:317559
## 036 :1 Palmerston North City:1 Median :422018
## 037 :1 Rangitikei District :1 Mean :400731
## 038 :1 Ruapehu District :1 3rd Qu.:492256
## 039 :1 Stratford District :1 Max. :646203
## (Other):3 (Other) :3
## SHAPE_Area
## Min. :3.946e+08
## 1st Qu.:2.163e+09
## Median :2.567e+09
## Mean :3.457e+09
## 3rd Qu.:4.484e+09
## Max. :6.964e+09
##
class(prov)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
#Plotting shapefiles
#using spplot, define the sp-object, the variable to be displayed
#optionally you can specify scales, north arrows etc
spplot(prov,c("TA2014_NAM"),scales=list(draw=TRUE,cex=0.75))
#in order to change the colors you can specify them yourself or use predefined colors specified in the package RColorBrewer
library(RColorBrewer)
#1. Sequential palettes are suited to ordered data that progress from low to high. Lightness steps dominate the look of these schemes, with light colors for low data values to dark colors for high data values.
#2. Diverging palettes put equal emphasis on mid-range critical values and extremes at both ends of the data range. The critical class or break in the middle of the legend is emphasized with light colors and low and high extremes are emphasized with dark colors that have contrasting hues.
#3. Qualitative palettes do not imply magnitude differences between legend classes, and hues are used to create the primary visual differences between classes. Qualitative schemes are best suited to representing nominal or categorical data
display.brewer.all(n=NULL, type="seq", select=NULL, exact.n=TRUE)
display.brewer.all(n=NULL, type="div", select=NULL, exact.n=TRUE)
display.brewer.all(n=NULL, type="qual", select=NULL, exact.n=TRUE)
spplot(prov,c("TA2014_NAM"),col.regions=brewer.pal(9,"PuBuGn"),scales=list(draw=TRUE,cex=0.75))
spplot(prov,c("TA2014_NAM"),col.regions=brewer.pal(10,"RdBu"),scales=list(draw=TRUE,cex=0.75))
spplot(prov,c("TA2014_NAM"),col.regions=brewer.pal(10,"Set3"),scales=list(draw=TRUE,cex=0.75))
#manual selection of the colors
colcodes<-c("red","green","yellow","purple","blue","lightblue","cyan","darkgrey","lightgrey","white")
spplot(prov,c("TA2014_NAM"),col.regions=colcodes,scales=list(draw=TRUE,cex=0.75))
#you can also draw the shapefiles using the plot() commands. Plot can overlay the different shapefiles, but spplot cannot
plot(prov,col=brewer.pal(10,"Set3"),axes=T,cex.axis=0.75)
reg<-readShapePoly("Regional_Council_Boundaries_2011NZTM",proj4string=CRS("+proj=tmerc +lat_0=0.0 +lon_0=173.0 +k=0.9996 +x_0=1600000.0 +y_0=10000000.0 +datum=WGS84 +units=m"))
## Warning: use rgdal::readOGR or sf::st_read
plot(reg,add=T)
#add some text and legend
text(coordinates(prov),labels=prov$TA2014_NAM,cex=.75)
legend("bottomleft",legend=prov$TA2014_NAM,fill=brewer.pal(10,"Set3"),cex=0.75,bty="n")
