python 调用gdal 处理dem数据

原作者：Sam Johnstone，斯坦福大学地理学博士 ——翻译:橙子
再见Arcgis，你好Python -第1部分-简介和大纲
原作者官网链接作者官网链接
这是第一系列的文章，我将致力于写一个教程，用python做一些基本的（可能更高级的）DEM处理。因为近期没有arcgis,完成一些任务备受挫折，所以我渴望通过锻炼使我自己渴望更加熟悉Python。我希望其他从事这类工作的人可以利用这一点开发自己的技术，我们可以更容易地拥有更容易、更透明、更开放的分析教程。
……………………….
………………………..
第一步：安装gdal,作者之前是用matlab 处理的，matlab读取dem,matlab中有一个映射工具箱，读取ascii码，至于gdal的概念这里不再多余阐述。
教程目录：
1、 导入dem并投影
2、 创建一个斜坡网格，绘制斜坡
3、 创建一个阴影，阴影覆盖的斜坡
4、 计算网格的曲率值
5、 存数据
6、 计算流域
7、 填坑
（上面的工具就足以制作一些简单的景观演化模型）
8、 模拟坡面
9、 坡面+通道
10、 测量长剖面（纵断面）
11、 数据滤波和光谱分析
可惜的是博士更新到第三章后就消失了，遗憾！！！！！！！！
第一小节：读入dem并投影为utm
笔者在地理空间数据云下载dem,如下:

利用博士的代码，有一部分报错，笔者将其修改后，代码如下，投影后的结果如下：

# -*- coding: utf-8 -*-
from osgeo import gdal
import osr
import numpy as np
from matplotlib import pyplot as plt
from gdalconst import *
from matplotlib import cm
gdal.AllRegister()
filename = r"C:\Users\Y\Desktop\12\ASTGTM2_N00E011\ASTGTM2_N00E011_num.tif"
data = gdal.Open(filename, GA_ReadOnly)
def convertToUTM(dataset, dx, utmZone):
    datum = 'WGS84'
    oldRef = osr.SpatialReference()  #初始化空间坐标

    # 从数据集中复制投影坐标系
    oldRef.ImportFromWkt(dataset.GetProjectionRef()) 
    newRef = osr.SpatialReference()
    newRef.SetUTM(abs(utmZone), utmZone > 0)
    # 创建一个坐标转换对象
    transform = osr.CoordinateTransformation(oldRef, newRef) 
    tVect = dataset.GetGeoTransform()  
    nx, ny = dataset.RasterXSize, dataset.RasterYSize 
    #转换左上角坐标,x,y
    (ulx, uly, ulz ) = transform.TransformPoint(tVect[0], tVect[3])
    #转换左下角坐标
    (lrx, lry, lrz ) = transform.TransformPoint(tVect[0] + tVect[1]*nx+ny*tVect[2],
                                              tVect[3] + tVect[5]*ny+nx*tVect[4])
    #分别是左上角x，每个像素宽，旋转，分别是左上角y,旋转，像素高
    newtVect = (ulx, dx, tVect[2], uly, tVect[4], -dx)
    return newtVect
nx, ny = data.RasterXSize, data.RasterYSize 
dx = 10 #格网间隔
utmZone = 10  #南半球为负,utm 带号=（经度整数位/6）的整数部分+31（东经为正值，西经为负值）

gridNew = data.ReadAsArray().astype(np.float)#所有点高程
(upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = convertToUTM(data, dx, utmZone)
xllcenter = upper_left_x + dx/2  # 左下角像素中心x坐标（utm)
yllcenter = upper_left_y - (ny-1)*dx - dx/2 # 左下角像素中心y坐标
###yllcenter 之所以是<- (ny-1)*dx - dx/2>,是此时坐标系为笛卡尔，
##东西为x,南北朝上为y,左上角坐标减去(ny-1)*dx可得左下角像素坐标，再减去0.5个格网，就得到左下角像素的中心坐标
#每个像素的坐标
xcoordinates = [x*dx + xllcenter for x in range(nx)]
ycoordinates = [y*dx + yllcenter for y in range(ny)]

#创建2d格网描述x,y
X,Y = np.meshgrid(xcoordinates, ycoordinates) 

plt.figure()
plt.contourf(X, Y, gridNew, levels=np.linspace(np.amin(gridNew[gridNew > 0]),
                   np.amax(gridNew), 50),cmap=cm.gist_earth)
plt.show()
data=None

  投影后的

转载自：https://blog.csdn.net/qq_15642411/article/details/79123677