【独家发布】Displaying Time Series, Spatial, and Space-Time Data with R

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Displaying Time Series, Spatial, and Space-Time Data with R

https://github.com/oscarperpinan/spacetime-vis
http://www.crcpress.com/product/isbn/9781466565203

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：Time Series displaying Display Playing Spatial Series

1. airStations <- read.csv2('data/airStations.csv')
   
2. coordinates(airStations) <- ~ long + lat
   
3. ## Geographical projection
   
4. -proj4string(airStations) <- CRS("+proj=longlat +ellps=WGS84")
   
5. +proj4string(airStations) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
   
6. 
   
7. ## Measurements data
   
8. airQuality <- read.csv2('data/airQuality.csv')

复制代码

1. library(rgdal)
   
2. writeOGR(NO2sp, 'data/NO2.geojson', 'NO2sp', driver='GeoJSON')
   
3. 
   
4. +library(rgdal)
   
5. +writeOGR(NO2sp, dsn='NO2_mean.kml', layer='mean', driver='KML')
   
6. +
   
7. +library(plotKML)
   
8. +plotKML(NO2sp["mean"], points_names=NO2sp$codEst)
   
9. +
   
10. ##################################################################
    
11. ## Additional information with tooltips and hyperlinks
    
12. ##################################################################

复制代码

1. ##################################################################
   
2. ## Source code for the book: "Displaying time series, spatial and
   
3. ## space-time data with R"
   
4. 
   
5. ## Copyright (C) 2013-2012 Oscar Perpiñán Lamigueiro
   
6. 
   
7. ## This program is free software you can redistribute it and/or modify
   
8. ## it under the terms of the GNU General Public License as published
   
9. ## by the Free Software Foundation; either version 2 of the License,
   
10. ## or (at your option) any later version.
    
11. 
    
12. ## This program is distributed in the hope that it will be useful, but
    
13. ## WITHOUT ANY WARRANTY; without even the implied warranty of
    
14. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    
15. ## General Public License for more details.
    
16. 
    
17. ## You should have received a copy of the GNU General Public License
    
18. ## along with this program; if not, write to the Free Software
    
19. ## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
    
20. ## 02111-1307, USA.
    
21. ####################################################################
    
22. 
    
23. ##################################################################
    
24. ## Initial configuration
    
25. ##################################################################
    
26. ## Clone or download the repository and set the working directory
    
27. ## with setwd to the folder where the repository is located.
    
28. 
    
29. ##################################################################
    
30. ## Air Quality in Madrid
    
31. ##################################################################
    
32. 
    
33. ## codeStations.csv is extracted from the document
    
34. ## http://www.mambiente.munimadrid.es/opencms/export/sites/default/calaire/Anexos/INTPHORA-DIA.pdf,
    
35. ## table of page 3.
    
36. 
    
37. codEstaciones <- read.csv2('data/codeStations.csv')
    
38. codURL <- as.numeric(substr(codEstaciones$Codigo, 7, 8))
    
39. 
    
40. ## The information of each measuring station is available at its own webpage, defined by codURL
    
41. URLs <- paste('http://www.mambiente.munimadrid.es/opencms/opencms/calaire/contenidos/estaciones/estacion', codURL, '.html', sep='')
    
42. 
    
43. ##################################################################
    
44. ## Data arrangement
    
45. ##################################################################
    
46. 
    
47. library(XML)
    
48. library(sp)
    
49. 
    
50. ## Access each webpage, retrieve tables and extract long/lat data
    
51. coords <- lapply(URLs, function(est){
    
52.   tables <- readHTMLTable(est)
    
53.   location <- tables[[2]]
    
54.   ## Clean the table content and convert to dms format
    
55.   ub2dms <- function(x){
    
56.     ch <- as.character(x)
    
57.     ch <- sub(',', '.', ch)
    
58.     ch <- sub('O', 'W', ch) ## Some stations use "O" instead of "W"
    
59.     as.numeric(char2dms(ch, "o", "'", "'' "))
    
60.   }
    
61.   long <- ub2dms(location[2,1])
    
62.   lat <- ub2dms(location[2,2])
    
63.   alt <- as.numeric(sub(' m.', '', location[2, 3]))
    
64. 
    
65.   coords <- data.frame(long=long, lat=lat, alt=alt)
    
66. 
    
67.   coords
    
68. })
    
