# Simple Data Science Of Global Warming In KDnuggets

Would love to get a post from you for KDnuggets (Gregory Piatetsky, KDnuggets President)

Some days ago, Gregory Piatetsky invited me to write a post for KDnuggets. I couldn’t say no. He suggested to me some topics and I decided to experiment around climate change to demonstrate how easy is to see some evidences of this alarming menace. You can read the post here.

This is the code I wrote to do this experiment:

```require(sqldf)
require(ggplot2)
require(ggthemes)
setwd("YOUR WORKING DIRECTORY HERE") #This line doen's work until you type a valid path
files=list.files(pattern = "TG_STAID")
results=data.frame(staid=character(0), trnd=numeric(0), nobs=numeric(0))
#Loop to calculate linear trends
for (i in 1:length(files))
{
table=table[table\$Q_TG==0,]
table\$date=as.Date(as.character(table\$DATE), "%Y%m%d")
results=rbind(data.frame(staid=files[i], trnd=lm.fit (x = matrix(table\$date), y = table\$TG)\$coefficients, nobs=nrow(table)), results)
}
write.csv(results, file="results.csv", row.names = FALSE)#Save your work
#Remove outliers
results=results[!results\$trnd %in% boxplot.stats(results\$trnd, coef = 4)\$out,]
#Histogram
hist(x=results\$trnd, breaks = 50,
col = "orange",
border = "pink",
freq=TRUE,
xlim=c(-0.03, 0.03),
ylim=c(0, 300),
xlab="Historical trend of daily mean temperatures",
ylab="Number of stations",
main="Evolution of temperatures in Europe and the Mediterranean",
sub="Source: European Climate Assessment & Dataset project")
results\$staid2=as.numeric(gsub("[^0-9]","",results\$staid)) #To join results with geographical coordinates
#Right tail of the distribution
hotpoints=sqldf("SELECT a.staid, a.staid2, a.trnd, a.nobs, b.staname, b.lat, b.lon
FROM results a INNER JOIN staids b ON (a.staid2=b.staid) WHERE TRND >= 0.02")
#Convert degrees:minutes:seconds to decimal coordinates
hotpoints=within(hotpoints, {
dms=do.call(rbind, strsplit(as.character(LAT), ":"))
lat=sign(as.numeric(dms[,1]))*(abs(as.numeric(dms[,1]))+(as.numeric(dms[,2]) + as.numeric(dms[,3])/60)/60);rm(dms)
})
hotpoints=within(hotpoints, {
dms=do.call(rbind, strsplit(as.character(LON), ":"))
lon=sign(as.numeric(dms[,1]))*(abs(as.numeric(dms[,1]))+(as.numeric(dms[,2]) + as.numeric(dms[,3])/60)/60);rm(dms)
})
#To make readable tha name of the station in the map
hotpoints\$staname=gsub("^\\s+|\\s+\$", "", hotpoints\$STANAME)
#To prepare coordinates to gvis function
hotpoints\$LatLong=with(hotpoints, paste(lat, lon, sep=":"))
#The amazing gvisMap function of googleVis package
hotpointsMap=gvisMap(hotpoints, "LatLong" , "STANAME",
options=list(showTip=TRUE,
showLine=TRUE,
enableScrollWheel=TRUE,
mapType='terrain',
useMapTypeControl=TRUE))
plot(hotpointsMap)
#The plot one of this hot stations