By: Júlio Hoffimann

Re-posted from: https://juliohm.github.io/science/california-earthquake/

January 4, 2018. It was late night at Stanford when an Earthquake gently woke me up.
For a moment, it felt like a dream…

In this post, I would like to share how I used open data and
GeoStats.jl to estimate Earthquake
magnitude at my home address.

Data

The data provided by USGS
comprises measurements from 519 seismic stations in California. Besides the measurements, USGS
provided estimates of Earthquake magnitude and uncertainty as raster images:

Their estimation procedure, which is explained in the
metadata,
was designed to preserve ground motion information near monitoring stations to the best extent possible.
Locations that were far away from the epicenter at Berkeley were penalized, and as a result,
many “bubble” artifacts were printed on the map.

I wasn’t particularly happy with the visualization, and so I decided to create my own map using
open source software. :ok_hand:

Preprocessing

Raw measurements from all 519 stations were extracted from the KML files, and the data was saved into a single CSV file:

$ head stations.csv
longitude,latitude,magnitude,station
-116.711900,33.714100,1.0,"Keenwild Fire Station, Mountain Center, CA, USA"
-116.459400,33.611700,1.0,"Pinyon Flats Observatory, CA, USA"
-116.847800,33.630000,1.0,"Red Mountain, Riverside Co, CA, USA"
-122.952698,38.823601,1.0,"Hale Ranch"
-122.702583,38.775879,1.0,"Engles Strong Motion"
-122.235580,37.876220,2.9,"Byerly Seismographic Vault, Berkeley, CA, USA"
-122.243180,37.874920,4.1,"LBL Building 67, Berkeley, CA, USA"
-122.254290,37.877150,4.2,"LBL Building 88, Berkeley, CA, USA"
-120.386510,38.034550,1.0,"Columbia College, Columbia, CA, USA"

Problem definition

The data was loaded into a GeoDataFrame:

using GeoStats

geodata = readtable("stations.csv", coordnames=[:longitude,:latitude])

# 519×4 GeoDataFrame (longitude and latitude)
# longitude	latitude	magnitude	station
# 1	-116.712	33.7141	1.0	Keenwild Fire Station, Mountain Center, CA, USA
# 2	-116.459	33.6117	1.0	Pinyon Flats Observatory, CA, USA
# 3	-116.848	33.63	1.0	Red Mountain, Riverside Co, CA, USA
# 4	-122.953	38.8236	1.0	Hale Ranch
# 5	-122.703	38.7759	1.0	Engles Strong Motion
# ...

and the estimation problem was defined for the area of interest:

minlon, maxlon = -123.5068, -121.0068
minlat, maxlat = 37.0219, 38.6885

domain  = RegularGrid([300,200], [minlon,minlat], [.0083,.0083])
problem = EstimationProblem(geodata, domain, :magnitude)

# 2D EstimationProblem
#   data:      519×4 GeoDataFrame (longitude and latitude)
#   domain:    300×200 RegularGrid{Float64,2}
#   variables: magnitude (Float64)

Solution

The problem was solved using Universal Kriging with a spherical variogram model.
The model is simplistic for this
cyclostationary process,
but it serves as an example:

# variogram model from variogram analysis
γ = SphericalVariogram(sill=2.90, range=2.00, nugget=0.10)

# Universal Kriging with quadratic drift
solver = Kriging(:magnitude => @NT(variogram=γ,degree=2))

solution = solve(problem, solver)

# 2D EstimationSolution
#   domain: 300×200 RegularGrid{Float64,2}
#   variables: magnitude

The solution was then plotted on top of a basemap provided by
Mapbox via the
MapTiles.jl package:

using MapTiles

provider = MapTiles.MapBoxProvider(accesstoken=ENV["MAPBOX_AUTH"], id="dark")
basemap  = MapTiles.fetchmap(minlon, minlat, maxlon, maxlat, provider=provider)

plot(basemap)
heatmap!(solution.mean[:magnitude], alpha=.4)

Great! I like this visualization better. :+1: Now let’s go back to the original question:
what is the magnitude of the Earthquake at a speficic location?

To answer the question, we change the domain from a RegularGrid to a PointCollection
containing a few locations of interest, and solve the problem:

# locations of interest
UCBerkeley = [-122.268093, 37.872051]
Stanford   = [-122.169525, 37.428890]

domain  = PointCollection([UCBerkeley Stanford])
problem = EstimationProblem(geodata, domain, :magnitude)

solution = solve(problem, solver)

solution.mean[:magnitude]

# 2-element Array{Float64,1}:
#  4.62765
#  3.73319

I felt a magnitude of $\approx 3.73$ whereas people in UC Berkeley felt a more intense
Earthquake as expected with a magnitude of $\approx 4.63$.

By: Júlio Hoffimann

Re-posted from: https://juliohm.github.io/dataviz/where-are-the-julians/

Data

The data was extracted from METADATA. It only includes members
of the community that have contributed to a registered Julia package (e.g. issues, pull requests) up
until 27-June-2017.

The Jupyter notebook used for data extraction is available in our
JuliaGraphsTutorials repository.

Say hello to a Julian near you. #JuliansInTheGlobe

By: Júlio Hoffimann

Re-posted from: https://juliohm.github.io/dataviz/DAG-of-Julia-packages/

Instructions

This interactive visualization made with D3 shows the Directed Acyclic Graph (DAG) of all registered Julia packages up until 27-June-2017.

The size of a node represents its influence (i.e. out degree) in the DAG. The color represents the required Julia version.

Hover the mouse over the elements to get more information.

Data

The data was extracted from METADATA with the following script:

UPDATE: An updated version of the script can be found in the JuliaGraphsTutorials repository in the form of a Jupyter notebook. Additional information about the Julia package system is extracted therein with the GitHub.jl package.

using JSON
using LightGraphs
using ProgressMeter

# find all packages in METADATA
pkgs = readdir(Pkg.dir("METADATA"))
filterfunc = p -> isdir(joinpath(Pkg.dir("METADATA"), p)) && p ∉ [".git",".test"]
pkgs = filter(filterfunc, pkgs)

# assign each package an id
pkgdict = Dict{String,Int}()
for (i,pkg) in enumerate(pkgs)
  push!(pkgdict, pkg => i)
end

# build DAG
G = DiGraph(length(pkgs))
@showprogress 1 "Building graph..." for pkg in pkgs
  children = Pkg.dependents(pkg)
  for c in children
    add_edge!(G, pkgdict[pkg], pkgdict[c])
  end
end

# find required Julia version
juliaversions = String[]
for pkg in pkgs
  versiondir = joinpath(Pkg.dir("METADATA"), pkg, "versions")
  if isdir(versiondir)
    latestversion = readdir(versiondir)[end]
    reqfile = joinpath(versiondir, latestversion, "requires")
    juliaversion = string(get(Pkg.Reqs.parse(reqfile), "julia", "NA"))
    push!(juliaversions, juliaversion)
  else
    push!(juliaversions, "BOGUS")
  end
end

# construct JSON
nodes = [Dict("id"=>pkgs[v],
              "indegree"=>indegree(G,v),
              "outdegree"=>outdegree(G,v),
              "juliaversion"=>juliaversions[v]) for v in vertices(G)]
links = [Dict("source"=>pkgs[u], "target"=>pkgs[v]) for (u,v) in edges(G)]
data = Dict("nodes"=>nodes, "links"=>links)

# write to file
open("DAG-Julia-Pkgs.json", "w") do f
  JSON.print(f, data, 2)
end