Topographic and tectonic stress inversions with halfspace: Preparing a DEM

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()