暂无图片
暂无图片
暂无图片
暂无图片
暂无图片
首页 / 使用Python进行栅格数据处理系列推文（1）：栅格重分类小技巧

使用Python进行栅格数据处理系列推文（1）：栅格重分类小技巧

我得学城 2022-07-28
571
得先点我得学城关注我哦～
导  语
本文介绍的是在GIS中常用的栅格重分类功能，用于栅格数据的处理，如在双评价或构建最小阻力模型中对评价指标进行重分类。将诸如高程、GDP、人口、NDVI等连续变量转化为离散变量，或是在土地利用分类数据中将二级地类重新划分为一级地类等多种场景。
 铝合金 | 供稿
Sherry | 编辑
 
01

前置准备

使用本文介绍的方法需要rasterio和numpy，前者为读取栅格数据，后者为重分类处理，建议使用conda建立独立的虚拟环境进行相关操作。
下附conda建立环境命令
conda create -n pygeo38 python=3.8 -y && conda activate pygeo38 && conda install -c conda-forge geopandas rasterio && python -m ipykernel install --user --name pygeo38 --display-name "GeoPython 3.8"
成功后进入该环境进行后续实践。
conda activate pygeo38

02
连续变量重分类——重庆市高程重分类
2.1 导入相关库
    from pathlib import Path # 处理文件路径
    import rasterio as rio
    from rasterio.plot import show
    import numpy as np
    import matplotlib.pyplot as plt
      data_path = Path(r'.\data')
      dem_file = data_path  "dem.tif"
        with rio.open(dem_file) as f:
            nodata = f.nodata
            meta = f.meta
            data = f.read(1)
            data[data == nodata] = np.nan
            show(f, cmap='terrain')
        ‍

        2.2 重庆市高程直方图
        当前的数据是格网形式的二维数据，制作直方图需将其压平为一维数据，使用flatten方法将numpy的数组压平。
          data_flat = data.flatten()
          plt.hist(data.flatten(), 20);
          重庆市高程如上图所示，我们以300、600、900、1200、1500、1800、2100、2400 m作为断点将高程数据重分类，划分为9个区间。
          2.3 重分类断点构建
            breaks = np.arange(300, 2401, 300)
              array([ 300,  600,  900, 1200, 1500, 1800, 2100, 2400])
              2.4 条件判断构建方法
              简单条件判断
              先手动写一个简单的条件判断，使用np.where对数组进行分类赋值，第一个参数为条件，第二个为条件为真时的返回值，第三个参数为条件为假的返回值。
              先做个小例子：
              建立一个1到5的数组，当数组内为偶数时（数值除以2，余数为0），数值本身乘-1，当为奇数时，数值本身乘2。
                array1 = np.array([1, 2, 3, 4, 5])
                np.where(array1%2 == 0, array1*-1, array1**2)
                  array([ 1, -2,  9, -4, 25])
                  嵌套条件判断
                  接下来我们做一个嵌套的条件判断。
                  先做个小例子：
                  建立一个1到10的数组，当数值小于等于3，赋值为1，大于3小于等于8赋值为2， 大于8赋值为3。
                    array1 = np.arange(1, 11)
                    np.where(array1 <= 3, 1, 
                       np.where(array1 <= 8, 2, 3) 
                    )
                      array([1, 1, 1, 2, 2, 2, 2, 2, 3, 3])
                      2.5 栅格数据重分类
                      上面的例子的逻辑为先判断是否小于3，如果为真，赋值为1，为假则进行下一步的判断，即是否小于8，以此类推。。‍‍‍‍‍‍‍‍
                      将这个思路应用于我们的使用场景，首先将空值赋值为0，当高程小于300时，赋值为1，300-600赋值为2，以此类推。
                        data_reclass = np.where(np.isnan(data), 0, 
                            np.where(data<300, 1, 
                            np.where(data<600, 2, 
                            np.where(data<900, 3, 
                            np.where(data<1200, 4, 
                            np.where(data<1500, 5, 
                            np.where(data<1800, 6, 
                            np.where(data<2100, 7, 
                            np.where(data<2400, 8, 9)
                        )))))))).astype(np.int8)
                        两种其他形式的写法
                        当然上述的代码写起来颇为冗长，这里也提供两个其他的简短写法，供大家参考。
                          breaks = np.arange(300, 2401, 300) # 不包含最小值
                          array_name = "data"
                          data_reclass = eval(f"np.where(np.isnan({array_name}), np.nan, \n" + "\n".join([f"np.where({array_name} < {b}, {v}, " for v, b in enumerate(breaks, start=1)]) + f"{len(breaks)+1}" + ")" * (len(breaks)+1))
                            breaks = np.arange(0, 2401, 300) # 包含最小值
                            data_reclass = np.sum(np.dstack([(data >= b) for b in breaks]), axis=2)
                            2.6 重分类结果一览
                              x_, y_ =np.unique(data_reclass, return_counts=True)
                              plt.bar(x_[1:], y_[1:])
                                f, axes = plt.subplots(1, 2, figsize=(20, 10))
                                show(data, transform=meta["transform"], cmap='terrain', ax=axes[0])
                                show(data_reclass, transform=meta["transform"], cmap='terrain', ax=axes[1])
                                for ax in axes:
                                    ax.set_xticklabels([])
                                    ax.set_yticklabels([])


                                plt.subplots_adjust(wspace=0, hspace=0)
                                 
                                 python栅格数据分类数据data数据处理
                                   
                                文章转载自我得学城，如果涉嫌侵权，请发送邮件至：contact@modb.pro进行举报，并提供相关证据，一经查实，墨天轮将立刻删除相关内容。
                                   
                                评论