1、第三章 栅格空间分析3.1 栅格分析基础3.2 栅格计算与重分类3.3 距离制图3.4 密度制图3.5 表面分析3.6 邻域分析3.7 统计分析3.8 水文分析CDIO教学改革教案3.1 栅格分析基础 关于栅格数据实际地面栅格图像3.1 栅格分析基础 栅格数据的空间位置表示2222244444555554444455555LeftRightTopbottomY cellsizeX cellsize3.1 栅格分析基础 ArcGIS的Grid离散型Grid连续型Grid3.1 栅格分析基础 栅格数据的Nodata3.1 栅格分析基础 栅格分析环境输出坐标系统输出空间范围CellSize和掩膜3.
2、1 栅格分析基础 掩膜原理分析范围有值Nodata掩膜3.2 栅格计算与重分类 问题：提供区域TM遥感数据，获取下表植被指数植被指数面积面积0-0.2XXXX0.2-0.4XXXX0.4-1XXXX3.2.1 栅格计算 原理 实现4 3 2 2 2 1 2 1 3 3 2 4 4 2 1 2 2 3 4 3 2 1 1 2 4 1 2 0 y=f(x1,x2,.)=x1+x2x1x25 5 5 6 5 4 2 4 5 6 5 6 2 3.2.2 重分类 原理 实现344614242647213407172428391612423824953319316231746301628278152035
3、0-10：110-20：220-30：330-40：Nodata40-50：552335311233225311213252331233.2.3 问题解决 流程植被指数植被指数面积面积0-0.2XXXX0.2-0.4XXXX0.4-1XXXX3.2.3 问题解决 空间分析建模 关键技术：Add Data Or Tool、Make Parameters模型保存到GeoDatabase后使用3.2.3 问题解决 问题扩展：只统计县域范围内的植被情况*3.2.4 基于GDAL的栅格数据操作 参见PPT：基于开源GIS库和C#3.3 距离制图 问题1：城市所有位置上的商业繁华度指标商业繁华度指标：Fi
4、=fi（1-di/d）；F=max（Fi）其中：fi为某级商业中心功能分；d为某级商业中心的服务半径；di为格网中心到该设施直线距离；F为所有商业中心对该格网作用分最大值。现需要计算一级商业中心（fi取100；d取3500）对所有格网的F值？3.3 距离制图 问题2：假设某地新建了一个小学。考虑到学生上学放学交通问题，打算修建一条公路与已有公路连接起来。在修建公路的时候需要考虑地形坡度问题和公路经过地区的土地利用情况，从而减少成本和工程难度公路起点学校拟修建公路3.3 距离制图 3.3.1 直线距离 3.3.2 成本距离加权 3.3.3 最短路径 3.3.4 区域分配3.3.1 直线距离原理：
5、计算每一个格网中心到最近设施(矢量或栅格)的直线距离，以该距离作为该栅格的Cell Value探索：如果设置离散型Grid的分辨率为10米，以图中两个点为设施，A、B栅格的值分别是？02214 02822202228362214101422322010 010202222141014101428222010 0103632221410143.3.1 直线距离A03.3.1 直线距离 算法实现2822202228362214101422322010 010202222141014101428222010 010363222141014A：在计算每个栅格时，遍历每个设施基本上不能实现（速度）B：在
6、计算每个栅格时，逐步扩大的半径搜索设施勉强可以C：先计算与设施相交栅格值，再以栅格扩散最优秀算法3.3.1 直线距离 实现：Euclidean Distance依据选择集3.3.1 直线距离 应用举例：对“省会城市”图层做0-100公里的栅格缓冲区分析Euclidean DistanceRaster Calculator或Reclassify3.3.1 直线距离 应用举例2：对“省会城市地理坐标”图层做0-100公里的栅格缓冲区分析设置显示坐标系统为投影坐标系统设置栅格分析坐标系统Same As Display,以及设置栅格分析范围和CellSize图层按显示坐标系统Export DataEu
7、clidean DistanceReclassify3.3.1 直线距离 问题1 解决设置显示坐标系统地理坐标数据依据显示坐标系统Export DataEuclidean Distance设置栅格分析环境（空间范围、CellSize等）RasterCalculator3.3.2 成本距离加权 问题来源沼泽3.3.2 成本距离加权 原理3.3.2 成本距离加权 实现：Path Distance 或 Cost Distance 关于backlink raster:方向矩阵3.3.3 最短路径 关于最短路径公路起点学校拟修建公路公路起点基本过程：学校成本GridCost Distance路径代价方向
8、矩阵Cost Path输入产生3.3.3 最短路径 关于成本GRID：坡度*0.6+土地利用*0.4Agriculture 4Barren land 6Brush/transitional 5Built up 9Forest 8Water 10Wetland Nodata0-515-10210-15 315-2562510坡度土地利用3.3.3 最短路径 问题2解决3.3.4 区域分配 原理：确定所有格网中心的最近设施，对于矢量格式的设施，以某一整型字段值或FID作为Cell Value；对于离散型栅格的设施，以像元值或整型字段值作为Cell ValueshapeIDPoint1Point23
9、.3.4 区域分配 实现：Cost Allocation Path Distance Allocation Euclidean Allocation3.3.4 区域分配 应用案例:导入Bank.shp,得到某银行（Banks_ID=559）在服务区域（最近直线距离）内拥有的顾客数Euclidean Allocation3.4 密度制图 问题：商业连锁店选址 人口密度必须在每平方公里1000人以上 据最近商业点直线距离在500米以外3.4 密度制图 功能描述300500100400A300w1+100w2+500w3除以搜索半径的圆面积距离近权重大权重相等3.4 密度制图 密度制图与空间插值的差
10、异Why map density?Density surfaces show where point or line features are concentrated.For example,you might have a point value for each town representing the total number of people in the town,but you want to learn more about the spread of population over the region.Since all the people in each town
11、do not live at the population point,by calculating density,you can create a surface showing the predicted distribution of the population throughout the landscape.When added together,the population values of the cells equal the sum of the population of the original point layer.3.4 密度制图 问题解决密度功能Raster
12、 Calculator直线距离人口统计数据商业网点数据3.5 表面分析 3.5.1 空间插值 3.5.2 表面分析功能3.5.1 空间插值 空间插值：由已知点数据来估算其他点数据 空间插值分类 全局插值、局部拟合 精确插值、非精确插值 确定性插值、随机性插值精确插值在样点处估算值与实际值相同确定性插值不提供预测值的误差检验，随机性插值则用估计变异提供估计误差的评价3.5.1 空间插值 空间插值分类 按插值计算方法：分块内插、部分内插和单点移面内插3.5.1 空间插值 ArcGIS空间插值方法插值方法插值方法描述描述反距离加权局部拟合、确定性、精确、单点移面泰森多边形局部拟合、确定性、精确、部分
13、内插自然邻点局部拟合、确定性、精确、部分内插趋势面全局拟合、确定性、非精确克里格局部拟合、随机性、精确、单点移面样条函数局部拟合、确定性、精确、分块内插3.5.1 空间插值 1）反距离加权（IDW）S0S1S2S3S4S53.5.1 空间插值 2）泰森多边形插值（最近邻点插值）使用最邻点值填充栅格像元值Euclidean Allocation（依据FID）待插值矢量点Raster To Polygon泰森多边形Join（属性表连接，分别依据GRIDCODE和FID）Polygon To Raster3.5.1 空间插值 3）自然邻点插值The natural neighbors of any
14、point are those associated with neighboring Voronoi(Thiessen)polygons.Initially,a Voronoi diagram is constructed of all the given points,represented by the olive-colored polygons.A new Voronoi polygon,beige color,is then created around the interpolation point(red star).The proportion of overlap betw
15、een this new polygon and the initial polygons is then used as the weights.3.5.1 空间插值 4）趋势面插值 ArcGIS的Trend工具用于全局趋势面插值，使用n阶（1至13）多项式方程模拟光滑曲面1阶：Z=b0+b1*x+b2*y2阶3.5.1 空间插值 5）克里格插值 不仅考虑临近点的距离，还考虑临近点的空间分布h变程拱高块金常数块金常数是由于观测误差造成的，代表随机性部分。越小表示空间相关性越强拱高+块金常数=基台值，表示数值最大变化幅度。当基台值一定时，拱高越大，表示结构性变化比重越大变程：空间相关性变化的范
16、围3.5.1 空间插值 5）克里格插值 实际的变异函数很难得到，常使用模型（球面、指数、线性、高斯等）3.5.1 空间插值 5）克里格插值 步骤：样点空间结构量化分析（半变异函数分析）；利用半变异函数、样点空间分布及样点数据值对未知点预测 分类 普通克里格：区域化变量数学期望E(Z(s)为常量 简单克里格：假定Z(s)=u+(s)，其中u为常数，(s)为随机性误差 范克里格：Z(s)=u+(s)，u不为常数5）克里格插值 普通克里格插值法举例（1）半变异函数分析，使用球状拟合模型邻域点编号邻域点编号X坐标坐标Y坐标坐标观测值观测值113105215100334105445100551115求待
17、估计点（2.75，2.75）的估计值5）克里格插值 普通克里格插值法举例（2）计算观测点距离样点对样点对距离距离样点对样点对距离距离1,22.0002,43.0001,32.2362,55.6571,43.6053,41.4141,54.4723,53.6062,32.2364,54.1245）克里格插值 普通克里格插值法举例（3）建立伽马矩阵如计算矩阵中样点1和样点3的半变异函数值：5）克里格插值 普通克里格插值法举例（4）求解伽马矩阵的逆矩阵5）克里格插值 普通克里格插值法举例（5）求解预测点与观测点的拟合半变异g从预测点到观测点从预测点到观测点距离距离11.76822.85031.275
18、42.57452.850从预测点到观测点从预测点到观测点拟合半拟合半变异变异132.097249.936323.390445.584549.93661.000计算同前面的伽马矩阵各元素5）克里格插值 普通克里格插值法举例（6）计算样点权系数：伽马逆矩阵*g（7）按样点权计算估计值以样点1为例：5）克里格插值 ArcGIS的克里格插值3.5.1 空间插值 6）Spline样条函数 重要参数：曲面类型（REGULARIZED、TENSION）；权重（决定光滑度或粗糙度）；插值点数按一定规则，用相邻数据点连线将块分割成若干个多边形分片分片，通过每分片上的全部数据点，展铺一张光滑的数学曲面数学曲面，并
19、使相邻分片间保持连续光滑的拼接。(连续分布)3.5.2 表面分析功能 问题1：漂亮的底图 问题2：退耕还林还草决策支持分析，导入LandUse和Elevation，得到分析区内坡度大于25度且Landuse为Agriculture的区域3.5.2 表面分析功能 1）地形因子的提取 2）等值线 3）山体阴影 4）应用举例1）地形因子的提取 坡向：斜坡的面对的方向。给定点位的坡向定义为斜面法线在平面上的投影和正北方向之间夹角1）地形因子的提取 坡度：指斜坡的倾斜程度。给定点位的坡度定义为斜面法线和垂直方向之间夹角。1）地形因子的提取 坡度和坡向Horn算法原理e1e2e3e4C0e5e6e7e8d
20、坡度S、坡向D计算公式：1）地形因子的提取 地面曲率、平面曲率、剖面曲率 地形表面上一点的扭曲变化程度的度量。在垂直和水平方向分量分别称为剖面曲率和平面曲率 剖面曲率是对地面坡度的沿最大坡降方向地面高程变化率的度量 平面曲率指在地形表面上，具体到任何一点，指用过该点的水平面沿水平方向切地形表面所得的曲线在该点的曲率值。也就是该点所在的地面等高线的弯曲程度1）地形因子的提取 地面曲率、平面曲率、剖面曲率地面曲率剖面曲率平面曲率2）等值线 算法思想 实现：Contour 和Contour List等值线间距初始等高线3）山体阴影 原理：给定太阳高度角和方位角，模拟地表光照值为避免产生反立体，请勿使
21、用偏南方向的太阳光3）山体阴影 计算方法：山体阴影的应用举例：图层融合H=255.0*(cos(Az)*cos(S)+sin(Az)*sin(S)*cos(Al-As)其中：Az为太阳入射角，即太阳高度角的余角；S为格网的坡度；Al为入射方向，等于2*Pi-太阳方位角+Pi/2图层透明度50%4）应用举例 问题2解决：退耕还林还草决策支持分析，导入LandUse和Elevation，得到分析区内坡度大于25度且Landuse为Agriculture的区域坡度分析ElevationLandUseRaster Calculator坡度数据3.6 邻域分析 问题1：提供区域建筑面矢量图层，模拟计算城
22、市任一点的建筑密度 问题2：求取满足或不满足在500米范围内至少有5个顾客点的银行3.6 邻域分析 原理：以计算像元为中心的窗口范围内计算窗口统计值，并以该统计值作为中心像元的分析结果值 统计值的类型：Mean，Maximum、Minimun、Median、Sum、Range、Majority、Variety、Standard Dev3.6 邻域分析 栅格邻域分析：Focal Statistics23 23 23 2825 25 25 2828 28 28 2830 28 28 23计算栅格图层新栅格图层统计最大值3.6 邻域分析 矢量点邻域分析：与栅格数据类似，差异在于只统计窗口内矢量点的某
23、字段值A300500100400如统计Sum，则A栅格值为9003.6 邻域分析 栅格邻域分析举例：平滑等值线均值滤波DEM平滑DEM等值线平滑等值线3.6 邻域分析 问题1解决：建筑密度的计算归一化居民地：建筑内像元为1重分类居民地：建筑外像元为0背景栅格：所有像元为1领域栅格：矩形窗口21*21，统计Sum建筑密度=居民地邻域/背景领域*1003.6 邻域分析 问题2解决投影后的顾客图层设置显示坐标系统为投影坐标系统银行和顾客点图层按显示坐标系统Export Data投影后的银行图层Point Statistics(500米;Circle)增加ID字段，取值全为1顾客邻域Raster To
24、 Polygon顾客领域多边形3.7 统计分析 3.7.1 单元统计 3.7.2 分区统计 分区统计几何特征 分区汇总统计 分区直方图统计 交叉列表统计3.7.1 单元统计 原理：同区域Grid图层同位置做统计3.7.1 单元统计 应用举例 1）植被指数十天合成 2）商业繁华度指标计算最终实现商业繁华度指标：Fi=fi（1-di/d）F=max（Fi）其中：fi为某级商业中心功能分d为某级商业中心的服务半径di为格网中心到该设施直线距离F为所有商业中心对该格网作用分最大值级别fid1100350026520003.7.1 单元统计 商业繁华度指标计算最终实现投影后的商业点图层设置显示坐标系统为
25、投影坐标系统商业点图层按显示坐标系统Export Data选择1级商业中心Euclidean DistanceRaster Calculator选择2级商业中心Euclidean DistanceRaster Calculator1级分值2级分值Cell Statistics:最大值3.7.2 分区统计 1）分区统计几何特征 Zonal Geometry Zonal Geometry As TableValueAreaPerimeter.05.014.15.014.23.08.42.06.CellSize=1.03.7.2 分区统计 2）分区汇总统计 原理：以一个离散型Grid作为分区，统计同
26、区域另一个Grid每个区域的Area、Mean等分区离散型Grid被统计的Grid.51/106821000102.9/623160061.MeanRangeMaxMinAreaCountValue3.7.2 分区统计 2）分区汇总统计 实现：Zonal Statistics、Zonal Statistics As TableZonal Statistics：最小值3.7.2 分区统计 2）分区汇总统计 应用举例：统计如下坡度范围内农用地的面积0-5度度5-15度度15-25度度 25度度面积DEMSlope坡度数据Reclassify坡度分级土地利用Zonal Statistics As T
27、able3.7.2 分区统计 2）分区汇总统计 综合应用案例：不同粮食产区（低产区、中产区和高产区）土壤有机质统计粮食产量采样点土壤有机质采样点3.7.2 分区统计 3）分区直方图统计：Zonal Histogram以离散型作为Grid分组，以另一个Grid的当前图例分组创建直方图。3.7.2 分区统计 3）分区直方图统计 应用案例：根据DEM做不同粮食产量区的坡向直方图DEM坡向数据Aspect粮食产量采样点空间插值粮食产量Reclassify粮食产量分区分区直方图分析3.7.2 分区统计 4）交叉列表统计：Tabulate Area以一个离散Grid为行，另一个同区域离散Grid为列，统计
28、交叉面积。3.7.2 分区统计 4）交叉列表统计 应用案例：根据DEM做不同粮食产量区不同坡向的面积粮食产量采样点空间插值粮食产量Reclassify粮食产量分区DEMAspect坡向数据Reclassify坡向分区交叉列表统计3.8 水文分析 引入：唐家山堰塞湖汶川大地震后形成的最大堰塞湖，地震后山体滑坡，阻塞河道形成的唐家坝堰塞湖位于涧河上游距北川县城约6公里处，是北川灾区面积最大、危险最大的一个堰塞湖。库容为1.45亿立方米。坝体顺河长约803米，横河最大宽约611米，顶部面积约30万平方米，由石头和山坡风化土组成，湖上游集雨面积3550平方公里。3.8 水文分析 3.8.1 水文分析基
29、本功能 3.8.2 堰塞湖水文特征分析 *3.8.3 地形因子的提取3.8.1 水文分析基本功能 水文分析基本功能DEM水流方向矩阵流水累积量水网水流长度计算河网连接填充洼地输 入流域盆地河网分级3.8.1 水文分析基本功能 1）水流方向矩阵 水流方向：水流离开每一栅格单元的指向 水流方向计算：中心栅格与领域网格最大距离落差 水流方向编码与计算787269746756695344水流方向编码326412816x1842+高程矩阵水流方向3.8.1 水文分析基本功能 1）水流方向矩阵3.8.1 水文分析基本功能 2）填充洼地（水文分析起点）目标：每个网格水流都能够流到流域出口 基本算法：填平处理
30、、平地微调增高767269744058695344767269744458695344填平处理3.8.1 水文分析基本功能 3）流水累积量 基本原理：假设研究区域每个格网有一个单位1的水量，按照水流方向累积计算每个格网的上游来水量。水流方向矩阵流水累积量矩阵3.8.1 水文分析基本功能 3）流水累积量ArcGIS支持权重栅格，如果未设置权重栅格，假定各栅格水量为13.8.1 水文分析基本功能 4）水网长度计算Flow Length 原理：根据水流方向Grid，计算每个格网距上游源头或下游出口的流水长度UpStreamDownStream3.8.1 水文分析基本功能 5）流域盆地 由分水岭分割而
31、成的汇水区域。它通过对水流方向数据的分析确定出所有相互连接并处于同一流域盆地的栅格3.8.1 水文分析基本功能 6）水网提取（山谷线提取）流水累积量Reclassify栅格河网（非河网栅格为Nodata）水流方向矩阵水网提取矢量河网探索：区域山脊线如何提取？3.8.1 水文分析基本功能 7）河网连接与分级Stream LinkStream Order流水累积量Reclassify栅格河网（非河网栅格为Nodata）水流方向矩阵水网提取河网连接河网分级3.8.2 堰塞湖水文特征分析 捕捉出水口Snap Pour Point 集水流域WaterShed源点捕捉的出水口捕捉半径Snap Pour P
32、ointWaterShed3.8.2 堰塞湖水文特征分析 堰塞湖集水区域计算集水流域栅格作为掩膜，对原DEM实施Raster Calculator,获取高程小于等于水位线的所有栅格3.8.3 地形因子的提取本节讲述坡长提取、沟壑密度提取3.8.3 地形因子的提取 1、沟壑密度 在特定区域内，地表单位面积沟壑的总长度水网提取计算河流长度沟壑总长/总面积3.8.3 地形因子的提取 2、坡长的提出 坡长：山脊水流沿坡向下水流长度Raster CalculatorDEM负地形水网水文分析思考题 1、地形起伏度是指，在所指定的分析区域内所有栅格中最大高程与最小高程的差。以3*3矩形窗口进行计算，给出计算流程。2、提供城市建筑多边形图层和DEM图层，请分析城市建筑密度与坡度和坡向的关系 3、给出堰塞湖水域面积、库容计算的流程
