Archivio tag: r

Monitorare le acque del Golfo del Messico con R

4 anni fa

In seguito all’incidente alla piattaforma petrolifera Deepwater Horizon, l’EPA ha avviato un programma di monitoraggio straordinario delle acque del Golfo del Messico, in particolare presso le coste di Louisiana, Mississippi, Alabama e Florida, per valutare la presenza di petrolio e altri rifuiti oleosi (sezione “Waste Management“), sostanze inquinanti associate sia al petrolio che ai prodotti utilizzati come dispersanti (sezione “Water data“), inquinanti in aria (sezione Air data) e sedimenti (sezione Sediments data).
Nei dataset dei risultati delle analisi che vengono pubblicati periodicamente sul sito sono incluse anche le coordinate di ciascun punto di monitoraggio, in modo da poterne visualizzarne le variazioni su una mappa della zona.

Per chi volesse farlo con R, è possibile utilizzare le mappe messe a disposizione dalla libreria “maps” seguendo questo post su R-Chart.

Io ho preferito usare il pacchetto “RgoogleMaps“, che fornisce una comoda interfaccia per recuperare mappe statiche da Google ed usarle come sfondo per altri grafici realizzati con R.
I dati che ho utilizzato sono quelli relativi alle analisi sulle acque, nelle quali nessuna delle sostanze testate è stata fino ad ora rilevata in concentrazioni significative.
Il grafico indica quindi solo le posizioni dei punti di monitoraggio (il fine del post è di mostrare l’utilità delle funzioni “maps” e “RgoogleMaps”), ma nel caso in cui si desideri visualizzare anche la presenza degli inquinanti è sufficiente sostituire le variabili relative alle coordinate con quelle relative alle sostanze (nel formato “data$variabile”).

library(RgoogleMaps)
data=read.table('http://www.epa.gov/bpspill/data/water_sampling_update.csv', header=TRUE, skip=1,sep=",")
myMap <- GetMap(markers="28.738889,-88.589722,black", center=c(29.108889,-89.289722), zoom = 8, destfile = "myMap.png", maptype = "mobile")
PlotOnStaticMap(myMap, lon=data$LONGITUDE, lat=data$LATITUDE, pch=21, col='grey', bg="red", verbose=0, add=F, cex=0.8)

Il marker nero, che indica la posizione della piattaforma, è aggiunto selezionando latitudine e longitudine del punto mediante l’attributo “markers”.
Il tipo di mappe è definito dall’attributo “maptype”.
Usando come nel codice di esempio l’opzione “mobile” (maptype = "mobile"), il grafico che si ottiene è:

Scegliendo invece l’opzione “satellite” (maptype = "satellite"):

Ricordo che per utilizzare il pacchetto “RgoogleMaps” e la “Google Maps API” è necessario richiedere una “Maps API key“.

L’elefante e la farfalla – interpretata da Excel e R

4 anni fa

Prendendo spunto da questo post di Gaspar in cui presenta una macro da lui programmata in Excel per risolvere un problema di ordinamento di alcuni ordini di un’azienda, sono rimasto sconcertato dal numero spropositato di comandi necessario per effettuare un’operazione così “semplice”.

Partendo da una tabella di ordini come questa,

è stato richiesto agli impiegati dell’azienda di ristrutturarla e di ottenere una lista come questa (perchè?).

Ho pensato a come si potesse fare in R nel minor numero di righe di codice possibile, e dopo una decina di minuti di prove questo è il risultato:

library(reshape)
ordini <- read.csv("ordini.csv")
ordini2 <- melt(ordini, id="Prodotto")
names(ordini2) <- c("Prodotti","Taglia","Quantità")
ordini3 <- ordini2[order(ordini2$Prodotti), ]

Lungi dal voler appiccare un flame sul confronto Excel vs R, i cui campi di applicazione e il target d’utenza non sempre sono i medesimi per la natura stessa di questi software, sono convinto che R potrebbe essere in diverse situazioni molto più utile di Excel anche in ambito aziendale.
E paradossalmente anche più semplice da utilizzare, dopo averne scalato un po’ la (più ripida) curva d’apprendimento.

Nota: le immagini le ho realizzate utilizzando il pratico package “xtable”, che ha convertito l’output del codice di R in un documento in LaTeX, riconvertito poi da me in pdf.
Il tutto con sole 3 linee in più di codice.

Il Metodo D’Hondt e la salmo trutta

5 anni fa

