xxxxxxxxxx
 
logo(); logo!(GMTjulia=1.5, x_off=5, show=true)

2.9 s

Working with ISC Data Set using the GMT Julia wrapper

Note: The examples in this Notebook are a reproduction of MATLAB® for Analyzing and Visualizing Geospatial Data availabe in Matlab's File Exchange.

In this example, we will load an historical data set, earthquake hypocenters from the ISC-GEM Catalogue

Preview the data

Having looked in advance, I know that this file has a lot of lines of explanation and header information up front, followed by tabular data. The visual file inspection also sows that the column names in CVS file start at line 93. We need to transmit this info to the CSV reader. And also want to load the data in a DatFrame and since we will be dealing with dates we need to import a couple packages. Let's do it and look at data. Note that we are providing the full file name to make this independent on your file location (but this must be adapted for user case, off course)

4.7 ms

x
 
using GMT, Dates, DataFrames, CSV

36.3 s

isc

	_date	lat	lon	smajax	sminax	strike	q	depth
	DateTime	Float64	Float64	String	String	String	String	Float64
1	1904-04-04T10:02:34.560	41.802	23.108	" 8.6 "	" 6.6 "	" 164.2 "	" B "	15.0
2	1904-04-04T10:26:00.880	41.758	23.249	" 8.3 "	" 6.9 "	" 15.2 "	" B "	15.0
3	1904-06-25T14:45:39.140	51.424	161.638	" 33.6 "	" 18.7 "	" 116.2 "	" C "	15.0
4	1904-06-25T21:00:38.720	52.763	160.277	" 28.6 "	" 14.6 "	" 43.1 "	" C "	30.0
5	1904-12-20T05:44:20.440	8.962	-84.042	" 26.3 "	" 13.2 "	" 74.1 "	" C "	10.0
6	1905-02-14T08:46:34.490	51.89	-176.702	" 30.3 "	" 19.5 "	" 157.5 "	" C "	35.0
7	1905-02-17T11:41:07.820	23.689	97.17	" 27.1 "	" 18.9 "	" 143.7 "	" C "	15.0
8	1905-02-19T04:34:55.800	-14.123	169.457	" 54.7 "	" 12.1 "	" 3.4 "	" C "	15.0
9	1905-03-04T15:59:53.260	-6.639	151.225	" 82.1 "	" 33.3 "	" 161.3 "	" C "	15.0
10	1905-03-04T23:17:22.030	-4.635	149.105	" 22.4 "	" 18.4 "	" 155.1 "	" C "	15.0
more
39160	2016-12-30T12:56:12.430	-30.546	-177.749	" 6.7 "	" 4.9 "	" 108.0 "	" A "	25.1

xxxxxxxxxx
 
isc = CSV.File("c:/v/isc-gem-cat.csv", header=93, normalizenames=true, dateformat=Dates.DateFormat("yyyy-mm-dd HH:MM:SS.sss")) |> DataFrame

15.0 s

Select a subset of information.

We now select the subset of columns that we want to explore from these historical data. You can see that we am choosing earthquake location, magnitude, moment, and quality information about the data.

6.9 μs

	_date	lon	lat	depth	q_1	mw	q_2	mo	mo_auth
	DateTime	Float64	Float64	Float64	String	Float64	String	String	String
1	1904-04-04T10:02:34.560	23.108	41.802	15.0	" C "	6.84	" C "	" 2.30 "	" bibliog "
2	1904-04-04T10:26:00.880	23.249	41.758	15.0	" C "	7.02	" C "	" 4.20 "	" bibliog "
3	1904-06-25T14:45:39.140	161.638	51.424	15.0	" C "	7.5	" C "	" 2.60 "	" bibliog "
4	1904-06-25T21:00:38.720	160.277	52.763	30.0	" C "	7.7	" C "	" 5.10 "	" bibliog "
5	1904-12-20T05:44:20.440	-84.042	8.962	10.0	" C "	7.29	" C "	" "	" "
6	1905-02-14T08:46:34.490	-176.702	51.89	35.0	" C "	7.59	" C "	" "	" "
7	1905-02-17T11:41:07.820	97.17	23.689	15.0	" C "	7.26	" C "	" "	" "
8	1905-02-19T04:34:55.800	169.457	-14.123	15.0	" C "	7.33	" C "	" "	" "
9	1905-03-04T15:59:53.260	151.225	-6.639	15.0	" C "	6.96	" C "	" "	" "
10	1905-03-04T23:17:22.030	149.105	-4.635	15.0	" C "	7.01	" C "	" "	" "
more
39160	2016-12-30T12:56:12.430	-177.749	-30.546	25.1	" A "	5.48	" A "	" 2.09 "	" gcmt "

