So an earthquake has happened... you have see some seismograms using the excellent Python Obsy toolkit and some simple scripts... now you want to plot a map or two.
 |
| Simple location map created in Python using the Folium library |
import folium
# setting the boundaries of the map and choosing a map style
m = folium.Map(location=[50.067, -5.187],zoom_start=10,tiles='Stamen Terrain')
# put a marker on the map for the reported BGS earthquake location
folium.Marker(location=[50.067,-5.187],popup='quake',icon=folium.Icon(color='red')).add_to(m)
# next add markers for the locations of all the seismic monitoring stations in Cornwall
folium.Marker(location=[50.1867,-5.2273], popup='BGS_CCA1', icon=folium.Icon(color='green')).add_to(m)
folium.Marker(location=[50.1486,-5.0945], popup='FL_RB30C', icon=folium.Icon(color='green')).add_to(m)
folium.Marker(location=[50.117, -5.535], popup='PZ_RB5E8', icon=folium.Icon(color='green')).add_to(m)
folium.Marker(location=[50.2344, -5.2384], popup='RD_RD93E', icon=folium.Icon(color='green')).add_to(m)
folium.Marker(location=[50.2596, -5.1027], popup='RL_R82BD', icon=folium.Icon(color='green')).add_to(m)
folium.Marker(location=[50.2609,-5.0434], popup='TS_R7FA5', icon=folium.Icon(color='green')).add_to(m)
# save the map to a n html file that can be viewed in any browser
m.save('simple_cornwall_map.html')
Perpendicular Bisectors
When we plotted the seismograms for these stations it was easy to see which station they arrived at first (station CCA1). It follows that the event location must be closer to this station than any other. If we consider two stations, e.g. CCA1 and PZ (penzance) the range of possible locations for the earthqauke must be in all the places that are closer to CCA1 than PZ. This space is bounded by the perpendicular bisector of a line that joins the two stations. If we carry out this procedure for a range of stations then we can narrow down the range of possible locations to a polygon.
 |
| Estimating location using perpendicular bisectors |
import folium, utm, math
map = folium.Map(location=[50.187, -5.23],zoom_start=10,tiles='Stamen Terrain')
folium.Marker(location=[50.067,-5.187], popup='BGS_quake_location',
icon=folium.Icon(color='red')).add_to(map)
# station statistics
seislist=['RB30C','RB5E8','RD93E','CCA1','R82BD','R7FA5']
locations=['FL','PZ','RD','CM','RL','TS']
latitudes=[50.1486,50.117,50.2344,50.1867,50.2596,50.2609]
longitudes=[-5.0945,-5.535,-5.2384,-5.2273,-5.1027,-5.0434]
arrival=[ 7.770, 8.46, 8.51,7.377, 9.413, 9.799]
easting=[]
northing=[]
zone=[]
band=[]
for i in range (len(seislist)):
x, y, z, b = utm.from_latlon(latitudes[i], longitudes[i])
easting.append(x)
northing.append(y)
zone.append(z)
band.append(b)
#maxstation = len(seislist) # determines which bisectors get plotted
maxstation=len(seislist)-2 #PDS edit miss off CM RL and TS
for stat1 in range(maxstation-1):
for stat2 in range(stat1+1, maxstation):
if arrival[stat1]<arrival[stat2]:
s1 = stat1
s2 = stat2
else:
s1 = stat2
s2 = stat1
# join the stations with a line
folium.PolyLine([(latitudes[s1],longitudes[s1]), (latitudes[s2],longitudes[s2])],
color="blue", weight=1, opacity=0.5).add_to(map)
# gradient of perpendicular bisector = -1/gradient of line joining stations
if northing[s1]-northing[s2] != 0:
m = (easting[s2]-easting[s1])/(northing[s1]-northing[s2])
else:
m = (easting[s2]-easting[s1])/0.0001
# calculate utm coordinate of midpoint and equation of line
y = (northing[s1]+northing[s2])/2
x = (easting[s1]+easting[s2])/2
c = y - m * x
# fix length of extrapolated line in x or y direction, depending on m
d = 30000
# calculate distance of shading from line
diff = math.sqrt(15000*15000/(1+m*m))
if northing[s1]>northing[s2]:
dy = -1 * diff
else:
dy = diff
if easting[s1]<easting[s2]:
dx = diff * math.sqrt(m*m)
else:
dx = -1 * diff * math.sqrt(m*m)
if m*m < 1:
x1 = x-d
y1 = m * (x-d) + c
x2 = x+d
y2 = m * (x+d) + c
else:
x1 = ((y-d)-c)/m
y1 = y-d
x2 = ((y+d)-c)/m
y2 = y+d
# plot the perpendicular bisector
folium.PolyLine([utm.to_latlon(x1, y1,zone[s1],band[s1]), utm.to_latlon(x2, y2,zone[s1],band[s1])],
color="black", weight=1, opacity=0.5).add_to(map)
# plot the shaded polygon
folium.Polygon([utm.to_latlon(x1, y1,zone[s1],band[s1]), utm.to_latlon(x2, y2,zone[s1],band[s1]),
utm.to_latlon(x2+dx, y2+dy,zone[s1],band[s1]), utm.to_latlon(x1+dx, y1+dy,zone[s1],band[s1])],
color='black', fill=True, weight=0, opacity=1).add_to(map)
for x in range (len(seislist)):
folium.Marker(location=[latitudes[x],longitudes[x]], popup=[locations[x]], icon=folium.Icon(color='green')).add_to(map)
map.save('cornwall_map4.html')
Intersecting Circles
If it is possible to identify P and S waves from the seismograms you can estimate how far away the earthquake happened from each station by usinge teh delay time between P and S wave arrivals and knowweledge fo P and S wave evlocities
 |
| Estimating the distance from each seismic station using P-S delay times |
import folium
map = folium.Map(location=[50.187, -5.23],zoom_start=10,tiles='Stamen Terrain')
folium.Marker(location=[50.067,-5.187],
popup='BGS_quake_location',
icon=folium.Icon(color='red')
).add_to(map)
# station statistics
seislist=['CCA1','RB30C','RB5E8','RD93E','R82BD','R7FA5']
locations=['CM','FL','PZ','RD','RL','TS']
latitudes=[50.1867,50.1486,50.117,50.2344,50.2596,50.2609]
longitudes=[-5.2273,-5.0945,-5.535,-5.2384,-5.1027,-5.0434]
P_arrival=[7.653, 7.885, 8.441, 8.491, 9.422, 9.802]
S_arrival=[9.57,9.99,10.93,11,12.7,13.5]
radius=[0,0,0,0,0,0]
colours=['red', 'orange', ' beige', 'green', 'blue', 'purple']
for x in range(6):
## estimate range for delay time between P and S arrival (P=5.76km/sec S=3.25km/sec)
radius[x]=(S_arrival[x]-P_arrival[x])*7.48
folium.Marker(location=[latitudes[x],longitudes[x]],
popup=seislist[x], icon=folium.Icon(color='blue')).add_to(map)
folium.Circle(
radius=radius[x]*1000, location=[latitudes[x],longitudes[x]],
popup=locations[x], color=colours[x], fill=False, weight=1
).add_to(map)
map.save('cornwall_map_circles.html')
Comments
Post a Comment