rm(list=ls())
gc(verbose=TRUE, reset=TRUE)

library(stringi)
library(plyr)
library(dplyr)
library(ggplot2)

# Function to get spots for a given month.  Downloads zip file from WSPR web site and saves it as local
# RData if necessary, otherwise reads local RData file.  Performs minor cleanup on downloaded file.

GetSpots <- function(year, month) {
    
    if ( ! dir.exists("./spots") ) {
        dir.create("./spots")
    }
    month <- sprintf("%02d", month)
    rdata <- stri_join("./spots/wsprspots-", year, "-", month, ".RData")
    if ( file.exists(rdata) ) {
        message(stri_join("Loading: ", rdata))
        load(rdata)
    } else {
        url <- stri_join("http://wsprnet.org/archive/wsprspots-", year, "-", month, ".csv.zip")
        zip <- stri_join("./spots/wsprspots-", year, "-", month, ".csv.zip")
        csv <- stri_join("./spots/wsprspots-", year, "-", month, ".csv")
        if ( ! file.exists(csv) ) {
            message(message(stri_join("Downloading: ", zip)))
            download.file(url, zip, method = "curl")
            unzip(zip, exdir = "./spots")
            file.remove(zip)
        }
        message(stri_join("Reading: ", csv))
        spots <- read.csv(csv,
                          header = FALSE,
                          stringsAsFactors = FALSE,
                          col.names = c("id", "ts", "dst.call", "dst.grid", "snr", "freq", "src.call",
                                        "src.grid", "pwr", "drift", "km", "az", "band", "ver", "code"),
                          colClasses = c("integer", "integer", "character", "character", "integer", "numeric",
                                         "character", "character", "integer", "integer", "integer", "integer",
                                         "integer", "character", "integer"))
        # Generate POSIXlt from decimal.
        spots$ts <- as.POSIXlt(spots$ts, origin = "1970-01-01", tz = "GMT")
        # UCASE call signs and grids.
        spots$src.call <- toupper(spots$src.call)
        spots$src.grid <- toupper(spots$src.grid)
        spots$dst.call <- toupper(spots$dst.call)
        spots$dst.grid <- toupper(spots$dst.grid)
        # Dump zero-distance spots.
        spots <- spots[spots$km > 0,]
        # Save as RData for future quick load.
        save(spots, file = rdata)
        file.remove(csv)
    }
    spots
}

bandname <- list()
bandname[[28]] <- "10M"
bandname[[14]] <- "20M"
bandname[[7]]  <- "40M"
rdata <- "./spotinfo.RData"
if ( file.exists(rdata) ) {
    message(stri_join("Loading: ", rdata))
    load(rdata)
} else {
    spotinfo <- NULL
    year <- 2014
    for ( month in c(7, 1) ) {
        spots <- GetSpots(year, month)
        for ( band in c(28, 14, 7) ) {
            # Trim to band of interest and remove ground wave (~100 miles/160km) spots.
            tmp <- spots[spots$band == band & spots$km > 160, c("ts", "src.call", "km")]
            # We will soon rely on the fact that all spots are transmitted on the top of the
            # minute - i.e. at YYYY-MM-MM HH:MM:00. (Note 00 seconds.)  Verify.
            stopifnot(0 == sum(spots$ts$sec != 0))
            # Convert timestamp to character as dplyr can't sort on POSIXlt.
            tmp$tsc <- as.character(tmp$ts)
            tmp$ts <- NULL
            # Group by src.call and timestamp.  I.e. group all receptions of a single transmission.
            g <- group_by(tmp, src.call, tsc)
            # Keep only the most distant reception for each transmission.
            s <- summarise(g, km.max = max(km, na.rm = TRUE))
            # Build data frame such that we can use ggplot2 facets.
            spotinfo <- rbind(spotinfo, data.frame(Month = month.name[month],
                                                   Band = bandname[[band]],
                                                   vals = s$km.max))
        }
    }
    save(spotinfo, file = rdata)
}

g <- ggplot(spotinfo, aes(x = vals)) +
    geom_density(data = spotinfo,
                 alpha = 0.3,
                 aes(x = vals, group = Month, colour = Month)) +
    facet_grid(Band ~ .) + 
    scale_color_brewer(palette="Set1") +
    # Title not required, use caption in Wikipedia to describe.
    # ggtitle("Density Distribution of WSPR Spots\nJanuary 2014 vs July 2014\n(using only most distant reception per spot)") +
    xlab("Kilometers") +
    ylab("Density") +
    theme(axis.title.x = element_text(size = rel(1.25))) +
    theme(axis.title.y = element_text(size = rel(1.25))) +
    theme(strip.text.y = element_text(angle = 0, face = "bold", size = rel(1.5)))
g

svg("./WSPR 2014.svg", width = 6, height = 5)
g
dev.off()

Io, detentore del copyright su quest'opera, dichiaro di pubblicarla con la seguente licenza:

Questo file è disponibile in base alla licenza Creative Commons Attribuzione-Condividi allo stesso modo 4.0 Internazionale

Tu sei libero:

• di condividere – di copiare, distribuire e trasmettere quest'opera
• di modificare – di adattare l'opera

Alle seguenti condizioni:

• attribuzione – Devi fornire i crediti appropriati, un collegamento alla licenza e indicare se sono state apportate modifiche. Puoi farlo in qualsiasi modo ragionevole, ma non in alcun modo che suggerisca che il licenziante approvi te o il tuo uso.
• condividi allo stesso modo – Se remixi, trasformi o sviluppi il materiale, devi distribuire i tuoi contributi in base alla stessa licenza o compatibile all'originale.

https://creativecommons.org/licenses/by-sa/4.0CC BY-SA 4.0 Creative Commons Attribution-Share Alike 4.0 truetrue

 

Questo grafico in grafica vettoriale è stato creato con R.