Message Boards Message Boards

Bárðarbunga earthquake swarm

GROUPS:

Data credit: Icelandic Met Office

Animated results (using techniques shown below) and how to generate them.

In 2010 the volcano Eyjafjallajökull erupted and shut down (among other things) transatlantic air travel for about six days. This time it is the volcano Bárðarbunga which appears to be ready to erupt. The Icelandic API web site http://apis.is is providing a useful service, by tracking all measured earthquakes and formatting them as JSON output, which is easily imported into the Wolfram Language. In this post I will make use of this data to visualize these earthquakes.

The first step is to import the dataset and to transform it into a Wolfram Language data structure aptly named Dataset.

dataset = Dataset[Map[Association, "results" /. Import["http://apis.is/earthquake/is", "JSON"]]];

We can inspect the content of this earthquake data by simply evaluating the dataset symbol:

dataset

dataset

And we can also sort the data on any column, by using the SortBy function on the depth column. In this case I use -depth to sort from deepest earthquake to shallowest:

dataset[SortBy[-#depth &]]

dataset sorted

Now, let's make a geographical map of the area around Bárðarbunga. Here, we can use the natural language input which can be accessed by pressing the Ctrl key and the = key simultaneously. This brings up the natural language input field:

natural language input

Now type in that field and hit the Enter key (followed by Shift-Enter to evaluate the input cell):

natural language result

We can do this anywhere in the input, so now we can create a map of the area near Bárðarbunga. We use the GeoGraphics function and place a GeoMarker at the center of the volcano, using a iconic image for a volcano. We set the GeoRange to 100km to see the surrounding area including a piece of Iceland's coastline. For the map background we use the contour map styling:

background map

To get all the earthquake positions as little round disks, we extract the latitude and longitude from the dataset as values and then wrap them in GeoDisk with a radius of 1km:

earthquakes = Normal[ dataset[ All, {"latitude", "longitude"} /* Values /* (GeoDisk[#, Quantity[1, "Kilometers"]] &)]];

And then plot them as red disks on the previous map:

map

To make a slightly more refined image we can composite a relief map with a contour map and overlay the earthquake data (here volcanoimage is the image used above to indicate the volcano location):

reliefmap = GeoGraphics[
Entity["Volcano", "Bardarbunga"], 
  GeoRange -> Quantity[100, "Kilometers"], 
  GeoBackground -> GeoStyling["ReliefMap"], ImageSize -> Large];

contourmap = 
 GeoGraphics[Entity["Volcano", "Bardarbunga"], 
  GeoRange -> Quantity[100, "Kilometers"], 
  GeoBackground -> GeoStyling["ContourMap"], ImageSize -> Large];

overlay = 
 GeoGraphics[{{Red, earthquakes}, 
   GeoMarker[Entity["Volcano", "Bardarbunga"], volcanoimage]}, 
  GeoRange -> (GeoRange /. Options[contourmap, GeoRange]), 
  GeoBackground -> None, 
     ImageSize -> Large];

Overlay[{ImageCompose[reliefmap, {contourmap, .5}], overlay}]

final map

POSTED BY: Arnoud Buzing
Answer
2 years ago

Wow!

POSTED BY: Marco Thiel
Answer
2 years ago

And Double-Wow!

POSTED BY: David Keith
Answer
2 years ago

This is truly great - thanks for posting!

The Wikipedia page says:

Seismic activity surrounding the Bárðarbunga volcano has been gradually increasing since 2007 with a brief respite during the eruption at Grímsvötn in 2011. Activity has now (2014) reached a level similar to that just before the Grímsvötn eruption.

Is it possible to plot time series of average magnitude of the earthquakes since 2007? I tried but data seem to include just few last days:

DateListPlot[
 Transpose[
  {Normal[dataset[All, {"timestamp", "size"}]][[All, 1]],
   Normal[dataset[All, {"timestamp", "size"}]][[All, 2]]}
  ],
 PlotRange -> All]

enter image description here

Also - shouldn't Wolfram Alpha has all these data too somehow? I do not see this "gradually increasing since 2007"

enter image description here

POSTED BY: Sam Carrettie
Answer
2 years ago

Very nice, indeed!

Would it be possible to animate this map so that the earthquakes pop up over time?

POSTED BY: Nils Rune Bodsberg
Answer
2 years ago

Yes, using the techniques described above but saving image for every minute: https://vimeo.com/104637028

And how to generate the animation.

POSTED BY: Arnoud Buzing
Answer
2 years ago

Just for fun, here's 3d version:

reliefmap = 
  GeoGraphics[Entity["Volcano", "Bardarbunga"], 
   GeoRange -> Quantity[100, "Kilometers"], 
   GeoBackground -> GeoStyling["ReliefMap"]];

latlong = 
  GeoBoundingBox[
   GeoDisk[Entity["Volcano", "Bardarbunga"], 
    Quantity[100, "Kilometers"]]];
alt = Reverse[
   QuantityMagnitude[
    GeoElevationData[latlong, UnitSystem -> "Metric"]]];

bardabunga = 
  ListPlot3D[alt, 
   PlotStyle -> Directive[Opacity[.7], Texture[reliefmap]], 
   Mesh -> 10, MeshFunctions -> {#3 &}, 
   TextureCoordinateFunction -> ({#1, #2} &), 
   MeshStyle -> GrayLevel[.5], Lighting -> "Neutral", 
   MaxPlotPoints -> 120, Boxed -> False, Axes -> False, 
   DataRange -> Reverse[Transpose[latlong[[All, 1]]]]];

data = "results" /. Import["http://apis.is/earthquake/is", "JSON"];

data3d = Map[{"size", {"longitude", "latitude", -100 "depth"}} /. # &,
    data];
minmax = Through[{Min, Max}[data3d[[All, 1]]]];

earthquakes = {Point[#2, 
      VertexColors -> (Hue[0, #] & /@ 
         Rescale[#1, minmax, {.1, 1}])]} & @@ Transpose[data3d];

prange = bardabunga // PlotRange;

prange[[3, 1]] = Min[data3d[[All, 2, 3]]];

Show[{bardabunga, Graphics3D[{PointSize[.004], earthquakes}]}, 
 ImageSize -> 720, PlotRange -> prange, BoxRatios -> {1, 1, .7}, 
 Background -> Black]

enter image description here

POSTED BY: Jaebum Jung
Answer
2 years ago

Group Abstract Group Abstract