栅格、矢量数据处理、统计 函数自用

一些自用的处理遥感数据小工具

遥感数据大致分为栅格（TIF）和矢量(shp）数据。TIF数据可以简化描述为辅助信息（info）和图像数组（array）；shp数据可以简化描述为对点集的形态描述。

1.栅格数据

可使用gdal和rasterio第三方库对tif进行操作

1.1 数据读取

tif = rasterio.open(tif_file)
基本信息：
height = tif.height
width = tif.width

1.2 转为数组

读取第一波段为数组，rasterio中读取数组是numpy的array格式
samples = tif.read(1)

1.3 波段计算

以计算ndvi为例：

red = tif.read(3).astype("float64")
nir_red = tif.read(4).astype("float64")
ndvi = (nir_red - red) / (nir_red + red)

栅格数据通常有无数据值填充，计算栅格特征值时需要去除无数据值的影响，如果无数据值是特定值，可以将其改为nan。

def tif_max(tif_arr):
    arr_max = np.nanmax(tif_arr)
    return arr_max


def tif_min(tif_arr):
    arr_min = np.nanmin(tif_arr)
    return arr_min


def tif_mean(tif_arr):
    arr_mean = np.nanmean(tif_arr)
    return arr_mean

替换-9999为nan

a[a == -9999] = np.nan

1.4 拼接tif

使用merge函数，此方法对超大数据无法使用，后期更新该算法。

import rasterio
from rasterio.merge import merge
from rasterio.enums import Resampling
import os

def concatenate_images(input_folder, output_file):
    image_files = [f for f in os.listdir(input_folder) if f.endswith('.tif')]
    if not image_files:
        print("文件夹中没有.tif文件。")
        return

    src_files_to_mosaic = []
    for image_file in image_files:
        src = rasterio.open(os.path.join(input_folder, image_file))
        src_files_to_mosaic.append(src)

    mosaic, out_trans = merge(src_files_to_mosaic, resampling=Resampling.nearest)

    # 更新输出变换矩阵，使其适应新的尺寸
    out_trans = rasterio.transform.Affine.translation(0.0, 0.0) * out_trans * rasterio.transform.Affine.scale(1.0, 1.0)

    # 更新输出元数据
    out_meta = src.meta.copy()
    out_meta.update({"driver": "GTiff",
                     "height": mosaic.shape[1],
                     "width": mosaic.shape[2],
                     "transform": out_trans,
                     "compress": "lzw"})

    # 保存结果图像
    with rasterio.open(output_file, "w", **out_meta) as dest:
        dest.write(mosaic)

# 使用示例
input_folder = ''  # 输入文件夹路径
output_file = '.tif'  # 输出文件路径

concatenate_images(input_folder, output_file)

1.5 统计重分类后占比

遥感影像在重分类后会使用几个像素值分类地物，每种地物的像素数占比等于分类的面积占比。
例如：分类后像素值为1-5，分别统计1-5的的像素值各有多少像素。

def cal_percent(tif_arr, star: int, stop: int, sum):
    # 存放每一个像素值个数的列表
    count_list = []
    for i in range(star, stop + 1):
        count = np.sum(tif_arr == i)
        count_list.append(count)
    count_sum = np.sum(count_list)
    percent_list = [a / count_sum for a in count_list]
    area = [i * sum for i in percent_list]
    return percent_list, area

2 矢量数据

矢量数据通常是后期生成的，使用不同字段存放矢量数据信息，处理也主要集中在对字段的操作。
使用geopandas对shp数据操作。

import geopandas as gpd

def cal_field(shp_path:str, field:str):
    # 读取shp文件
    shapefile = shp_path
    data = gpd.read_file(shapefile)

    # 指定要求和的字段名称
    field_name = field

    # 计算指定字段的总和
    field_sum = data[field_name].sum()
    print(f"字段 {field_name} 的总和为: {field_sum}")
    
    return field_sum