26、使用NumPy和SciPy进行地图代数运算

使用NumPy和SciPy进行地图代数运算

1. 地图代数概述

地图代数是一种使用代数运算来处理栅格数据集的方法，这些代数运算包括加法、减法等，只不过这里使用的是两个或多个栅格，而非单纯的数字。通过地图代数，可以以各种方式处理栅格数据，从简单到复杂，比如对数据集进行润色使其在地图上显示得更好，或者从一个或多个现有数据集中创建全新的数据集。

地图代数主要有四种类型，适用于不同类型的分析：
- 局部分析 ：对单个像素进行操作。
- 焦点分析 ：使用目标像素周围的像素来计算值。
- 区域操作 ：对具有相同值的像素进行操作。
- 全局分析 ：对整个数组进行操作。

1.1 创建新栅格的函数

在深入细节之前，我们先定义一个函数，用于创建与现有数据集具有相同维度、地理变换和投影的输出GeoTIFF文件。以下是函数代码：

def make_raster(in_ds, fn, data, data_type, nodata=None):
    """Create a one-band GeoTIFF.
    in_ds     - datasource to copy projection and geotransform from
    fn        - path to the file to create
    data      - NumPy array containing dat