Continuing our previous post on overlaying the earthquake location data on topographic map, this time we will overlay the earthquake data on a three-dimensional map. The three-dimensional map will be very similar to the one we built using PyGMT in one of our previous post. Read the data fileTaiwan is located at a compressive tectonic boundary between the Eurasian Plate and the Philippine Sea Plate. The island straddles the compressional plate boundary between the Eurasian Plate to the west and the north and the Philippine Sea Plate to the east and south. The present convergence rate is about 7 cm/year in a NW-SE direction. The data for this post was downloaded from Central Weather Bureau of Taiwan, and saved as We can read the tabular data using pandas’ import pandas as pd# dataformat: 1 1995-07-05 17:33:48.600 24.8313 122.1158 27 4.87data = pd.read_csv('evtlist.dat', names=['SN','year','time','latitude','longitude','depth_km','magnitude'], delimiter='\s+') Load the topographic dataWe load the topographic data directly using the PyGMT minlon, maxlon = 119.5, 124.0minlat, maxlat = 19.8, 26.0# Load sample earth relief datagrid = pygmt.datasets.load_earth_relief(resolution="03s", region=[minlon, maxlon, minlat, maxlat]) Define three-dimensional perspectiveWe first specify the cpt files for the topographic colors, and then we specify the three-dimensional perspective view. The view is set to have the azimuth and elevation angle of the viewpoint to be 150 and 30 degrees. The cover type of the grid is taken as “surface”. The other options are mesh, waterfall, image, etc. See the We define the region by the bounds of latitudes, longitudes and depths. The projection was taken to be “Mercator” with width of 15cm. The vertical height was taken to be 4cm. pygmt.makecpt(cmap='geo',series=f'-{maxdep}/4000/100',continuous=True)fig = pygmt.Figure()fig.grdview(grid=grid,region=[minlon, maxlon, minlat, maxlat, -maxdep, 4000],perspective=[150, 30],frame=frame,projection="M15c",zsize="4c",surftype="i",plane=f"-{maxdep}+gazure",shading=0,# Set the contour pen thickness to "1p" contourpen="1p",) Plot the data on a topographic mapSince the depths of the event ranges up to 180 km, we had to constrain it to show it within the range of our chosen maximum depth for the three-dimensional map. We project the depths in the range of 0 to the maximum depth of the displayed topography. Another good choice for displaying the actual depth for each earthquake events is displaying on a cross-sectional map. The advantage of displaying the events on the topographic map is understanding its relation with the tectonics (topography or bathymetry). depthData = data.depth_km.valuesdepthData = maxdep*(depthData - data.depth_km.min())/(data.depth_km.max()-data.depth_km.min()) #normalize the depth to the maximum of maxdep# colorbar colormappygmt.makecpt(cmap="jet", series=[ depthData.min(), depthData.max()])fig.plot3d(x=data.longitude,y=data.latitude,z=-1*depthData,size=0.05*data.magnitude, #scale the magnitude color=depthData,cmap=True,style="cc",pen="black",perspective=True,) Add colorbar to the figureAt last, we can add the colorbar to show the values of actual depths of the event. pygmt.makecpt(cmap="jet", series=[ data.depth_km.min(), data.depth_km.max()])fig.colorbar(frame='af+l"Depth (km)"', perspective=True,) Complete scriptReferences
|
|
来自: LibraryPKU > 《制图》