Topographic and tectonic stress inversions with halfspace: Preparing a DEM
Richard Styron
In order to complete the tutorial, please first clone/download the github repo which comes with the DEM, coseismic slip model, and some auxiliary files, all with the correct pathnames.
(Or just download this notebook)
In [1]:
# inline plotting
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
In [2]:
from __future__ import division
import gdal
import numpy as np
import matplotlib.pyplot as plt
import json
import pyproj as pj
Open GIS raster file and get metadata¶
In [3]:
# open dataset in gdal
dem_dataset = gdal.Open('../test_data/dem/baloch_500m_wide_utm41n.tif')
The dem_dataset is a container class, with various methods for reading data and metadata. No data has been read yet.
We will use the GetGeoTransform() method to read in the metadata as a tuple.
In [4]:
dem_dataset.GetGeoTransform()
Out[4]:
(398285.44218398223, 445.9049868175545, 0.0, 3329810.0066404096, 0.0, -445.9049868175547)
In [5]:
dem_dataset.RasterXSize, dem_dataset.RasterYSize
Out[5]:
(1371, 1516)
Write GeoTransform parameters to dictionary¶
Python dictionaries can be saved as JSON text files, which are widely used and also much easier to work with than similar files like XML. Additionally, it is very easy to read JSON files back into Python dictionaries, which makes them ideal for metadata and configuration files.
In [6]:
# write GeoTransform parameters to dictionary
dem_transform = {'upper_left_x': dem_dataset.GetGeoTransform()[0],
'x_res_m' : dem_dataset.GetGeoTransform()[1],
'x_rotation': dem_dataset.GetGeoTransform()[2],
'upper_left_y': dem_dataset.GetGeoTransform()[3],
'y_rotation': dem_dataset.GetGeoTransform()[4],
'y_res_m': dem_dataset.GetGeoTransform()[5],
'n_cols': dem_dataset.RasterXSize,
'n_rows': dem_dataset.RasterYSize}
dem_transform['y_res_m'] *= -1 # correct for the upper-left origin thing
print(dem_transform)
{'y_res_m': 445.9049868175547, 'n_cols': 1371, 'upper_left_x': 398285.44218398223, 'upper_left_y': 3329810.0066404096, 'x_rotation': 0.0, 'x_res_m': 445.9049868175545, 'n_rows': 1516, 'y_rotation': 0.0}
Now let's add some other metadata which will come in handy.
In [7]:
dem_transform['east_min'] = dem_transform['upper_left_x']
dem_transform['east_max'] = (dem_transform['x_res_m'] * dem_transform['n_cols']
+ dem_transform['east_min'])
dem_transform['north_max'] = dem_transform['upper_left_y']
dem_transform['north_min'] = (-1 * dem_transform['y_res_m'] * dem_transform['n_rows']
+ dem_transform['north_max'])
In [8]:
# get latitude and longitude for the corners of the array, just 'cause.
wgs84 = pj.Proj(init='epsg:4326')
utm41n = pj.Proj(init='epsg:32641')
dem_transform['lon_min'], dem_transform['lat_min'] = pj.transform(utm41n, wgs84,
dem_transform['east_min'],
dem_transform['north_min'])
dem_transform['lon_max'], dem_transform['lat_max'] = pj.transform(utm41n, wgs84,
dem_transform['east_max'],
dem_transform['north_max'])
In [9]:
# save metadata dictionary to file
with open('../test_data/dem/baloch_dem_metadata.json', 'w') as f:
json.dump(dem_transform, f)
Read in DEM data¶
In [10]:
# read DEM data as numpy array
dem = dem_dataset.GetRasterBand(1).ReadAsArray()
Let's look at things.¶
First let's make some indices for the array, which is helpful for plotting.
In [11]:
lat_index = np.linspace(dem_transform['lat_min'], dem_transform['lat_max'],
num=dem_transform['n_rows'])
lon_index = np.linspace(dem_transform['lon_min'], dem_transform['lon_max'],
num=dem_transform['n_cols'])
In [12]:
# look at dem
plt.figure(figsize=(14,12))
plt.pcolormesh(lon_index, lat_index, dem)
plt.colorbar()
plt.xlim([dem_transform['lon_min'], dem_transform['lon_max']])
plt.ylim([dem_transform['lat_min'], dem_transform['lat_max']])
plt.xlabel('longitude')
plt.ylabel('latitude')
plt.show()
Hmmm, south is up.¶
In [13]:
# turn that frown upside-down
dem = np.flipud(dem)
In [14]:
# let's look again, with a better color scheme
# look at dem
plt.figure(figsize=(14,12))
plt.pcolormesh(lon_index, lat_index, dem,
cmap=plt.cm.gist_earth, vmin=-1000)
plt.colorbar()
plt.xlim([dem_transform['lon_min'], dem_transform['lon_max']])
plt.ylim([dem_transform['lat_min'], dem_transform['lat_max']])
plt.xlabel('longitude')
plt.ylabel('latitude')
plt.show()
In [15]:
# save numpy array of dem
np.save('../test_data/dem/baloch_dem_utm41n_445m.npy', dem)
Table of Contents for the halfspace tutorial
- Installation and configuration
- Preparing a Digital Elevation Model file
- Preparing the coseismic slip model
- Calculating the topographic stress field
- Calculating the topographic stresses on the fault
- Bayesian inversion of tectonic stress pt. 1: Rake inversion
- Bayesian inversion of tectonic stress pt. 2: Conditions at failure