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.

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

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_da

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

主题: XMLPythonTI
分页:12
转载请注明
本文标题:Topographic and tectonic stress inversions with halfspace: Preparing a DEM
本站链接:http://www.codesec.net/view/481062.html
分享请点击:


1.凡CodeSecTeam转载的文章,均出自其它媒体或其他官网介绍,目的在于传递更多的信息,并不代表本站赞同其观点和其真实性负责;
2.转载的文章仅代表原创作者观点,与本站无关。其原创性以及文中陈述文字和内容未经本站证实,本站对该文以及其中全部或者部分内容、文字的真实性、完整性、及时性,不作出任何保证或承若;
3.如本站转载稿涉及版权等问题,请作者及时联系本站,我们会及时处理。
登录后可拥有收藏文章、关注作者等权限...
技术大类 技术大类 | 开发(python) | 评论(0) | 阅读(33)