一些自用的处理遥感数据小工具
遥感数据大致分为栅格(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