Lecture notes, lectures 6 - 10 PDF

Preview

Summary

L06 – Generic Mapping Tools (GMT) - Part 21. Plotting Fields - GrdimageLast time we plotted a couple of maps. This was great, but there wasn’t a lot of information on those maps other than where the coastlines are located. The next example shows how to add some additional information. In this case, ...

Description

Geophysical Computing

L06-1

L06 – Generic Mapping Tools (GMT) - Part 2 1. Plotting Fields - Grdimage Last time we plotted a couple of maps. This was great, but there wasn’t a lot of information on those maps other than where the coastlines are located. The next example shows how to add some additional information. In this case, we are plotting seismic S-wave velocities as derived from the tomographic inversion by Jeroen Ritsema (U. Michigan). The values of S-wave velocity are given (in % difference from a standard reference model – PREM) in the file s20rts.grd (this file is supplied on the webpage along with the data sets for use in the homework). Grab the files available for the homework and try the following example. #!/bin/csh # Generate a color palette table makecpt –Ccelsius –T-5/5/.1 –Z –I >! color.cpt # Plot the gridded image grdimage s20rts.grd –Ccolors.cpt –R120/280/-55/40 –JX6.4i/3.8i \ -B20g10000f10/10g10000nSeW –P –E300 –K >! plot.ps # Add the coastline information pscoast –R –JX6.4d/3.8d –Dc –W1/2/255/255/255 –P –O –K –A10000 \ –N1/1/255/255/255 >> plot.ps # Add a scale bar of the colors psscale –D2.0i/4.75i/3.5i/.3ih –O –Ccolors.cpt –B2.5 >> plot.ps rm colors.cpt gs –sDEVICE=x11 plot.ps The resulting image you should get is shown here:

So, what did we do in creating this image? I am assuming you are expert in searching the man pages by now to determine what all of the flags mean. But, let’s discuss a couple of the main points in creating such an image:

Geophysical Computing

L06-2

1) We need data in a specific type of gridded format. For this example we specifically needed our data in longitude, latitude, S-wave velocity format. In this case it was specified in a special kind of file called a netCDF or grid file. This is the supplied file called s20rts.grd. Later on in this lecture we will talk about how to create this type of file. 2) We needed a file that provided a mapping between seismic wave velocities and the color used to plot them. This is called a color palette table or .cpt file. Here we used the file colors.cpt. The next section of this lecture shows how to generate such a table. 3) Once we have a .grd and .cpt file we can create a colored image using the GMT command grdimage. 4) Note that we can now plot coastline information with pscoast, but that it has to overlay the color information or else we won’t be able to see it. 5) Finally, we should always add a scale bar to our plot so we know what the colors mean. This is done with GMTs psscale command.

2. Color Palette Tables Now, that we’ve seen an example of plotting a gridded image let’s discuss how the coloring of this image is done. Note in the above example that we specified a file: colors.cpt. The .cpt extension is used to denote a color palette table. Here is an example of a cpt file designed specifically for plotting the topography of Scotland: # scotland.cpt # # by Eric Gaba -1750 121 178 222 -1500 132 185 227 -1250 141 193 234 -1000 150 201 240 -750 161 210 247 -500 172 219 251 -200 185 227 255 -100 200 235 255 -50 216 242 254 0 172 208 165 50 148 191 139 100 168 198 143 200 189 204 150 400 209 215 171 600 239 235 192 800 222 214 163 1000 202 185 130 1200 192 154 83 B 0 0 0 F 255 255 255 N 255 0 0

-1500 -1250 -1000 -750 -500 -200 -100 -50 0 50 100 200 400 600 800 1000 1200 1400

121 132 141 150 161 172 185 200 216 172 148 168 189 209 239 222 202 192

178 185 193 201 210 219 227 235 242 208 191 198 204 215 235 214 185 154

222 227 234 240 247 251 255 255 254 165 139 143 150 171 192 163 130 83

What does this mean? Looking at the first line we have:

Geophysical Computing

L06-3

-1750 121 178 222 -1500 121 178 222

What this says is that we want to color elevations between -1750 to -1500 (so this is actually below sea level) with the RGB color: R=121, G=178, B=222. This color looks like: . The next line is: -1500 132 185 227 -1250 132 185 227

Which simply states that we want to color elevations between -1500 and -1250 with the color: R=132, G=185, B=227. So the cpt file is just a mapping between elevations and the color used to represent those elevations. The final three lines of the cpt file are special: B 0 0 0 F 255 255 255 N 255 0 0