Leggendo il libretto d’istruzioni per le operazioni di elezione dei consiglieri comunali ho scoperto che per l’attribuzione dei seggi alle diverse liste si usa una procedura chiamata “Metodo D’Hondt“.
Finora sia come presidente di seggio che come scrutatore non mi ero mai posto il problema di come venissero ripartiti i seggi tra le liste, ma visto che domani farò parte dell’ufficio elettorale incaricato di farlo mi sono documentato un po’ sull’applicazione del metodo.

Il “Metodo D’Hondt” prevede che il totale dei voti assegnati a ciascuna lista venga diviso per 1, 2, 3, . , n (dove n è il numero totale di seggi da assegnare).
I risultati di ciascuna divisione vengono inseriti in una tabella le cui righe rappresentano gli n numeri interi per cui il totale di lista è stato diviso, e le colonne le diverse liste a cui i seggi devono essere assegnati.
Si selezionano poi gli n migliori risultati, e si contano quanti di essi sono corrispondenti a ciascuna lista.
Questo numero corrisponderà al numero di seggi assegnato a quella lista.

Si capisce meglio con un esempio, spiegato con una funzione in R trovata qui.

Alle elezioni del Consiglio delle Trote, composto da 10 seggi, si presentano cinque partiti, con i seguenti risultati:

Voti validi: 28.500 schede
Partito Salmo ferox: 10.000 preferenze
Partito Salmo nigripinnis: 7.800 preferenze
Partito Salmo stomachicus: 5.800 preferenze
Partito Salmo marmoratus: 3.000 preferenze
Partito Salmo ezenami: 1.900 preferenze

I dati risultanti dalle divisioni sono:

Numeri divisori Salmo ferox Salmo nigripinnis Salmo stomachicus Salmo marmoratus Salmo ezenami
1 10.000 7.800 5.800 3.000 1.900
2 5.000 3.900 2.900 1.500 950
3 3.333,3 2.600 1.933,3 1000 633,3
4 2.500 1950 1450 750 475
5 2000 1560 1160 600 380

I numeri in grassetto corrispondono ai risultati migliori forniti dalla funzione:

dHont <- function( candidates, votes, seats ){
    tmp <- data.frame(
                candidates = rep( candidates, each = seats ),
                scores     = as.vector(sapply( votes, function(x) x /
1:seats ))
                )
    tmp <- tmp$candidates[order( - tmp$scores )] [1:seats]
    table(tmp)
}

La funzione viene applicata così:

> votes <- c(10000,7800,5800,3000,1900)
> dHont(letters[1:5], votes, 10)
tmp
a b c d e 
4 3 2 1 0

Il Consiglio delle Trote sarà quindi costituito da:

4 consiglieri del Partito Salmo ferox
3 consiglieri del Partito Salmo nigripinnis
2 consiglieri del Partito Salmo stomachicus
1 consiglieri del Partito Salmo marmoratus
0 consiglieri del Partito Salmo ezenami

Monitorare il clima con R

5 anni fa

Il sito processtrends.com presenta degli utili script in R per generare dei grafici relativi ad alcuni degli indicatori più usati nel monitoraggio dello stato del clima (sono quelli più legati al riscaldamento globale).
I dati utilizzati sono sempre i più recenti disponibili, quindi il codice va modificato solo nel caso si desideri specificare intervalli temporali diversi o includere altre variabili oltre a quelle di default.

Gli indicatori disponibili (con relativi script) sono:

  • aumento della temperatura globale, con i dati GISS, UAHRSS, Hadley e NOAA sulle anomalie delle temperature superficiali
  • aumento del livello medio del mare
  • diminuzione dell’estensione della calotta artica
  • andamento del livello di CO2 in atmosfera (Mauna Loa)
  • irradianza solare totale
  • fenomeni El Niño e La Niña
  • “amplificazione polare”

Nel caso si preferisca usare un altro software (anche Excel, se proprio dovete) sono a disposizione dei file di testo con i dati più recenti.

Frattali con R

5 anni fa

Tra le tante cose che si possono fare con R non potevano mancare i frattali.
Limitandosi all’insieme di Mandelbrot, gli esempi che si trovano in rete sono tanto semplici quanto efficaci.

Il primo, che permette di ottenere sia immagini statiche che gif animate, l’ho trovato grazie alla mailing list R-help.

# written by Jarek Tuszynski, PhD.
# http://tolstoy.newcastle.edu.au/R/help/05/10/13198.html

library(fields)  # for tim.colors
library(caTools) # for write.gif
m = 400          # grid size
C = complex( real=rep(seq(-1.8,0.6, length.out=m), each=m ), imag=rep(seq(-1.2,1.2, length.out=m), m))
C = matrix(C,m,m)

