Crust1.0 Data

Using Python to access Crust 1.0 data files

Gabi Laske, Zhitu Ma, Guy Masters and Michael Pasyanos at Lawrence Livermore National Lab have created a crustal model for the Earth at 1.0 intervals. See this website for more information https://igppweb.ucsd.edu/~gabi/crust1.html and to download the crustal data files.

The Crust1.0 website provides only FORTRAN code to access the data. This notebook presents a method to access the data for crustal density and dimension using Python. A similar approach could be easily adapted to access the compressional and shear wave velocities.

import pandas as pd

# Read in the dateset from csv

import pandas as pd

vp = pd.read_csv('/Users/erhepp/fortran/crust1.0/crust1.vp', delim_whitespace=True, header=None, \
                 names=['top_water', 'bottom_water', 'bottom_ice', 'bottom_sed1', 'bottom_sed2', 'bottom_sed3', \
                       'bottom_crust1', 'bottom_crust2', 'bottom_crust3'])

vs = pd.read_csv('/Users/erhepp/fortran/crust1.0/crust1.vs', delim_whitespace=True, header=None, \
                 names=['top_water', 'bottom_water', 'bottom_ice', 'bottom_sed1', 'bottom_sed2', 'bottom_sed3', \
                       'bottom_crust1', 'bottom_crust2', 'bottom_crust3'])

rho = pd.read_csv('/Users/erhepp/fortran/crust1.0/crust1.rho', delim_whitespace=True, header=None, \
                 names=['top_water', 'bottom_water', 'bottom_ice', 'bottom_sed1', 'bottom_sed2', 'bottom_sed3', \
                       'bottom_crust1', 'bottom_crust2', 'bottom_crust3'])

bnds = pd.read_csv('/Users/erhepp/fortran/crust1.0/crust1.bnds', delim_whitespace=True, header=None, \
                 names=['top_water', 'bottom_water', 'bottom_ice', 'bottom_sed1', 'bottom_sed2', 'bottom_sed3', \
                       'bottom_crust1', 'bottom_crust2', 'bottom_crust3'])

print vp.shape, vs.shape, rho.shape, bnds.shape

(64800, 9) (64800, 9) (64800, 9) (64800, 9)

rho.head()

	top_water	bottom_water	bottom_ice	bottom_sed3	bottom_crust1	bottom_crust2	bottom_crust3
0	1.02	0.92	1.93	2.55	2.85	3.05	3.34
1	1.02	0.92	1.93	2.55	2.85	3.05	3.34
2	1.02	0.92	1.93	2.55	2.85	3.05	3.34
3	1.02	0.92	1.93	2.55	2.85	3.05	3.34
4	1.02	0.92	1.93	2.55	2.85	3.05	3.34

bnds.head()

	bottom_water	bottom_ice	bottom_sed1	bottom_sed2	bottom_sed3	bottom_crust1	bottom_crust2	bottom_crust3
0	-3.69	-3.69	-4.99	-4.99	-4.99	-5.67	-7.15	-11.75
1	-3.66	-3.66	-4.96	-4.96	-4.96	-5.64	-7.13	-11.74
2	-3.64	-3.64	-4.94	-4.94	-4.94	-5.62	-7.11	-11.72
3	-3.62	-3.62	-4.92	-4.92	-4.92	-5.60	-7.09	-11.71
4	-3.61	-3.61	-4.91	-4.91	-4.91	-5.59	-7.08	-11.70

# Return the closest crustal model grid point values for any given set of latitude longitude coordinates. 
# Parameters are latitude, longitide and either rho or bnds, to access density or dimension, respectivey

def get_crust (qlat, qlong, crust_param):
    # Prevent any float representation uncertainty from pushing latitude and longitude out of range
    if qlat >= 89.999: qlat = 89.999
    if qlat <= -89.999: qlat = -89.999
    if qlong >= 179.999: qlong = 179.999    
    if qlong <= -179.999: qlong = -179.999 
    
    # Find the closest crust1.0 grid locatio to the requested latitude and longitude
    if qlat >= 0:
        rlat = round(qlat - 0.5) + 0.5
    else:
        rlat = round(qlat + 0.5) - 0.5
   
    if qlong >= 0:
        rlong = round(qlong - 0.5) + 0.5
    else:
        rlong = round(qlong + 0.5) - 0.5

    # Determine the index in the crust1.0 files that corresponds to the closest grid location  
    # Information from the crust1.0 readme file
    #   The model is defined from 89.5 to -89.5 deg latitude and -179.5 to 179.5 deg
    #   longitude. Longitudes are the inner loop, i.e. all longitudes are stored
    #   for each latitude, then the next latitude is given. The model starts at
    #   89.5 N and 179.5 W.
    
    dist_from_north = 89.5 - rlat
    dist_from_west = 179.5 + rlong
    idx = int((dist_from_north)*360 + dist_from_west)
    
    # Calculate the latitude and longitude of the crust1.0 datapoint, and return the index loacation
    # latitude, longitude and, for each crustal layer, the density data.
    # Information from the crust1.0 readme file on order and meaning of values
    #    1) top of water
    #    2) bottom of water
    #    3) bottom of ice
    #    4) bottom of sediments 1
    #    5) bottom of sediments 2
    #    6) bottom of sediments 3
    #    7) bottom of cryst. crust 1
    #    8) bottom of cryst. crust 2
    #    9) bottom of cryst. crust 3 = Moho (depth to Moho, not crustal thickness!)

    # For this project, we will only use the density values for sediment and crystalline crust layers
    
    flat = 90 - dist_from_north-0.5
    flon = -180. + dist_from_west+0.5
    if crust_param == 'rho':
        return [idx, flat, flon, round(rho.iloc[idx,0],2), round(rho.iloc[idx,1],2), \
                                             round(rho.iloc[idx,2],2), round(rho.iloc[idx,3],2), \
                                             round(rho.iloc[idx,4],2), round(rho.iloc[idx,5],2), \
                                             round(rho.iloc[idx,6],2), round(rho.iloc[idx,7],2), \
                                             round(rho.iloc[idx,8],2)]
    if crust_param == 'bnds':
        return [idx, flat, flon, round(rho.iloc[idx,0],2), round(bnds.iloc[idx,1],2), \
                                             round(bnds.iloc[idx,2],2), round(bnds.iloc[idx,3],2), \
                                             round(bnds.iloc[idx,4],2), round(bnds.iloc[idx,5],2), \
                                             round(bnds.iloc[idx,6],2), round(bnds.iloc[idx,7],2), \
                                             round(bnds.iloc[idx,8],2)]

    
    
    
    
print get_crust(43.4,-110.7,'rho') 

print get_crust(43.4,-110.7,'bnds')

[16629, 43.5, -110.5, 1.02, 0.92, 2.11, 2.46, 0.0, 2.74, 2.78, 2.95, 3.28]
[16629, 43.5, -110.5, 1.02, 2.47, 2.47, 1.47, 0.97, 0.97, -17.25, -33.4, -40.43]

Written on October 3, 2017