impact crater data plot

iwantpngfiles <- TRUE

waitreturn <- function(prompt="Press Return") {
    if (interactive() && !iwantpngfiles) readline(prompt)
}

#install.packages("RJSONIO")
library(RJSONIO)
j <- fromJSON("Impact_database_2010_1.json") # or wherever you put it
y <- data.frame(matrix(unlist(lapply(j, function(xx) { row <- c(xx$age, xx$diameter); if (length(row) == 1) { row <- c(row, NA) }; return(row) })), ncol=2, byrow=T))

if (iwantpngfiles) {
    png("impact_craters_%01d.png", 500)
}

par(mar=c(5, 4, 1, 1))

plot(y, xlab="Age [Myr]", ylab="Crater diameter [km]", log="xy", pch="★", xlim=c(2e7, 1e-5), ylim=c(0.01, 1e3), axes=FALSE)
grid()
axis(1, at=10^(-5:7), label=c(expression(10^-5, 10^-4, 10^-3, 10^-2, 0.1, 1, 10, 100, 10^3, 10^4, 10^5, 10^6, 10^7)))
axis(2, at=10^(-2:3), label=c(expression(10^-2, 0.1, 1, 10, 100, 10^3)))

waitreturn()

plot(y, xlab="Age [Myr]", ylab="Crater diameter [km]", log="xy", pch="★", xlim=c(3e3, 1e-5), ylim=c(0.01, 1e3), axes=FALSE)
grid()
axis(1, at=10^(-5:3), label=c(expression(10^-5, 10^-4, 10^-3, 10^-2, 0.1, 1, 10, 100, 1000)))
axis(2, at=10^(-2:3), label=c(expression(10^-2, 0.1, 1, 10, 100, 10^3)))

waitreturn()

plot(y, xlab="Age [Myr]", ylab="Crater diameter [km]", log="xy", pch="★", xlim=c(3e3, 1e-5), ylim=c(0.01, 1e3), axes=FALSE)
grid()
axis(1, at=10^(-5:3), label=c(expression(10^-5, 10^-4, 10^-3, 10^-2, 0.1, 1, 10, 100, 1000)))
axis(2, at=10^(-2:3), label=c(expression(10^-2, 0.1, 1, 10, 100, 10^3)))
curve((1e3*x)^0.5, add=TRUE, lty="dashed")
text(0.01, 200, "anthropic shadow?", cex=2)
curve((x/1e3)^0.5, add=TRUE, lty="dashed")
text(200, 0.03, "too boring\nto record?\neroded away?", cex=1.5)

if (iwantpngfiles) {
    dev.off()
}