Note that the color palette table in our example scotland.cpt is valid for elevations between -1750 and 1400 m. These three lines state: 1) B – Elevations are less than the defined range. In this example if the elevation occurs in the map with a value < -1750 m, then color that space as black (R=0, G=0, B=0). 2) F – Elevations are greater than the defined range. In this example if the elevation is > 1400 m then color that space as white (R=255, G=255, B=255). 3) N – If no elevation data is given for that space (i.e., elevation is given as NaN - Not-ANumber) then color that space as red (R=255, G=0, B=0). This is an example of a categorical type of cpt. What this means is that we do not interpolate color for a given block of data. That is, in our example, for the range of elevation between -1750 and -1500 m, we use a single color. Hence, a scale bar for this type of cpt looks like discrete blocks as shown here:

However, we can also make a continuous cpt, where color is linearly interpolated between data ranges. To do this, just consider the first line of our cpt file: -1750 121 178 222 -1500 121 178 222

Now, what would happen if we wrote: -1750 121 178 222 1400 192 154 83

This line by itself could be a cpt file. Try it out by creating a cpt file with just that line, and plotting a color bar using psscale.

Geophysical Computing

L06-4

As you can see we can make cpt files by just plain hand editing. But, luckily GMT comes with a number of preloaded cpts. To see their names just type makecpt. The makecpt utility is great for generating cpts with custom bounds. For example, what if I wanted to make a continuous color palette table using the GMT base file relief for topography, but having the elevations limited to the range between +/- 1000 feet. We can do this easily by: >> makecpt –Crelief –Z –T-1000/1000/25 >! color.cpt A couple of flags are of special note: 1) –Z: This flag says to make a continuous cpt file. 2) –I: This flag says reverse the colors As a final note, there is an excellent resource called CPT City located at: http://soliton.vm.bytemark.co.uk/pub/cpt-city/ There are literally hundreds of color palette tables contributed to this web site, many by professional cartographers, so it’s often useful to peruse the site for interesting color schemes. Note that the colors here are given in a variety of formats, but the cpt format is the same as that used in GMT.

3. Generating Grid Files (xyz2grd) In general our data sets aren’t arranged in .grd files. But luckily GMT comes supplied with a nice command to create them: xyz2grd. To make a .grd file one typically needs to have data ordered in a table of x, y, and z values. For example, if I have global data set of temperatures on the Earths surface at a specific time, then I would want the table to look like: lon 1 lon 2 lon 3 … lon N

lat 1 lat 2 lat 3 … lat N

T1 T2 T3 … TN

If I named this input file as temperature.xyz, and my longitude and latitude points were spaced on a 1° × 1° interval, then I could generate a .grd file as follows: >> xyz2grd temperature.xyz –Gtemps.grd –R0/360/-90/90 –I1/1 –V Here, the –G flag says that we want to create a grid file named temps.grd. The –I flag tells us the spacing between x and y grid points. In our example we suggested a 1° × 1° interval which is why this flag takes on its current values. The –V flag states we want xyz2grd to spit out as much

Geophysical Computing

L06-5

information as possible while generating the .grd file. This can be really useful in trying to debug problems.

4. Global Datasets The National Oceanic and Atmospheric Administration (NOAA) manage a number of really nice global datasets, which can often be downloaded in netCDF grid file format. The Etopo series of surface topography and bathymetry are probably the most useful data set you will find on-line. http://www.ngdc.noaa.gov/mgg/global/global.html

Other interesting data sets include ocean crustal ages: http://www.ngdc.noaa.gov/mgg/ocean_age/ocean_age_2008.html

Geophysical Computing

L06-6

and Geomagnetic data: http://www.ngdc.noaa.gov/geomag/geomag.shtml

5. What are .grd files? In the older GMT examples, the authors always referred to their gridded data sets as .grd files, or grid files. However, I notice that now they more commonly refer to the files as .nc files. There is no change in the file format, they are still of the form called netCDF. NetCDF (Network Common Data Form) is a set of software libraries and self-describing, machine-independent data formats that support the creation, access, and sharing of array-oriented scientific data. http://www.unidata.ucar.edu/software/netcdf/ This data format is embraced much more in the atmospheric sciences discipline than it is in the geophysical community. But, in the geophysical community everyone still thinks they should design their own data format for every thing they do. Many of the old timers still like to cling to their data as well. But, if we were to adopt a common format for sharing data, netCDF would probably be the format. If you really want details on the netCDF format see my other handout that is on the web page. Here is a reiteration of the important points: NetCDF is an array based data structure for storing multidimensional data. A netCDF file is written with an ASCII header and stores the data in a binary format. The space saved by writing binary files is an obvious advantage, particularly because one does not need worry about the byte order. Any byte-swapping is automatically handled by the netCDF libraries and a netCDF binary file can thus be read on any platform. Some features associated with using the netCDF data format are as follows:

