栅格转矢量:
栅格转矢量是遥感数据处理中常见且重要的步骤,特别是在需要对特定区域或特定值范围进行分析时。传统的工具,如 ArcGIS 中的“栅格转面”工具,虽然提供了基础的栅格转矢量功能,但灵活性较差。在实际应用中,我们可能只需要将某个阈值或阈值以上的栅格值转为矢量,而不是整个栅格数据都转为矢量,这时 ArcGIS 的转化方式常常会导致操作卡顿,并且生成的矢量数据往往无法满足我们的需求。
相比之下,ENVI 提供了通过设置阈值生成 roi的功能,虽然操作上较为灵活,但导出为 shp 时仍然存在一些麻烦。为了提高效率,我提供了一个更为便捷的栅格转矢量方法,能够根据不同的需求更精准地提取所需区域,避免了冗余数据的生成和处理过程中的性能问题。
栅格分块:
另一方面,当我们面对非常大的栅格数据时,通常会发现直接导入到 ArcGIS 进行可视化时会出现加载缓慢甚至卡顿的现象。因此,分块处理大规模栅格数据成为必不可少的步骤。通过将栅格数据分块,我们可以更加高效地导入、处理和分析这些数据。接下来,我将提供一个简单易用的栅格分块代码,帮助大家快速将大栅格数据分割成更小的部分,便于后续处理和可视化。
1.栅格转矢量
在遥感数据处理过程中,栅格转矢量是常见的操作。通常,我们需要将某个特定值、某个阈值以上或以下的栅格值转换为矢量,以便进行进一步分析。使用 rasterio.features.shapes 函数,我们可以提取栅格数据的几何形状,从而实现栅格转矢量。
为了加快数据处理的效率,我们先将符合条件的栅格值转化为布尔值,再将其转换为整数值,如 1。这种方法不仅能将连续区域的同一值合并为一个矢量要素,还能避免逐像素处理带来的性能问题。如果我们直接对原始影像进行栅格转矢量操作,每个像元都会生成一个独立的矢量要素,尽管每个要素仍然可以包含其对应的栅格值并存储在属性中,但这种方式会显著增加计算和存储。因此,通过先转换为布尔值或整数值再进行矢量化处理,能够提高效率,且保留了原始栅格数据的关键信息。
注意这段代码还可以生成栅格外廓线哦,比如说一个不规则的栅格有效值都是大于0,背景值等于0,那么将mask写成mask = image > 0即可。

import rasterio
from rasterio.features import shapes
import geopandas as gpd
# 读取数据
with rasterio.open(r'C:\Users\HP\Desktop\test011\night2019.tif.tif') as src:
image = src.read(1)
# 创建一个掩膜,DN 值大于阈值的区域为 True
mask = image > 2
# 使用 shapes 函数提取多边形,读取每个矢量shapes() 是读取矢量
results = (
{'properties': {'raster_val': v}, 'geometry': s}
for i, (s, v)
in enumerate(shapes(mask.astype(rasterio.uint8), mask=mask, transform=src.transform)))
#rasterio.features.shapes(source, mask=None, connectivity=4, transform=Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0))
# shapes返回两个值一个是pologon一个是value,函数中介绍value是来自影像,如果这里source填的是image,返回的是各像元转矢量,value就是对应像元的值
# mask.astype(rasterio.uint8)是将该大于阈值的为1,小于阈值的为0,这能将同一值的像元合并,避免单个像元导出为矢量
# 将结果转换为 GeoDataFrame
geoms = list(results)
# geoms = list(results)会将所有字典转换为一个列表,最终传递给 GeoDataFrame.from_features() 方法
gdf = gpd.GeoDataFrame.from_features(geoms)
gdf.crs = src.crs
# 导出矢量文件
gdf.to_file(r'C:\Users\HP\Desktop\test011\output_boundary02.shp')</pre>
2.栅格分块
对于大规模栅格数据,尤其是在导入 GIS 软件(如 ArcGIS)进行可视化和分析时,我们常常遇到加载缓慢甚至卡顿的问题。为了提升数据处理效率,栅格分块成为了一个重要的操作。通过将大栅格数据分成多个小块,我们可以将这些小块分别导入和处理,从而避免内存溢出或处理速度过慢的问题。
在实际操作中,我们可以根据大栅格的大小来确定分割的块数。例如,如果原始栅格文件大小为 30GB,而 ArcGIS 在导入 1GB 栅格数据时比较流畅,那么我们可以将栅格数据分割成 25 个小块。通过计算大栅格的长宽,除以横纵向的块数,即可得出每个小块栅格的尺寸。然后,通过获取每个小块的起始像素位置,并更新相应的仿射变换信息,我们就能够精确地进行栅格数据分块。

import rasterio
from rasterio.windows import Window
from rasterio.transform import from_origin
import os
input_ras=r"C:\Users\HP\Desktop\image_process\projected\ldst_rgb.tif"
output_folder=r"E:\image_download\yuxian2013"
#准备将大栅格纵向分成几块,横向分成几块
x=6
y=6
def split_raster_rio(input_raster, output_dir, rows=x, cols=y):
with rasterio.open(input_raster) as src:
# 得到影像的宽度和宽度,与横向、纵向分块数相处,就能得到每块栅格的尺寸
width = src.width
height = src.height
transform = src.transform
profile = src.profile
# 每块的宽度和高度
tile_width = width // cols
tile_height = height // rows
for row in range(rows):
for col in range(cols):
# 计算窗口起始像素位置
x_off = col * tile_width
y_off = row * tile_height
# 处理边界(防止最后一块越界),这里也就是最后一块就是总宽度减去最后一块起点位置
win_width = tile_width if col < cols - 1 else width - x_off
win_height = tile_height if row < rows - 1 else height - y_off
window = Window(x_off, y_off, win_width, win_height)#创建窗口
#这一行代码将根据Window对象来计算栅格图像的仿射变换,用来更新小块的transform
transform_tile = rasterio.windows.transform(window, transform)
tile_data = src.read(window=window)
# 更新原数据 profile
tile_profile = profile.copy()
tile_profile.update({
'height': win_height,
'width': win_width,
'transform': transform_tile})
# 输出路径
filename = os.path.splitext(os.path.basename(input_raster))[0]
out_path = os.path.join(output_dir, f"{filename}_tile_{row}_{col}.tif")
with rasterio.open(out_path, 'w', **tile_profile) as dst:
dst.write(tile_data)
split_raster_rio(input_ras,output_folder)