Skip to content

How to make gridded data sets to use with PlateFlex

Jump to bottom
paudetseis edited this page Sep 26, 2019 · 8 revisions

Packaged with this software is a set of topography, gravity and crustal structure data that are easily imported into the PlateFlex package. Those data sets have been produced by processing global grids of data from the WGM2012, GEBCO and CRUST1.0 models using GMT.

  • The WGM2012 models are used to extract Topography, Bouguer and Free-air anomaly data over North America, and Free-air anomaly data over the NW Pacific.
  • The GEBCO model is used for the Bathymetry over the NW Pacific.
  • The CRUST1.0 models are used to extract the crustal thickness and crustal density (corresponding to upper crust, or layer #6) over both North America and the NW Pacific.

Since those models are all defined using a different registration and grid sampling, we use GMT to manipulate the grids and produce consistent data sets to use with PlateFlex. We used GMT 5.4.1 in the following examples. These commands should be copied to a bash script file and executed from the terminal.

WGM2012

These models are defined at a resolution of 2' with a gridline registration.

# Set region and projection for North America. Use a transverse mercator 
# with a central meridian of -95 degrees to minimize high-latitude distortions.
reg=-R-125/13/-28/58+r
proj=-JT-95/3i

# Make map of region
ps=Bouguer_NA.ps
gmt grdimage WGM2012_Bouguer_ponc_2min.grd $reg $proj -K -P \
            -B10WSne -CPALET_WGM_Bouguer_Global.cpt -X1.5i -Y1.25i > $ps
gmt psscale -CPALET_WGM_Bouguer_Global.cpt -DJRM+o0.25i/0+e+mc -R -J -O >> $ps

# Project onto Cartesian grid with sampling interval of 20 km
gmt grdproject WGM2012_Bouguer_ponc_2min.grd -GBouguer_NA.grd $proj $reg -D20 -Fk

# Save grid information as header for use with PlateFlex software
gmt grdinfo -C Bouguer_NA.grd > header.txt
sed 's/^/# /g' header.txt > Bouguer_NA.xyz 
gmt grd2xyz Bouguer_NA.grd >> Bouguer_NA.xyz

You can repeat these commands for the corresponding Free-air anomaly and Topography data.

GEBCO

This model is defined on a very fine grid (15" resolution) with a gridline registration. We first resample the grid onto a 2' grid.

# Set region and projection for NW Pacific - here we use a regular 
# Mercator since it straddles the equator.
reg=-R145/215/-15/45
proj=-JM3i

# Resample on a 2' grid
gmt grdsample GEBCO_data/GEBCO_2019.nc -GGEBCO_resampled_2m.grd -I2m -Rd

# Make map of region
ps=Bathy_PAC.ps
gmt grdimage GEBCO_resampled_2m.grd $reg $proj -K -P \
            -B10WSne -CPALET_WGM_ETOPO1_Global.cpt -X1.5i -Y1.25i > $ps
gmt psscale -CPALET_WGM_ETOPO1_Global.cpt -DJRM+o0.25i/0+e+mc -R -J -O >> $ps

# Project onto Cartesian grid
gmt grdproject GEBCO_resampled_2m.grd -GBathy_PAC.grd $proj $reg -D10 -Fk

# Dump grid onto ASCII file with header
gmt grdinfo -C Bathy_PAC.grd > header.txt
sed 's/^/# /g' header.txt > Bathy_PAC.xyz 
gmt grd2xyz Bathy_PAC.grd >> Bathy_PAC.xyz

CRUST1.0

The CRUST1.0 model is defined on a very coarse grid (1 degree). This is not a problem as we only use those fields in the inversion step, and not in the calculation of the wavelet transform. However, we still require an identical grid specification, so the first step is to resample them on a fine grid using spline interpolation.

Note: You first need to extract the proper layer within CRUST1.0. For crustal thickness, the weblink contains a global grid of crustal thickness (crsthk.xyz). For the density layer, you will have to use the provided fortran code getCN1xyz.f and rename the output file xyz-ro6 to something like crustal_density.xyz (for example). Layer 6 corresponds to the upper crustal layer. See the weblink for details.

# Set region and projection for North America
reg=-R145/215/-15/45
proj=-JM3i

# Create 2 .grd file from the ASCII data
gmt xyz2grd crsthk.xyz -Gcrust_thickness_1deg.grd -Rd -I1 -T

# Resample on a 2' grid using spline interpolation
gmt grdsample crust_thickness_1deg.grd -Gcrust_thickness_2m.grd -nb -Rg -I2m

# Make map of region
ps=Crust_thick_PAC.ps
gmt makecpt -Chaxby -T5/30/1 > crust.cpt
gmt grdimage crust_thickness_2m.grd $reg $proj -K -P \
            -B10WSne -Ccrust.cpt -X1.5i -Y1.25i > $ps
gmt psscale -Ccrust.cpt -DJRM+o0.25i/0+e+mc -R -J -O >> $ps

# Project onto Cartesian grid
gmt grdproject crust_thickness_2m.grd -Gcrust_thickness.grd $proj $reg -D10 -Fk

# Dump grid onto ASCII file with header
gmt grdinfo -C crust_thickness.grd > header.txt
sed 's/^/# /g' header.txt > crustal_thickness_PAC.xyz 
gmt grd2xyz crust_thickness.grd >> crustal_thickness_PAC.xyz
Clone this wiki locally