如何用Python模拟ArcGIS做坡度分析

———前面的话———

吴道长在算法的道路上越走越远……

有算法癖的同学，可以找吴道长好好聊唠嗑。

道长公众号的二维码：

———我是分割线———

大学时，总结老板所教的栅格分析方式，无非就是两种：叠加分析和邻域分析。

叠加分析，就是N个栅格数据在同一个像素位置上进行数学运算，获得新的栅格，这种方式就是上一篇文章所提到的NDVI计算，通过对影像的第三第四波段进行归一化而得到新的二值化的栅格。

而领域分析就是只在一个栅格上做文章，通过计算和周围像素的相关，来确定当前栅格的。

新值，比较典型就是九宫格模式，计算中心点周围八个方向的像素来给中心点新值，这里典型就是基于规则网格的坡度坡向计算。

算法原理参照是ESRI的帮助文档:

根据算法原理，最后得到e点的坡度值，用程序来实现：

  1. 模拟了一个九宫格一样的数据驱动，当识别到当前的像素坐标          时，给出一当前这个像素周围八个方向上的像素坐标索引；

  2. 将要计算的栅格用GDAL化作一个矩阵，在遍历行列上通过数据驱       动进行识别定位计算，这边我取巧了，九宫格计算在最边缘行或       列并不适用，我选择忽略。

将结果呈现方式，这次想做：结果数据和原始数据对比呈现。这边我写了两套代码，因为空间数据在Matplotlib展示上，市面上有两种常用的包，Cartopy和Basemap。之前我更常用的是Cartopy，他的处理方式是把包安装后，预先运行一段样例代码，该代码通过在开源的数据网站下载SHP数据在本地后进行数据展示，而我们需要就是该本地的路径，而后的操作就是我们将展示的数据直接存放在该路径下，直接用代码识别就行，操作比较简单。对于空间数据展示，核心就是投影的方式，Cartopy是将方法直接嵌入在Matplotlib中subplot()或者axes()中。

将做好的效果做出来，我并不太满意，因为第二张图的Y轴刻度和第一张图绞在一起了。如果去掉Y轴刻度，那么图中的经纬网就会没有纬线，显得不伦不类。这就丧失了我再做进一步美化的欲望。

而后开始尝试用Basemap，相对于Cartopy，Basemap功能是更强大，从下载的whl包就能知道Cartopy才8M，而Basemap达到了119M，然而Basemap基础包Matplotlib也才124M，Basemap直接给自己一个BASEMAP()的方法，在此制订了四至及投影，再此基础上直接用Matplotlib，因此Basemap显然和Matplotlib更加紧密。用Basemap确实能做出我想要的效果。

在这段代码中，labels=[0,0,0,0]，以0/1组合方式来设置刻度的位置及是否存在，相当方便。Basemap也提供多种栅格显示方式，这里我原始数据采用pcolormesh()，渲染的效果会比较柔和，而结果数据我则用contourf()，让渲染更尖锐一点，然而貌似效果也不很明显，是因为结果的极值差太小。

最后我装了Jupyter Notebook，因为没经验，胡乱用了Anaconda的，最后pip就装了，瞅瞅编译的环境效果，还真是舒服。

2017/6/8