69. 
    
70. airStations <- cbind(codEstaciones, do.call(rbind, coords))
    
71. 
    
72. ## The longitude of "El Pardo" station is wrong (positive instead of negative)
    
73. airStations$long[22] <- -airStations$long[22]
    
74. 
    
75. write.csv2(airStations, file='data/airStations.csv')
    
76. 
    
77. ## Fill in the form at
    
78. ## http://www.mambiente.munimadrid.es/opencms/opencms/calaire/consulta/descarga.html
    
79. ## to receive the Diarios11.zip file.
    
80. unzip('data/Diarios11.zip')
    
81. rawData <- readLines('data/Datos11.txt')
    
82. ## This loop reads each line and extracts fields as defined by the
    
83. ## INTPHORA file:
    
84. ## http://www.mambiente.munimadrid.es/opencms/export/sites/default/calaire/Anexos/INTPHORA-DIA.pdf
    
85. datos11 <- lapply(rawData, function(x){
    
86.   codEst <- substr(x, 1, 8)
    
87.   codParam <- substr(x, 9, 10)
    
88.   codTec <- substr(x, 11, 12)
    
89.   codPeriod <- substr(x, 13, 14)
    
90.   month <- substr(x, 17, 18)
    
91.   dat <- substr(x, 19, nchar(x))
    
92.   ## "N" used for impossible days (31st April)
    
93.   idxN <- gregexpr('N', dat)[[1]]
    
94.   if (idxN==-1) idxN <- numeric(0)
    
95.   nZeroDays <- length(idxN)
    
96.   day <- seq(1, 31-nZeroDays)
    
97.   ## Substitute V and N with ";" to split data from different days
    
98.   dat <- gsub('[VN]+', ';', dat)
    
99.   dat <- as.numeric(strsplit(dat, ';')[[1]])
    
100.   ## Only data from valid days
     
101.   dat <- dat[day]
     
102.   res <- data.frame(codEst, codParam, ##codTec, codPeriod,
     
103.                     month, day, year=2011,
     
104.                     dat)
     
105.   })
     
106. datos11 <- do.call(rbind, datos11)
     
107. write.csv2(datos11, 'data/airQuality.csv')
     
108. 
     
109. ##################################################################
     
110. ## Combine data and spatial locations
     
111. ##################################################################
     
112. 
     
113. library(sp)
     
114. 
     
115. ## Spatial location of stations
     
116. airStations <- read.csv2('data/airStations.csv')
     
117. coordinates(airStations) <- ~ long + lat
     
118. ## Geographical projection
     
119. proj4string(airStations) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
     
120. 
     
121. ## Measurements data
     
122. airQuality <- read.csv2('data/airQuality.csv')
     
123. ## Only interested in NO2
     
124. NO2 <- airQuality[airQuality$codParam==8, ]
     
125. 
     
126. NO2agg <- aggregate(dat ~ codEst, data=NO2,
     
127.                     FUN = function(x) {
     
128.                         c(mean=signif(mean(x), 3),
     
129.                           median=median(x),
     
130.                           sd=signif(sd(x), 3))
     
131.                         })
     
132. NO2agg <- do.call(cbind, NO2agg)
     
133. NO2agg <- as.data.frame(NO2agg)
     
134. 
     
135. library(maptools)
     
136. ## Link aggregated data with stations to obtain a SpatialPointsDataFrame.
     
137. ## Codigo and codEst are the stations codes
     
138. idxNO2 <- match(airStations$Codigo, NO2agg$codEst)
     
139. NO2sp <- spCbind(airStations[, c('Nombre', 'alt')], NO2agg[idxNO2, ])
     
140. save(NO2sp, file='data/NO2sp.RData')
     
141. 
     
142. ##################################################################
     
