forked from oscarperpinan/spacetime-vis
-
Notifications
You must be signed in to change notification settings - Fork 0
/
physical.R
151 lines (123 loc) · 5.98 KB
/
physical.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
##################################################################
## Source code for the book: "Displaying time series, spatial and
## space-time data with R"
##
## Copyright (C) 2013-2012 Oscar Perpiñán Lamigueiro
##
## This program is free software you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published
## by the Free Software Foundation; either version 2 of the License,
## or (at your option) any later version.
##
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.
####################################################################
##################################################################
## Initial configuration
##################################################################
## Clone or download the repository and set the working directory
## with setwd to the folder where the repository is located.
##################################################################
## Physical maps
##################################################################
library(raster)
library(rasterVis)
library(maptools)
library(rgeos)
library(latticeExtra)
library(colorspace)
## Longitude-Latitude projection
proj <- CRS(' +proj=longlat +ellps=WGS84')
##################################################################
## Retrieving data from DIVA-GIS, GADM and Natural Earth Data
##################################################################
old <- setwd(tempdir())
download.file('http://biogeo.ucdavis.edu/data/gadm2/shp/BRA_adm.zip',
'BRA_adm.zip')
unzip('BRA_adm.zip')
brazilAdm <- readShapePoly('BRA_adm1.shp', proj4string=proj)
Encoding(levels(brazilAdm$NAME_1)) <- 'latin1'
download.file('http://biogeo.ucdavis.edu/data/diva/alt/BRA_alt.zip',
'BRA_alt.zip')
unzip('BRA_alt.zip')
brazilDEM <- raster('BRA_alt')
## World Water lines (Natural Earth)
download.file('http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_rivers_lake_centerlines.zip',
'neRivers.zip')
unzip('neRivers.zip')
worldlRiv <- readShapeLines('ne_10m_rivers_lake_centerlines', proj4string = proj)
download.file('http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/raster/OB_LR.zip',
'neSea.zip')
unzip('neSea.zip')
worldSea <- raster('OB_LR.tif')
brazilSea <- crop(worldSea, brazilDEM)
setwd(old)
##################################################################
## Intersection of shapefiles and elevation model
##################################################################
## only those features labeled as "River" are needed
worlRiv<- worlRiv[worlRiv$featurecla=='River',]
## Define the extent of Brazil as a SpatialPolygons
extBrazil <- as(extent(brazilDEM), 'SpatialPolygons')
proj4string(extBrazil) <- proj
## and intersect it with worldRiv to extract brazilian rivers
## from the world database
brazilRiv <- gIntersection(worldRiv, extBrazil, byid=TRUE)
## and especially the famous Amazonas River
amazonas <- worldRiv[worldRiv$name=='Amazonas',]
##################################################################
## Labels
##################################################################
## Locations of labels of each polygon
centroids <- coordinates(brazilAdm)
## Location of the "Brazil" label (average of the set of polygons centroids)
xyBrazil <- apply(centroids, 2, mean)
admNames <- strsplit(as.character(brazilAdm$NAME_1), ' ')
admNames <- sapply(admNames,
FUN=function(s){
sep=if (length(s)>2) '\n' else ' '
paste(s, collapse=sep)
})
##################################################################
## Overlaying layers of information
##################################################################
blueTheme <- rasterTheme(region=brewer.pal(n=9, 'Blues'))
seaPlot <- levelplot(brazilSea, par.settings=blueTheme,
maxpixels=1e6, panel=panel.levelplot.raster,
margin=FALSE, colorkey=FALSE)
terrainTheme <- rasterTheme(region=terrain_hcl(15))
altPlot <- levelplot(brazilDEM, par.settings=terrainTheme,
maxpixels=1e6, panel=panel.levelplot.raster,
margin=FALSE, colorkey=FALSE)
amazonasLab <- label(amazonas, 'Amazonas')
png(filename="figs/brazil.png",res=300,height=2000,width=2000)
seaPlot + altPlot + layer({
## Rivers
sp.lines(brazilRiv, col='darkblue', lwd=0.2)
## Amazonas
sp.lineLabel(amazonas, amazonasLab,
lwd=1, col='darkblue', col.line='darkblue',
cex=0.5, fontfamily='Palatino')
## Administrative boundaries
sp.polygons(brazilAdm, col='black', lwd=0.2)
## Centroids of administrative boundaries ...
panel.points(centroids, col='black')
## ... with their labels
panel.pointLabel(centroids, labels=admNames,
cex=0.7, fontfamily='Palatino', lineheight=.8)
## Country name
panel.text(xyBrazil[1], xyBrazil[2], labels='B R A Z I L',
cex=1.5, fontfamily = 'Palatino', fontface=2)
})
dev.off()
png(filename="figs/rastersBrazil.png",res=300,height=2000,width=2000)
print(seaPlot, split=c(1, 1, 2, 1), more=TRUE)
print(altPlot, split=c(2, 1, 2, 1))
dev.off()