R script to create and publish interpolated temperature map of Israel

IMS_Map.R

###############################################################
# Code to publish interpolated temperature map                #
# based on current observations from the                      #
# Israel Meteorological Service (IMS) website                 #
# Author: Michael Dorman                                      #
###############################################################

library(XML)
library(httr)
library(plyr)
library(rgdal)
library(sp)
library(raster)
library(rgeos)
library(gstat)
library(automap)
library(plotGoogleMaps)
library(markdown)

# Fixed values and paths
url.name = "http://www.ims.gov.il/ims/PublicXML/observ.xml"

# Read XML contents
url.get = GET(url.name)
url.content = content(url.get, as = "text")
dat = xmlToList(url.content)

# Subset most recent observation (top of XML file)
time_positions = which(names(dat) == "date_selected")
time_raw = dat[time_positions][1]

# Process time data
time_raw = time_raw[[1]]
time_raw = gsub(" ", "", time_raw, fixed = TRUE)
time = as.POSIXct(time_raw, format = "%Y%m%d%H", tz = "GMT")
time = format(time, tz = "Asia/Jerusalem", usetz = TRUE)
time = as.POSIXct(time)
time_char = format(time, format = "%Y%m%d%H")

# Convert to table and pre-process
dat = dat[(time_positions[1]+1):(time_positions[2]-1)]
dat = ldply(dat, data.frame, stringsAsFactors = FALSE)
dat = cbind(dat[dat$.id == "surface_station", ],
            dat[dat$.id == "surface_observation", ])
dat[, !(colnames(dat) %in% c(".id", "station_name"))] = 
  colwise(as.numeric)(dat[, !(colnames(dat) %in% c(".id", "station_name"))])
dat = dat[, which(colnames(dat) != ".id")]
dat = dat[, !apply(dat, 2, function(x) all(is.na(x)))]
colnames(dat) = sapply(strsplit(colnames(dat), ".", fixed = TRUE), "[[", 1)
colnames(dat)[which(colnames(dat) == "station_height")] = "elevation"
dat$time = time

# Convert to point layer
coordinates(dat) = ~ station_lon + station_lat
proj4string(dat) = "+proj=longlat +datum=WGS84"

# Read elevation raster
grid = raster("C:\\Users\\Michael Dorman\\Google Drive\\R Code\\XML\\SRTM_Israel_900.tif")
dat = spTransform(dat, CRS(proj4string(grid)))

# Remove empty observations
dat = dat[!is.na(dat$TT), ]

# Apparently TT data are sometimes given in deg. and sometimes in 1/10 deg...
# This is a quick solution -
if(max(dat$TT) > 50) dat$TT = dat$TT/ 10

# Perform Universal Kriging interpolation
v = autofitVariogram(TT ~ elevation, dat)
g = gstat(
  formula = TT ~ elevation, 
  model = v$var_model, 
  data = dat)
names(grid) = "elevation"
z = interpolate(grid, g, xyOnly = FALSE)
z = mask(z, grid)

z = round(z)

# Convert raster to Polygons
pol = rasterToPolygons(z, n = 4, dissolve = TRUE)

# Prepare vector layers
pnt = dat[, c("station_name", "time", "TT")]
pol$observations = nrow(pnt)
pol$time = time
colnames(pol@data)[which(colnames(pol@data) == "layer")] = "TT"

# Create fixed color scale for 1-35 deg. range
cols = rev(rainbow(40)[1:35])
cols = cols[pol$TT+5]

# Rename variable
colnames(pol@data)[which(colnames(pol@data) == "TT")] = "T (deg. C)"
colnames(dat@data)[which(colnames(dat@data) == "TT")] = "T (deg. C)"

# Layer name
l_name = as.character(time)

# Create base map of interpolated values
m_pol = plotGoogleMaps(pol,
                       zcol = "T (deg. C)",
                       add = TRUE,
                       colPalette = cols,    		
                       strokeWeight = 0.1,
                       layerName = l_name,
                       mapTypeId = "SATELLITE",
                       legend = FALSE,
                       openMap = FALSE)

# Add observed values and save HTML
m_pnt = plotGoogleMaps(dat[, c("station_name", "T (deg. C)")],
                       zcol = "T (deg. C)",
                       previousMap = m_pol,
                       filename = "C:\\Test\\T_Israel.html",
                       iconMarker = iconlabels(dat$TT, height=10, icon = TRUE), 
                       strokeWeight = 0.01,
                       layerName = paste0("Stations (", nrow(dat), ")"),
                       legend = FALSE,
                       openMap = FALSE)

# Publish HTML to Rpubs
rpubsUpload("Current Temperature in Israel", "C:\\Test\\T_Israel.html",
	id = "https://api.rpubs.com/api/YOUR_RPUBS_DOCUMENT_ID_HERE")