Geophysical Computing

L06-7

Coordinate systems: support for N-dimensional coordinate systems. • X-coordinate (e.g., lat) • Y-coordinate (e.g., lon) • Z-coordinate (e.g., elevation) • Time dimension • Other dimensions Variables: Support for multiple variables. • E.g., S-wave velocity, P-wave velocity, density, stress components… Geometry: Support for a variety of grid types (implicit or explicit). • Regular grid (implicit) • Irregular grid • Points Self-Describing: Dataset can include information defining the data it contains. • Units (e.g., km, m/sec, gm/cm3,…) • Comments (e.g., titles, conventions used, names of variables (e.g., P-wave velocity), names of coordinates (e.g., km/sec),... As a special point, netCDF comes with a utility called ncdump that will allow you to dump the contents of a netCDF file into an ASCII readable format. It also allows you to check just the header information. Try it: >> ncdump –h netcdf_file GMT also comes with a similar utility: >> grdinfo netcdf_file

6. Homework 1) Write a C Shell script utility that one can use to visualize a color palette table. That is it should use the psscale command to generate a plot that one can quickly see what the .cpt file looks like. If we name the script viewcpt then we should be able to run it as: >> viewcpt cptfile where cptfile is the name of the color palette table we are interested in seeing. 2) Using the Etopo1 grid file create a world map to emulate the first plot shown in Section 4. 3) Make a global plot showing elevation of land masses using the Etopo1 elevations. Instead of using bathymetry for the oceans use the ocean ages and overlay the plate boundaries we used in the last homework. This problem will require clipping. Refer to the GMT Example #17 for hints on how to do this. Hint: Elevations will look best here as gray scale.

Geophysical Computing

L06-8

4) The webpage contains a file: Tomography_TXBW.xyz. This file has longitude, latitude, and variations of S-wave velocity based on the tomographic inversion of Steve Grand (University of Texas at Austin). The velocities are given for a layer just above the core mantle boundary. Convert this file into a .grd file and plot it in two panels centered on (a) the central Pacific Ocean region and (b) Africa. It is customary to plot the lowest values (low velocities) with red colors and the highest values (high velocities) with blue colors. Make the color palette table a continuous table. Also plot world hot spot locations as triangles (to emulate volcanoes) on the maps. Be sure to include a scale bar for the colors. Is there any relation between deep mantle seismic velocities and hot spots locations?

Geophysical Computing

L07-1

L07 – Generic Mapping Tools (GMT) - Part 3 1. Local Scale Mapping In our last lecture we showed some examples of global scale mapping. We also introduced a global elevation model Etopo1. This file has a single elevation value for each 1 minute × 1 minute grid point. This is an elevation value for grid sizes about 1.85 km × 1.85 km. This is a really fantastic grid size for plotting at global scales, or even on the statewide scale. Consider the next script which shows an example of plotting elevations for the state of Utah. #!/bin/csh makecpt -CColor_DEM.cpt -Z -T800/5000/100 >! utah.cpt grdimage ETOPO1_Ice_g_gmt4.grd -R-116/-108/36/43 -JX5i/6i \ -B2g10000nSeW -Cutah.cpt -V -P -K >! utah.ps pscoast -R -JX5d/6d -Df -V -P -O -K -W2 -N2 \ -S33/204/255 >> utah.ps # draw an inset box psxy -R -JX -W8/255/0/0 -O -P > utah.ps -111.833 40.75 -111.833 40.583 -111.5 40.583 -111.5 40.75 -111.833 40.75 eof gs -sDEVICE=x11 utah.ps The execution of which results in the following plot:

Geophysical Computing

L07-2

Note that I use a color palette table derived form Color_DEM that was obtained from the CPT City website and available in this lectures accompanying material on the web site. The resulting image on the scale of the state of Utah is quite nice. But notice the red inset box which approximately covers the area between Mill Creek and Big Cottonwood Canyons. If we zoom into this space what we get looks like this:

As you can see the resolution isn’t that pretty at this scale. So, what do we do? Fortunately there are a ton of good data sources for Digital Elevation Models (DEM), ranging from 90 meters × 90 meters (lateral grid size), down to 2 m × 2 m in certain locations. Data sources A fair amount of data is currently available online. However, it is usually on a state by state basis. For example, one can find DEM data for the state of Arizona at: http://sco.az.gov. But, what is available there compared with other states varies quite dramatically. Utah has one of the nicest repositories of DEM data I’ve seen yet. This is available at: http://gis.utah.gov To get to the data click on GIS Data & Resources > SGID, Utah GIS Data Note that here there are two important types of data resources: (1) Vector Data, and (2) Raster GIS Data. We will talk about using both of these file sources.

2. Generating netCDF files from DEM files. We will first look at Raster GIS Data, and work through collecting elevation models and converting this to a .grd file. Just follow the instructions below which are rather recipe like. •

Go to GIS Data & Resources > SGID, Utah GIS Data > Raster GIS Data > Elevation/Terrain Data

•

Next go to: 5 Meter Auto-Correlated Elevation Model (DEM)

•

Next go to: Retrieve 5 meter DEM via Interactive Map.

•

Let’s generate an image of the region between Mill Creek and Big Cottonwood Canyons that will look nicer than the above image. Zoom into the region around Salt Lake City, and notice that the map is broken up into squares. Each square represents a different DEM file.

Geophysical Computing

L07-3

Download the files: 12TVK200800.asc, 12TVK400800.asc, 12TVL200000.asc, and 12TVL400000.asc •

The .asc files we downloaded are not written in the most convenient manner for converting to a .grd file with GMT. In the material for this lecture there is a C Shell script called catascfiles.csh. This is a script written by myself that will combine all of the files we downloaded into one big file. The final output file is named: temp01, and is just a table of easting, northing, and elevation (in meters). The brunt of the work is done with a Fortran 90 code called Lidar_ASC2XYZ.f90. We will cover Fortran 90 later in the semester, but for now just assume this code I wrote works. So, let’s make one big file of x, y, z positions by typing: >> ./catascfiles.csh

•

Now, let’s generate a grid file. Here is an example file on how to do this:

#!/bin/csh # # Example script for generating .grd files # # Input file is 'temp01' written out from catascfiles.csh # #------------------------------------------------------------------# #------------------------------------------------------------------# # Set input parameters #------------------------------------------------------------------# set ifile = temp01 # input file output from catascfiles.csh set elev = f # output elevation in meters (m) or feet (f) set cell = 5 # cellsize (5 meters for the 5-m DEM) set ofile = Mt_Olympus # Prefix naming for output files #------------------------------------------------------------------# #------------------------------------------------------------------# # Determine boundaries of grid file #------------------------------------------------------------------# echo "Determiniming map bounds..." minmax -C $ifile >! temp.f1 set set set set

xmin xmax ymin ymax

= = = =

`awk `awk `awk `awk

'NR 'NR 'NR 'NR

== == == ==

1 1 1 1

{print {print {print {print

$1}' $2}' $3}' $4}'

temp.f1` temp.f1` temp.f1` temp.f1`

rm temp.f1 #------------------------------------------------------------------#

Geophysical Computing

L07-4

#------------------------------------------------------------------# # Convert elevations to feet if switch is set #------------------------------------------------------------------# if ($elev == 'f') then echo "Converting elevations from meters to feet..." awk '{print $1, $2, ($3*3.2808399)}' $ifile >! temp.f2 else cp $ifile temp.f2 endif #------------------------------------------------------------------# #------------------------------------------------------------------# # make .grd file #------------------------------------------------------------------# echo "Gridding file..." xyz2grd temp.f2 -R${xmin}/${xmax}/${ymin}/${ymax} \ -I${cell}/${cell} -V -G${ofile}.grd #------------------------------------------------------------------#

•

At this point we can generate a plot using grdimage to view our elevations. following is an example that will get us in the same vicinity as our image above.

The

#!/bin/csh makecpt -CColor_DEM.cpt -Z -T5000/10500/100 >! mtolympus.cpt grdimage Mt_Olympus.grd -R430000/459995/4488500/4502000 -JX6i/2.7i \ -B10000g100000/2000g1000000000nSeW -Cmtolympus.cpt -V -P -K \ -X1 -Y2 >! olympus.ps pscoast -R -JX3d/2d -Df -V -P -O -K -W2 -N2 >> olympus.ps gs -sDEVICE=x11 olympus.ps

This gives us the following image:

This has significantly more detail than our image generated using Etopo1. But, we can make it look better than this.

Geophysical Computing

L07-5

3. Intensity Files If we really want our elevation images to pop then we need to add some illumination. The easiest way to do this is with the grdgradient command. Grdgradient is used to calculate the directional derivative (or gradient) of a .grd file. Consider our example above where we generated a grid file for the...