Z = 0
X = array(0, c(m,m,20))
for (k in 1:20) {
  Z = Z^2 C
  X[,,k] = exp(-abs(Z))
} 

image(X[,,k], col=tim.colors(256)) # show final image in R
write.gif(X, "Mandelbrot.gif", col=tim.colors(256), delay=100) # drop "Mandelbrot.gif" file from current directory on any web browser to see the animation

mandelbrot-1

Il secondo esempio, ben più spettacolare (soprattutto se la dimensione delle immagini è grande), è quello di una funzione creata da Mario dos Reis.
Tutti i file necessari sono raccolti nell’archivio “mandelbrot.zip”, scaricabile dalla sua pagina dedicata ai frattali.

Una volta scompattato il file zip, al suo interno troverete il file “mandelc.R”, contenente il codice della funzione mandelbrot().
In essa viene richiamata una routine in C (contenuta nel file “mandelbrot.c”) per rendere più veloce la computazione.

In un terminale, compilare il file “mandelbrot.c” con il comando:

R CMD SHLIB mandelbrot.c

che creerà il file “mandelbrot.so”.
Aprire una sessione di R e dare:

source("mandelc.R")

Caricare l’oggetto “mandelbrot.so”:

dyn.load("mandelbrot.so")

Per una prima prova:

# graph A
image(mandelbrot())

Nel file “examples.R” ci sono due funzioni, “surf.colors” e “zoom”.

Dopo averle copiate ed incollate nel terminale dove avete aperto R (o nella console di R), le istruzioni per ottenere esempi di grafici più complessi sono:

frac <- mandelbrot(iter = 50)
fcol <- c(heat.colors(49), "black")

# graph B
image(frac, col = fcol)

mandelbrot-a-b

Il grafico successivo è in 3D:

# graph C
fcol <- surf.colors(frac, col = c(gray(seq(.5, 1, len = 49)), "black"))
frac$z <- -1/frac$z # nicer for plotting and coloring
persp(frac, col = fcol, border = NA, theta = 225, phi = 45, shade = .5)

mandelbrot-c

Modificando la funzione evidenziata in rosso (con logaritmi o potenze, ad esempio), è possibile ottenere forme sempre diverse:

# graph D

frac <- mandelbrot(x = c(-0.8438146, -0.8226294),
                   y = c(0.1963144, 0.2174996), iter = 500)

image(frac, col = fcol.gr)
image(-1/frac$z, col = fcol.gr)

mandelbrot-d

Per provare la funzione zoom, usate il codice seguente ed all’interno dell’immagine cliccare agli estremi opposti dell’area da ingrandire:

par(mfrow = c(2, 2))

gr <- gray(seq(1, 0, len = 255))
fcol.gr <- c(rep(c(gr, rev(gr)), 2), "black")

image(frac, col = fcol.gr)
zoom(col = fcol.gr) # Now click with your mouse two opposite corners of the
                    # area you wish to enlarge
zoom(col = fcol.gr) # do it again!
zoom(col = fcol.gr) # and again . 

Per chi volesse cimentarsi anche con un’insieme di Julia o un attrattore di Rössler, in questo blog (purtroppo non più aggiornato) si può trovare molto altro codice con cui divertirsi (sempre che non abbiate niente di meglio da fare :-)

Nota: tutto questo funziona in Linux, in Windows e Mac non so. Oltre ad R e ai packages aggiuntivi richiesti, è necessario un compilatore C.

Esempio di infomappa con R

5 anni fa

grafico-R-disoccupazione

Dopo questo articolo su FlowingData in cui si mostra come realizzare una carta tematica (nello specifico una choropleth map, ovvero una mappa colorata con tonalità diverse in base al valore assunto da una certa variabile) sullo stato della disoccupazione negli USA usando Python, alcuni utenti R hanno raccolto la sfida lanciata da David Smith (l’autore dell’ottimo blog Revolutions) e ne hanno riprodotto il risultato in R.

I risultati, ottenuti con approcci diversi, li trovate elencati qui.

C’è poi chi ha utilizzato dati relativi alla situazione occupazionale tedesca e, nel mio caso, italiana.

I dati più recenti (anno 2007) relativi all’occupazione li ho raccolti dal sito dell’Istat, mentre quelli spaziali dal sito GADM (direttamente come oggetti R!), plottati poi con la funzione spplot.
Per chi volesse curiosare nei dati, fate attenzione che la Puglia è chiamata Apulia (almeno fino a quando correggeranno l’errore), cosa che può sconvolgere le corrispondenze tra vettori, notoriamente abbastanza irascibili.

Il mio codice lo trovate sotto.

Continua a leggere