Group Abstract Group Abstract

Message Boards Message Boards

Visualizing Hawaii Seismic Activity with EarthquakeData

Null Null
Posted 8 years ago

If you've been following the news, you've probably seen that the Kilauea volcano in Hawaii has erupted, and there have been reports of many earthquakes. I thought I'd use the built-in EarthquakeData in the Wolfram Language to see if I could come up with any interesting visualizations to explore the seismic activity. I'm going to be using some very basic examples from the documentation.

Let's first grab some data from 5/1/18 to 5/6/18 and plot locations.

data = EarthquakeData[
    Entity["AdministrativeDivision", {"Hawaii", "UnitedStates"}], 
    4, {{2018, 5, 1}, {2018, 5, 6}}, "Position"]["Values"];

GeoGraphics[{GeoStyling["ReliefMap", MaxPlotPoints -> 300], Red, 
  data /. GeoPosition[{x_, y_}] :> Point[{y, x}]}]

enter image description here

Wolfram Language returned a message indicating there were duplicates, so some values were combined.

Let's use an example from the documentation and examine earthquake depth.

data2 = {#["Position"], #["Depth"]} & /@ 
   Values[EarthquakeData[
     Polygon[Entity[
       "AdministrativeDivision", {"Hawaii", "UnitedStates"}]], 
     4, {{2018, 5, 1}, {2018, 5, 6}}]];

Graphics3D[{Opacity[0.6], 
  Map[Append[Reverse[#], 0] &, 
   EntityValue[
     Entity["AdministrativeDivision", {"Hawaii", "UnitedStates"}], 
     "Polygon"] /. GeoPosition -> Identity, {-2}], Red, Opacity[1], 
  Line[{Append[Reverse[First[#1]], 0], 
      Append[Reverse[First[#1]], -QuantityMagnitude[#2]/10]} & @@@ 
    data2]}, Axes -> True]

enter image description here

Pretty impressive for just using some basic examples from the documentation. I'm not an expert in programming, data visualization or geocomputation, so I'm curious what some of you might be able to come up with!

Attachments:
POSTED BY: Null Null
4 Replies
Sort By:

Nice task! Here comes my late night first attempt:

First I wanted to see how in time and space (and strength) the shocks occur:

(* get "raw data" as Dataset:  *)
    eqd0 = Dataset@
   Values[EarthquakeData[Polygon[Entity["AdministrativeDivision", {"Hawaii", "UnitedStates"}]], 4, {{2018, 5, 1}, {2018, 5, 6}}]];
(* extract interesting info: *)
    eqd1 = eqd0[All, {"Depth", "Position", "Magnitude", "Period"}];
(* cartesian coordinates and strength, sorted by time:  *)
    data0 = Normal[Values@eqd1[All, {"Depth" -> ({0, 0, -#} &)}]];
data1 = Most /@ SortBy[{First@GeoPositionXYZ@GeoPositionENU[#1, #2], #3, #4} & @@@ data0, Last];
Graphics3D[{Arrow[#, 1500] & /@ Partition[First /@ data1, 2, 1], Sphere[#1, 300 #2] & @@@ data1}]

which gives:

enter image description here

(* relation to the surface: *)
pos = Normal@eqd1[All, "Position"];
Show[GeoListPlot[pos, 
  GeoServer -> "https://tile.opentopomap.org/`1`/`2`/`3`.png", 
  Joined -> True], 
 GeoGraphics[GeoDisk[GeoPosition[Mean[First /@ pos]], 30000]]]

enter image description here

It surely would be nice to see all in one picture -- regards!
POSTED BY: Henrik Schachner
Attachments:
POSTED BY: Kyle Keane

@Henrik, here you go:

enter image description here

$textureGeoZoomLevel = 10;
$elevationGeoZoomLevel = 7;

$area = {Transparent,    GeoCircle[Entity["Volcano", "Kilauea"],    Quantity[100, "Kilometers"]]};

elevationData =    GeoElevationData[$area, "Geodetic", "GeoPosition", 
  GeoZoomLevel -> $elevationGeoZoomLevel]

geoBounds = elevationData // Point // GeoBounds

eqData = EarthquakeData[$area[[2, 1]],   4, {{2018, 5, 1}, {2018, 5, 10}}];

texture = GeoGraphics[{$area, Point[#Position & /@ eqData]},
   GeoRangePadding -> None,
   GeoRange -> geoBounds,
   GeoProjection -> "Mercator",
   GeoServer -> "http://tile.stamen.com/terrain/`1`/`2`/`3`.jpg",
   GeoZoomLevel -> $textureGeoZoomLevel,
   ImageSize -> 800
 ];

eqPrim = Module[{pos, scaledMag},
     pos =       Append[{#2, #} & @@ #Position[[1]], -QuantityMagnitude[#Depth,         "Meters"]];
     scaledMag = Rescale[#Magnitude, {4, 7}];
     {Thick, Blend["TemperatureMap", scaledMag], 
      Tooltip[HalfLine[pos, {0, 0, 1}], Dataset@#], 
      AbsolutePointSize[20 scaledMag], Point[pos]}
 ] & /@ Values[eqData];

Show[
 ListPlot3D[{#2, #, #3} & @@@ Catenate@First@elevationData,
  PlotStyle -> Texture@SetAlphaChannel[Image@texture, .8],
  Lighting -> "Neutral",
  ImageSize -> 800,
  MeshFunctions -> {#3 &},
  Background -> Black,
  SphericalRegion -> True
  ],
 Graphics3D@eqPrim,
 PlotRange -> All
]

enter image description here
POSTED BY: Kuba Podkalicki

@Kuba Podkalicki : This is awesome - thanks for sharing! Nice idea to use Texture!

Regards -- Henrik
POSTED BY: Henrik Schachner
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard