Message Boards Message Boards

GROUPS:

Visualizing Hawaii Seismic Activity with EarthquakeData

Posted 4 months ago
866 Views
|
4 Replies
|
25 Total Likes
|

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:
4 Replies

It is an interesting topic for computational exploration Swede, I did a little 60 minute riff from your prompt and came up with something that I think points in an interesting direction of 3D visualization of Earthquake data - although it needs some beautification and better scaling. If I had some more time I would use the map of Hawaii as a texture and apply it to the BSplineSurface[] of the islands' GeoElevationData[], but alas other duties call.

This reminds me how fun Wolfram Language is for exploring the world around me. It is pretty darn cool that EarthquakeData[] seems to be pulling data nearly as fast as it is being reported. I made this so it pulls the last three days worth of earthquake data, but this could easily be spun into a FormPage[] that lets people explore the region during any time period they want.

And of course, lots of compassion to those being affected by the recent volcanic activity!

Here is the output of my exploration, it is a 3D topographical map of the Hawaiian Islands along with a bubble for each earthquake - radius of each bubble is determined by the magnitude of the quake: geodesic 3D plot of Hawaii and recent earthquake activity

Here is how I got there (notebook attached for those that want it)...

pull data from last three days about earthquake activity near the big island

data=EarthquakeData[ Entity["Island", "HawaiiIsland"], 2, {DatePlus[Now,{-3,"Days"}],Now}];

Do some cleaning, could be done more tersely using Query[], but I am old school in some ways

simpleData=Values/@data[[All,{"Period","Magnitude","Position","Depth"}]];

simpleXYZ = MapAt[First@GeoPositionXYZ@# &, simpleData, {All, 3}];

simpleXYZ=Values@MapAt[QuantityMagnitude@#&,simpleXYZ,{All,4}];

Make the visuals for the events (could have done something clever with "Depth", but ran out of time).

vizMe=Map[Function[Ball[#[[3]],#[[2]] 10^3]],simpleXYZ];

Next, I make a 3D plot of Hawaii's elevation

elevation = 
  First@GeoElevationData[Entity["AdministrativeDivision", {"Hawaii", "UnitedStates"}],  Automatic, "GeoPositionXYZ", UnitSystem -> "Metric"];

surface = BSplineSurface@elevation;

and voila, you can see a simple map of all the Earthquakes

Graphics3D[{surface, vizMe}]

geodesic 3D plot of Hawaii and recent earthquake activity

Attachments:

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!

@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

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

Regards -- Henrik

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract