GMT课程lesson6

葡萄

Lesson 6
Goals

Geological and geophysical analyses are often all about comparing disparate data sets. For example, a problem might involve understanding the relationships between volcano locations and faults, or might involve comparing variations in the gravity field and the distributions of faults and volcanoes. In potential field geophysics, we often compare the spatial distribution of anomalies with other mapped geologic or geophysical features. This kind of comparison helps sharpen interpretation of the potential field anomalies. An essential step in such comparisons is to make sure that the data are plotted together accurately. Obviously the quality of the interpretation depends on the quality of the mapping!

The goal of this lesson is to learn about plotting different types of geologic and geophysical together on the same map using GMT.

---

Seismic Tomographic Anomalies and Volcanoes

Over the last twenty years, geophysicsts have used seismic tomography to image variations in the Earth's crust and mantle in more detail than was previously possible. Seismic tomography works in a way similar to a CAT scan used to image the human body. Seismometers are deployed over a region of interest and these seismometers monitor either natural earthquakes or explosions (or other energy sources). Each seismometer senses the arrival times of seismic waves (for example the P (compressional) wave) associated with each seismic event. Using a simple model of the earth, one might predict the variation of P-wave arrival times among all the seismometers in the region (an exercise that is not too abstract when you realize that rate x time = distance, and that given the source and the seismometer location, it is possible to predict the travel path of the seismic wave). Of course, actual arrival times of P-waves at the seismometers will differ from arrival times predicted by a simple model.

In tomography, the earth is discretized into a series of cubic regions (voxels). The speed of P-waves within each vowel is assumed to be constant. Then using linear algebra techniques, the "best-fit" velocity distribution among all the voxels can be determined that best explains the variation in P-wave arrival times observed by the seismic network.

The variations in P-wave velocity are produced by changes in the elastic properties of rocks. One way to produce slow velocity regions is to develop partial melting along grain boundaries. Consequently, you might expect to see slow velocity regions correlating with evidence of active volcanism.

One of the best seismic tomographic images of the lower crust and mantle comes from northern Japan, where a dense seismic network is deployed. Hasegawa, Zhao, Nikajima, Tamura and others have been able to use these seismic tomographs to create realistic models of subduction zones. The following figure shows contours of P-wave velocity anomalies in northern Japan, using a tomographic model published by Zhao and others. The figure contours P-wave velocity variation along voxels located at a depth of approximately 40 km.

葡萄

On this map P-wave anomalies in voxels at 40 km depth are contoured and imaged with a color table, as was illustrated in Lesson 5. Volcanoes from the "Quaternary Volcano Database for Japan (1999) are superimposed on the map as black solid triangles. Note that although the correlation is not perfect, slow P-wave regions (blue) correlate with areas of abundant Quaternary volcanism, suggesting that volcanism is related to areas of partial melting in the mantle.


Making the Map in GMT
The tomographic data used to make the map are located in the file Vp.dat. Download this file and use the "minmax" program in GMT to discover the limits of the data. The "z" values of the x,y,z triplets in the file are the P-wave velocities, given in kilometers per second.
The file edited_vents.dat contains a partial list of Quaternary volcanoes from this part of Japan. The map is generated using the script: Vp_map.gmt". Play around with this script and convince yourself you know how it works. Try changing the color table. What does the "pscoast" application do to clip the contour lines of seissmic velocity at the coastline? In this case, the clipping makes sense because the tomographic model was constructed using seismometers located on land. Try turning the clipping off. Note that because all of the data are given in latitute longitude coordinates, it is a simple matter to plot them together in GMT. At some scales, sigificant errors can creep in if the data sets use different map datumes (e.g., WGS84, NAD27, etc.). It is a good idea tokeep track of map datums when possible, although sometimes they are not reported at all. Note on the map plotted above that the contour lines do not always correlate with changes in color. Experiment with the script to try to get these to match better, producing a more pleasing figure.
Assignment 6
Download the data file: grav_lat_long.xyz (4 Mb). This file contains the simple Bouguer gravity anomaly for the northern part of Japan, for the same area as the Vp data. First, search the web and learn a little about simple Bouguer gravity anomalies. Modify the script Vp_map.gmt to contour gravity anomalies rather than tomographic anomalies. Compare the volcano distribution to the distribution of volcanoes (plot them on the same map!). Do the same correlations hold? Copy your maps and scripts into your running word document.

Chuck Connor
Last modified: Thu Jan 12 11:25:19 EST 2006

===============================================================

dos批处理：

REM set the frame width to 1.0 point - not really needed

gmtset FRAME_PEN 2.0p
gmtset HEADER_FONT 5
gmtset LABEL_FONT 5
gmtset HEADER_FONT_SIZE 20
gmtset LABEL_FONT_SIZE 16
gmtset ANOT_FONT_SIZE 14
REM surface does the interpolation with a minimum curvature algorithm
REM -I2  means the grid spacing is 2 units in each direction
REM -Gfilename indicates the output file
REM -R specifies the west,east,south, and north bounds of the map
surface -I0.1 -V Vp.dat -GVp.grd -R138/142/37/42

REM makecpt creates a color table
REM -Cno_green specifies the basic color table
REM -T specifies the range and interval
makecpt -Chaxby -T7.2/8.2/0.05 -V > Vp.cpt

REM grdimage plots the image (color map)
REM -JX6.0i is the scale (must match the following)
REM -Ccn.cpt is the color scale created with makecpt
REM -P portait mode (must match the following)
REM -E is the dpi of the color shading
REM -K more postscript to be appended to cn.ps in the following
grdimage Vp.grd -R -JM6i -CVp.cpt -P -E100 -K -V -Y3.0 > Vp.ps

REM grdcontour draws the contours from the grid
REM -JX6.0i means the map will be 6 inches wide
REM -C250 means there is a 250 nT contour interval
REM -A500 means the contours are annotated every 500 nT
REM -B25a50f5/WSne draw the frame, 25 m tick with 50 m label, add 5 m ticks, label on
REM west and south side only
REM -W0.25p set line width
REM -P draw in portrait mode
REM -O overlay contours on the image
grdcontour Vp.grd -JM6i -C.05 -A0.1 -B1a1:"Longitude":/:"Latitude":WSne -W0.25p -P -V -O -K  >> Vp.ps

REM plots volcano locations
REM switch -: reverses the longitude/latitude order
psxy edited_vents.dat -JM6i -R -W1.0 -V -O -P -St0.085i -G0 -K >> Vp.ps

pscoast -JM6i -R -B2g1 -N1/5.0 -Dh -W5.0 -S255 -Lfx10/-1/10/100 -K -O -V -P >> Vp.ps

psscale -D3/-.8/5/0.4h -CVp.cpt -B0.20/:"km/s": -O -P -I -V >> Vp.ps

del .gmt*

效果图：

葡萄

后续作业
REM set the frame width to 1.0 point - not really needed
gmtset FRAME_PEN 2.0p
gmtset HEADER_FONT 5
gmtset LABEL_FONT 5
gmtset HEADER_FONT_SIZE 20
gmtset LABEL_FONT_SIZE 16
gmtset ANOT_FONT_SIZE 14
REM surface does the interpolation with a minimum curvature algorithm
REM -I2  means the grid spacing is 2 units in each direction
REM -Gfilename indicates the output file
REM -R specifies the west,east,south, and north bounds of the map
surface -I0.1 -V grav_lat_long.xyz -Ggra.grd -R138/142/37/42
REM makecpt creates a color table
REM -Cno_green specifies the basic color table
REM -T specifies the range and interval
grd2cpt gra.grd -Chaxby -Z -I -V > gra.cpt
REM grdimage plots the image (color map)
REM -JX6.0i is the scale (must match the following)
REM -Ccn.cpt is the color scale created with makecpt
REM -P portait mode (must match the following)
REM -E is the dpi of the color shading
REM -K more postscript to be appended to cn.ps in the following
grdimage gra.grd -R -JM6i -Cgra.cpt -P -E100 -K -V -Y3.0 > gra.ps
REM grdcontour draws the contours from the grid
REM -JX6.0i means the map will be 6 inches wide
REM -C250 means there is a 250 nT contour interval
REM -A500 means the contours are annotated every 500 nT
REM -B25a50f5/WSne draw the frame, 25 m tick with 50 m label, add 5 m ticks, label on
REM west and south side only
REM -W0.25p set line width
REM -P draw in portrait mode
REM -O overlay contours on the image
grdcontour gra.grd -JM6i -C5 -A25+p+g255 -S10 -B1a1:"Longitude":/:"Latitude":WSne -W0.25p -P -V -O -K  >> gra.ps
REM plots volcano locations
REM switch -: reverses the longitude/latitude order
psxy edited_vents.dat -JM6i -R -W1.0 -V -O -P -St0.085i -G0 -K >> gra.ps
pscoast -JM6i -R -B2g1 -N1/5.0 -Dh -W5.0 -S255 -Lfx10/-1/10/100 -K -O -V -P >> gra.ps
psscale -D3/-.8/7/0.4h -Cgra.cpt -B50/:"mgal": -O -P -I -V >> gra.ps
del .gmt*