xxxxxxxxxx
 
ds = isc[:, [:_date, :lon, :lat, :depth, :q_1, :mw, :q_2, :mo, :mo_auth]]

110 ms

Select a subset with only high quality data.

3.4 μs

xxxxxxxxxx
 
eqs = ds[ds.q_1 .== " A ", :];

262 ms

Let's see what we've got.

2.7 μs

xxxxxxxxxx
 
scatter(eqs.lon, eqs.lat, marker=:point, mc="#0072BD", show=true)

6.0 s

Let's make a histogram of the number of earthquakes vs. time.

But first let's see what the timespan is.

4.4 μs

1911-07-12T04:07:38.110

2016-12-30T12:56:12.430

xxxxxxxxxx
 
minDt, maxDt = extrema(eqs._date)

10.7 ms

Now our histogram.

3.3 μs

xxxxxxxxxx
 
histogram(eqs[!, :_date], xlabel=:year, ylabel="# quakes", title="Histogram of Quakes by Year", show=true)

3.4 s

Plot time vs. magnitude

4.1 μs

xxxxxxxxxx
 
scatter([eqs[!, :_date] eqs[!, :mw]], marker=:point, mc=:blue, xlabel=:Year, ylabel=:Magnitude, title="Quake magnitudes vs. time", show=true)

2.3 s

Looking at Earthquakes in Japan

We will look at data inside the region [135 150; 35 45]

6.6 μs

xxxxxxxxxx
 
eqsJ = eqs[(135 .< eqs.lon .< 150) .& (35 .< eqs.lat .< 45), :];

307 ms

Find the Largest Earthquake

9.1 μs

	_date	lon	lat	depth	q_1	mw	q_2	mo	mo_auth
	DateTime	Float64	Float64	Float64	String	Float64	String	String	String
1	2011-03-11T05:46:23.150	142.546	38.285	20.0	" A "	9.09	" A "	" 5.31 "	" gcmt "

xxxxxxxxxx
 
eqsJ[eqsJ.mw .== maximum(eqsJ.mw), :]

170 ms

Plot on a Map of Japan

We first plot the earthquakes on a map to see if they are located in a reasonable position.

We wil scale the symbols by some function of magnitude that allows to see the largests events. The base data for creating the topo/bathymetry map will be automatically downloaded from a GMT server.

5.4 μs

xxxxxxxxxx
 
ptSz = (eqsJ.mw).^2 ./10 .* (2.54/72);

81.5 ms

xxxxxxxxxx
 
begin
    grdimage("@earth_relief_30s", region=(135,150,35,45), cmap=:geo, J=:guess,
        title="Japan earthquakes", shade=true, coast=true, colorbar=true)
    scatter!(eqsJ.lon, eqsJ.lat, markersize=ptSz, mc=:yellow,
        markeredgecolor=:black, show=true)
end

5.4 s

Convert lat, lon, depth to x, y, z values and plot

We are interested in location as a function of depth so to plot the quakes, we first convert our locations to a local Cartesian coordinates. We assume geodetic coordinates relative to a reference ellipsoid (WGS84) to get our x,y,z (enu) values.

4.1 μs

xxxxxxxxxx
 
# Divide by 1000 to get data in km
xEast, yNorth, zDepth = geodetic2enu(eqsJ.lon, eqsJ.lat, -eqsJ.depth*1000, 145, 38, 0) ./ 1000;

18.8 ms

xxxxxxxxxx
 
scatter3(xEast, yNorth, zDepth, ms=ptSz, view=(180,90), zcolor=eqsJ.mw, colorbar=true,
    xlabel="East (km)", ylabel="North (km)", axes=:autoXYZg, aspect=:equal,
    title="Japan earthquakes - overhead view", show=true)

843 ms

Display image from published results to compare.

3.1 μs

xxxxxxxxxx
 
scatter3(xEast, yNorth, zDepth, ms=ptSz, view=(200,5), zcolor=eqsJ.mw, aspect=:equal,
    xlabel="East (km)", ylabel="North (km)", zlabel="Depth (km)", zsize=5,
    marker=:u, frame=:autoXYZg, colorbar=true, show=true)

830 ms

Compare to model

From: http://www.earthobservatory.sg/files/project/images/Tohoku2-blocdiagrammejapan_earthquakes.jpg)

7.8 μs

xxxxxxxxxx
 
imshow("http://www.earthobservatory.sg/files/project/images/Tohoku2-bloc_diagramme_japan_earthquakes.jpg")

23.0 s