143. ## Spanish General Elections
     
144. ##################################################################
     
145. 
     
146. dat2011 <- read.csv('data/GeneralSpanishElections2011.gz')
     
147. 
     
148. census <- dat2011$Total.censo.electoral
     
149. validVotes <- dat2011$Votos.válidos
     
150. ## Election results per political party and municipality
     
151. votesData <- dat2011[, 12:1023]
     
152. ## Abstention as an additional party
     
153. votesData$ABS <- census - validVotes
     
154. ## Winner party at each municipality
     
155. whichMax <- apply(votesData,  1, function(x)names(votesData)[which.max(x)])
     
156. ## Results of the winner party at each municipality
     
157. Max <- apply(votesData, 1, max)
     
158. ## OTH for everything but PP, PSOE and ABS
     
159. whichMax[!(whichMax %in% c('PP',  'PSOE', 'ABS'))] <- 'OTH'
     
160. ## Percentage of votes with the electoral census
     
161. pcMax <- Max/census * 100
     
162. 
     
163. ## Province-Municipality code. sprintf formats a number with leading zeros.
     
164. PROVMUN <- with(dat2011, paste(sprintf('%02d', Código.de.Provincia),
     
165.                                sprintf('%03d', Código.de.Municipio),
     
166.                                sep=""))
     
167. 
     
168. votes2011 <- data.frame(PROVMUN, whichMax, Max, pcMax)
     
169. write.csv(votes2011, 'data/votes2011.csv', row.names=FALSE)
     
170. 
     
171. ##################################################################
     
172. ## CM SAF
     
173. ##################################################################
     
174. 
     
175. library(raster)
     
176. 
     
177. tmp <- tempdir()
     
178. unzip('data/SISmm2008_CMSAF.zip', exdir=tmp)
     
179. filesCMSAF <- dir(tmp, pattern='SISmm')
     
180. SISmm <- stack(paste(tmp, filesCMSAF, sep='/'))
     
181. ## CM-SAF data is average daily irradiance (W/m2). Multiply by 24
     
182. ## hours to obtain daily irradiation (Wh/m2)
     
183. SISmm <- SISmm * 24
     
184. 
     
185. ## Monthly irradiation: each month by the correspondent number of days
     
186. daysMonth <- c(31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
     
187. SISm <- SISmm * daysMonth / 1000 ## kWh/m2
     
188. ## Annual average
     
189. SISav <- sum(SISm)/sum(daysMonth)
     
190. writeRaster(SISav, file='SISav')
     
191. 
     
192. library(raster)
     
193. ## http://neo.sci.gsfc.nasa.gov/Search.html?group=64
     
194. pop <- raster('875430rgb-167772161.0.FLOAT.TIFF')
     
195. ## http://neo.sci.gsfc.nasa.gov/Search.html?group=20
     
196. landClass <- raster('241243rgb-167772161.0.TIFF')

复制代码

1. ##################################################################
   
2. ## Source code for the book: "Displaying time series, spatial and
   
3. ## space-time data with R"
   
4. ##
   
5. ## Copyright (C) 2013-2012 Oscar Perpiñán Lamigueiro
   
6. ##
   
7. ## This program is free software you can redistribute it and/or modify
   
8. ## it under the terms of the GNU General Public License as published
   
9. ## by the Free Software Foundation; either version 2 of the License,
   
10. ## or (at your option) any later version.
    
11. ##
    
12. ## This program is distributed in the hope that it will be useful, but
    
13. ## WITHOUT ANY WARRANTY; without even the implied warranty of
    
14. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    
15. ## General Public License for more details.
    
16. ##
    
17. ## You should have received a copy of the GNU General Public License
    
18. ## along with this program; if not, write to the Free Software
    
19. ## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
    
20. ## 02111-1307, USA.
    
21. ####################################################################
    
22. 
    
23. ##################################################################
    
24. ## Initial configuration
    
25. ##################################################################
    
26. ## Clone or download the repository and set the working directory
    
27. ## with setwd to the folder where the repository is located.
    
28. 
    
29. ##################################################################
    
30. ## OpenStreetMap with Hill Shade layers
    
31. ##################################################################
    
32. 
    
33. ##################################################################
    
34. ## Retrieving data from OpenStreetMap
    
35. ##################################################################
    
36. 
    
37. library('osmar')
    
38. 
    
39. api <- osmsource_api()
    
40. ymax <- 43.7031
    
41. ymin <- 43.6181
    
42. xmax <- -8.0224
    
43. xmin <- -8.0808
    
44. box <- corner_bbox(xmin, ymin, xmax, ymax)
    
45. cedeira <- get_osm(box, source=api, full=TRUE)
    
46. 
    
47. summary(cedeira$nodes)
    
48. 
    
49. idxHighways <- find(cedeira, way(tags(k=='highway')))
    
50. highways <- subset(cedeira, way_ids=idxHighways)
    
51. idxStreets <- find(highways, way(tags(v=='residential')))
    
52. idxPrimary <- find(highways, way(tags(v=='primary')))
    
53. idxSecondary <- find(highways, way(tags(v=='secondary')))
    
54. idxTertiary <- find(highways, way(tags(v=='tertiary')))
    
55. idxOther <- find(highways,
    
56.                  way(tags(v=='unclassified' |
    
57.                           v=='footway' |
    
58.                           v=='steps')))
    
59. 
    
60. spFromOSM <- function(source, index, type='lines'){
    
61.   idx <- find_down(source, index)
    
62.   obj <- subset(source, ids=idx)
    
63.   objSP <- as_sp(obj, type)
    
64.   }
    
65. 
    
66. streets <- spFromOSM(cedeira, way(idxStreets))
    
67. primary <- spFromOSM(cedeira, way(idxPrimary))
    
68. secondary <- spFromOSM(cedeira, way(idxSecondary))
    
69. tertiary <- spFromOSM(cedeira, way(idxTertiary))
    
70. other <- spFromOSM(cedeira, way(idxOther))
    
71. 
    
72. idxPlaces <- find(cedeira, node(tags(k=='name')))
    
73. places <- spFromOSM(cedeira, node(idxPlaces), 'points')
    
74. 
    
75. nms <- subset(cedeira$nodes$tags, subset=(k=='name'), select=c('id', 'v'))
    
76. ord <- match(idxPlaces, nms$id)
    
77. nms <- nms[ord,]
    
78. places$name <- nms$v[ord]
    
79. 
    
80. ## Cedeira town will be printed differently
    
81. idxCedeira <- which(nms$v=='Cedeira') ##Main town
    
82. cedeiraCoords <- coordinates(places[idxCedeira,])
    
83. places <- places[-idxCedeira,]
    
84. 
    
85. ##################################################################
    
86. ## Hill Shading
    
87. ##################################################################
    
88. 
    
89. library(raster)
    
90. ## Galicia DEM
    
91. ## http://ide.unex.es/geonetwork/srv/es/main.search?any=MDE_Galicia
    
92. ## http://ide.unex.es:8180/geonetwork/srv/es/resources.get?id=21&fname=dem_gal.7z&access=private
    
93. 
    
94. old <- tempdir()
    
95. download.file('http://ide.unex.es:8180/geonetwork/srv/es/resources.get?id=21&fname=dem_gal.7z&access=private', 'dem_gal.7z')
    
96. unzip('dem_gal.7z')
    
97. demGalicia <- raster('dem_gal.asc')
    
98. setwd(old)
    
99. 
    
100. cedeiraSP <- as_sp(cedeira, 'points')
     
101. projCedeira <- projection(cedeiraSP)
     
102. ##extCedeira <- bbox(cedeiraSP)
     
103. ## or summary(cedeira$nodes)$bbox
     
104. extCedeira <- extent(-8.15, -7.95, 43.6, 43.75)
     
105. demCedeira <- crop(demGalicia, extCedeira)
     
106. projection(demCedeira) <- projCedeira
     
107. demCedeira[demCedeira <= 0] <- NA
     
108. 
     
109. slope <- terrain(demCedeira, 'slope')
     
110. aspect <- terrain(demCedeira, 'aspect')
     
111. hsCedeira <- hillShade(slope=slope, aspect=aspect,
     
112.                        angle=20, direction=30)
     
113. 
     
114. ##################################################################
     
115. ## Overlaying layers of information
     
116. ##################################################################
     
117. 
     
118. pdf(file="figs/cedeiraOsmar.pdf")
     
119. library(maptools)
     
120. library(latticeExtra)
     
121. library(colorspace)
     
122. library(rasterVis)
     
123. 
     
124. ##Auxiliary function to display the roads. A thicker black line in
     
125. ##the background and a thinner one with an appropiate color.
     
126. sp.road <- function(line, lwd=5, blwd=7,
     
127.                     col='indianred1', bcol='black'){
     
128.   sp.lines(line, lwd=blwd, col=bcol)
     
129.   sp.lines(line, lwd=lwd, col=col)
     
130. }
     
131. 
     
132. ## The background color of the panel is set to blue to represent the sea
     
133. hsTheme <- modifyList(GrTheme(), list(panel.background=list(col='skyblue3')))
     
134. ## DEM with terrain colors and semitransparency
     
135. terrainTheme <- modifyList(rasterTheme(region=terrain_hcl(n=15)),
     
136.                                 list(regions=list(alpha=0.6)))
     
137. ## Hill shade and DEM overlaid
     
138. levelplot(hsCedeira, maxpixels=ncell(hsCedeira),
     
139.           par.settings=hsTheme, margin=FALSE, colorkey=FALSE) +
     
140.   levelplot(demCedeira, maxpixels=ncell(demCedeira),
     
141.             par.settings=terrainTheme) +
     
142.   ## Roads and places
     
143.   layer({
     
144.     ## Street and roads
     
145.     sp.road(streets, lwd=1, blwd=2, col='white')
     
146.     sp.road(other, lwd=2, blwd=3, col='white')
     
147.     sp.road(tertiary, lwd=3, blwd=4, col='palegreen')
     
148.     sp.road(secondary, lwd=4, blwd=6, col='midnightblue')
     
149.     sp.road(primary, col='indianred1')
     
150.     ## Places except Cedeira town
     
151.     sp.points(places, pch=19, col='black', cex=0.4, alpha=0.8)
     
152.     sp.pointLabel(places, labels=places$name,
     
153.                       fontfamily = 'Palatino',
     
154.                       cex=0.6, col='black')
     
155.     ## Cedeira town
     
156.     panel.points(cedeiraCoords, pch=18, col='black', cex=1)
     
157.     panel.text(cedeiraCoords, labels='Cedeira', pos=2, offset=1)
     
158.     })
     
159. dev.off()

复制代码

1. ##################################################################
   
2. ## Source code for the book: "Displaying time series, spatial and
   
3. ## space-time data with R"
   
4. 
   
5. ## Copyright (C) 2012 Oscar Perpiñán Lamigueiro
   
6. 
   
7. ## This program is free software you can redistribute it and/or modify
   
8. ## it under the terms of the GNU General Public License as published
   
9. ## by the Free Software Foundation; either version 2 of the License,
   
10. ## or (at your option) any later version.
    
11. 
    
12. ## This program is distributed in the hope that it will be useful, but
    
13. ## WITHOUT ANY WARRANTY; without even the implied warranty of
    
14. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    
15. ## General Public License for more details.
    
16. 
    
17. ## You should have received a copy of the GNU General Public License
    
18. ## along with this program; if not, write to the Free Software
    
19. ## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
    
20. ## 02111-1307, USA.
    
21. ####################################################################
    
22. 
    
23. ##################################################################
    
24. ## SIAR
    
25. ##################################################################
    
26. 
    
27. ##################################################################
    
28. ## Daily data of different meteorological variables
    
29. ##################################################################
    
30. 
    
31. library(zoo)
    
32. 
    
33. aranjuez <- read.zoo("data/aranjuez.gz",
    
34.                      index.column = 3, format = "%d/%m/%Y",
    
35.                      fileEncoding = 'UTF-16LE',
    
36.                      header = TRUE, fill = TRUE,
    
37.                      sep = ';', dec = ",", as.is = TRUE)
    
38. aranjuez <- aranjuez[, -c(1:4)]
    
39. 
    
40. names(aranjuez) <- c('TempAvg', 'TempMax', 'TempMin',
    
41.                      'HumidAvg', 'HumidMax',
    
42.                      'WindAvg', 'WindMax',
    
43.                      'Radiation', 'Rain', 'ET')
    
44. 
    
45. 
    
46. summary(aranjuez)
    
47. 
    
48. aranjuezClean <- within(as.data.frame(aranjuez),{
    
49.   TempMin[TempMin>40] <- NA
    
50.   HumidMax[HumidMax>100] <- NA
    
51.   WindAvg[WindAvg>10] <- NA
    
52.   WindMax[WindMax>10] <- NA
    
53. })
    
54. 
    
55. aranjuez <- zoo(aranjuezClean, index(aranjuez))
    
56. 
    
57. save(aranjuez, file='data/aranjuez.RData')
    
58. 
    
59. ##################################################################
    
60. ## Solar radiation measurements from different locations
    
61. ##################################################################
    
62. 
    
63. ##################################################################
    
64. ## Unemployment in the United States
    
65. ##################################################################
    
66. 
    
67. unemployUSA <- read.csv('data/unemployUSA.csv')
    
68. nms <- unemployUSA$Series.ID
    
69. ##columns of annual summaries
    
70. annualCols <- 14 + 13*(0:12)
    
71. ## Transpose. Remove annual summaries
    
72. unemployUSA <- as.data.frame(t(unemployUSA[,-c(1, annualCols)]))
    
73. ## First 7 characters can be suppressed
    
74. names(unemployUSA) <- substring(nms, 7)
    
75. head(unemployUSA)
    
76. 
    
77. library(zoo)
    
78. 
    
79. Sys.setlocale("LC_TIME", 'C')
    
80. idx <- as.yearmon(row.names(unemployUSA), format='%b.%Y')
    
81. unemployUSA <- zoo(unemployUSA, idx)
    
82. 
    
83. isNA <- apply(is.na(unemployUSA), 1, any)
    
84. unemployUSA <- unemployUSA[!isNA,]
    
85. 
    
86. save(unemployUSA, file='data/unemployUSA.RData')
    
87. 
    
88. ##################################################################
    
89. ## Gross National Income and $CO_2$ emissions
    
90. ##################################################################
    
91. 
    
92. library(WDI)
    
93.   
    
94. CO2data <- WDI(indicator=c('EN.ATM.CO2E.PC', 'EN.ATM.CO2E.PP.GD',
    
95.               'NY.GNP.MKTP.PP.CD', 'NY.GNP.PCAP.PP.CD'),
    
96.           start=2000, end=2011,
    
97.           country=c('BR', 'CN', 'DE', 'ES',
    
98.               'FI', 'FR', 'GR', 'IN', 'NO', 'US'))
    
99. 
    
100. names(CO2data) <- c('iso2c', 'Country.Name', 'Year',
     
101.                     'CO2.capita', 'CO2.PPP',
     
102.                     'GNI.PPP', 'GNI.capita')
     
103. 
     
104. isNA <- apply(is.na(CO2data), 1, any)
     
105. CO2data <- CO2data[!isNA, ]
     
106. 
     
107. CO2data$Country.Name <- factor(CO2data$Country.Name)
     
108. 
     
109. save(CO2data, file='data/CO2.RData')

复制代码

谢谢楼主的介绍。

ReneeBK 发表于 2014-11-21 04:08
Displaying Time Series, Spatial, and Space-Time Data with R
https://github.com/oscarperpinan/spac ...

лл

楼主，怎么只有数据和代码，怎么没有对应的PDF书籍Displaying time series, spatial and space-time data with R？

楼主发的两本都是只有代码没有书，要是有书就好了