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 )

This is a brief step that simply involves reading and flipping a DEM. DEMs in GeoTIFF format are referenced differently, with the origin at the upper left. We will also get some of the projection metadata from the DEM.


# inline plotting %matplotlib inline %config InlineBackend.figure_format = 'retina'


from __future__ import division import gdal import numpy as np import matplotlib.pyplot as plt import json import pyproj as pj


# 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.




(398285.44218398223, 445.9049868175545, 0.0, 3329810.0066404096, 0.0, -445.9049868175547)


dem_dataset.RasterXSize, dem_dataset.RasterYSize


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

本文开发(python)相关术语:python基础教程 python多线程 web开发工程师 软件开发工程师 软件开发流程

主题: XMLPythonTI
本文标题:Topographic and tectonic stress inversions with halfspace: Preparing a DEM

技术大类 技术大类 | 开发(python) | 评论(0) | 阅读(17)