geo-rasterize: 一个针对地理空间应用的纯Rust 2D光栅化器

此crate旨在帮助那些有一些矢量数据（如 geo::Polygon）和一个光栅源（如可能使用 GDAL 打开的 GeoTiff）的人，他们希望生成一个表示多边形填充哪些光栅位的布尔数组。还有一个可用的 Python包装器。

此实现基于 GDAL 的 GDALRasterizeGeometries，允许您光栅化 geo-types 包支持的任何类型，包括

• 点
• 线
• 线字符串
• 多边形
• 矩形 和 三角形
• 多点、多线字符串 和 多边形
• 几何形状 和 几何形状集合

这些形状可以具有任何坐标和任何数值类型，只要它可以转换为 f64。

此crate在GDAL提供 ALL_TOUCHED=TRUE 选项时与GDAL的行为相匹配。因此，如果您只需要一个与GDAL兼容的光栅化器，则可以使用它作为GDAL的直接替换。此外，不支持GDAL的 BURN_VALUE_FROM=Z。但除此之外，此代码应该与GDAL光栅化器产生相同的结果——光栅化算法是直接移植的。我们使用 proptest 对GDAL执行随机差异比较，以提高我们对兼容性的信心。

动机：卫星影像数据分析

假设您对ESA的Sentinel-2卫星任务提供的免费10米分辨率数据感兴趣。[Sentinel-2]。您可能特别感兴趣的是农田随时间的变化，您有一批以多边形表示的农田。您从[AWS][https://registry.opendata.aws/sentinel-2-l2a-cogs/]获取了一些Sentinel-2数据。由于Sentinel-2瓦片非常大（超过1.1亿像素！）并且以[Cloud Optimized GeoTiffs][https://www.cogeo.org/]存储，因此您将多边形转换为窗口，并从图像瓦片中提取这些窗口；这样，您只需下载您关心的像素。

但现在您遇到了问题！您的多边形并不是与Sentinel-2瓦片系统对齐的完美矩形。因此，虽然您有小的农田芯片，但您不知道这些芯片的哪些部分对应于您的多边形。更糟糕的是，您的某些多边形有孔（例如，表示农田上的房屋或池塘）。这正是[geo-rasterize]发挥作用的地方！使用[geo-rasterize]，您可以将您的农田多边形转换为与Sentinel-2农田芯片类似的二值栅格。然后，您可以使用这些掩码芯片来选择您关心的Sentinel-2芯片的像素。通过在掩码上过滤，您现在可以生成影像的时间序列，确信您只检查您多边形内的像素！

二值栅格化

假设您要将多边形栅格化到4像素宽、5像素高的网格。为此，您只需使用[BinaryBuilder]构造一个[BinaryRasterizer]，使用[rasterize]调用您的多边形，然后使用[finish]调用以获得布尔值[Array2]。

# fn main() -> geo_rasterize::Result<()> {
use geo::polygon;
use ndarray::array;
use geo_rasterize::BinaryBuilder;

let poly = polygon![
    (x:4, y:2),
    (x:2, y:0),
    (x:0, y:2),
    (x:2, y:4),
    (x:4, y:2),
];

let mut r = BinaryBuilder::new().width(4).height(5).build()?;
r.rasterize(&poly)?;
let pixels = r.finish();

assert_eq!(
    pixels.mapv(|v| v as u8),
    array![
        [0, 1, 1, 0],
        [1, 1, 1, 1],
        [1, 1, 1, 1],
        [0, 1, 1, 1],
        [0, 0, 1, 0]
    ]
);
# Ok(()) }

...带有多个形状

但您想栅格化多个几何形状怎么办？这很简单！

# fn main() -> geo_rasterize::Result<()> {
use geo::{Geometry, Line, Point};
use ndarray::array;
use geo_rasterize::BinaryBuilder;

let shapes: Vec<Geometry<i32>> =
    vec![Point::new(3, 4).into(),
	     Line::new((0, 3), (3, 0)).into()];

let mut r = BinaryBuilder::new().width(4).height(5).build()?;
for shape in shapes {
    r.rasterize(&shape)?;
}

let pixels = r.finish();
assert_eq!(
    pixels.mapv(|v| v as u8),
    array![
        [0, 0, 1, 0],
        [0, 1, 1, 0],
        [1, 1, 0, 0],
        [1, 0, 0, 0],
        [0, 0, 0, 1]
    ]
);
# Ok(())}

标签（非二值栅格化）

到目前为止，我们一直在生成二进制数组；如果您想将不同的形状栅格化到相同的整数数组，存储每个像素对应于每个形状的不同值怎么办？为此，我们有一个[Rasterizer]，我们使用[LabelBuilder]来构建。当您使用[Rasterizer]烧录形状时，您不仅提供形状，还提供前景标签。但在您烧录任何内容之前，您必须指定一个背景标签，用于填充空的栅格数组。

# use geo_rasterize::{Result, LabelBuilder, Rasterizer};
# fn main() -> Result<()> {
use geo::{Geometry, Line, Point};
use ndarray::array;

let point = Point::new(3, 4);
let line = Line::new((0, 3), (3, 0));

let mut rasterizer = LabelBuilder::background(0).width(4).height(5).build()?;
rasterizer.rasterize(&point, 3)?;
rasterizer.rasterize(&line, 7)?;

let pixels = rasterizer.finish();
assert_eq!(
    pixels.mapv(|v| v as u8),
    array![
        [0, 0, 7, 0],
        [0, 7, 7, 0],
        [7, 7, 0, 0],
        [7, 0, 0, 0],
        [0, 0, 0, 3]
    ]
);
# Ok(())}

热图

如果两个形状接触相同的像素会发生什么？在上面的例子中，最后写入的形状获胜。但您可以通过指定不同的值来更改此行为，使用LabelBuilder::algorithm为[MergeAlgorithm]。实际上，使用MergeAlgorithm::Add，您可以轻松创建一个热图，显示形状密度，每个像素值都告诉您有多少形状落在该像素上！

# use geo_rasterize::{Result, LabelBuilder, Rasterizer, MergeAlgorithm};
# use geo::{Geometry, Line, Point};
# use ndarray::array;
# fn main() -> Result<()> {

let lines = vec![Line::new((0, 0), (5, 5)), Line::new((5, 0), (0, 5))];

let mut rasterizer = LabelBuilder::background(0)
    .width(5)
    .height(5)
    .algorithm(MergeAlgorithm::Add)
    .build()?;
for line in lines {
    rasterizer.rasterize(&line, 1)?;
}

let pixels = rasterizer.finish();
assert_eq!(
    pixels.mapv(|v| v as u8),
    array![
        [1, 0, 0, 0, 1],
        [0, 1, 0, 1, 1],
        [0, 0, 2, 1, 0],
        [0, 1, 1, 1, 0],
        [1, 1, 0, 0, 1]
    ]
);

# Ok(())}

两条线在中心相交，在那里您会发现2。请注意，[Rasterizer]不仅限于整数；任何可复制的类型都可以。[Rasterizer]提供了类似于[rasterio.readthedocs.io/][https://rasterio.readthedocs.io/en/latest/api/rasterio.features.html#rasterio.features.rasterize]的函数。

地理变换

到目前为止，我们所有的示例都假设我们的形状坐标是在图像空间中。换句话说，我们假设x坐标将在范围0到宽度之间，y坐标将在范围0到高度之间。但遗憾的是，情况往往并非如此！

对于卫星影像（或广义上的遥感影像），在元数据中几乎总是指定了坐标参考系统（CRS）和仿射变换。有关详细信息，请参阅rasterio的地理配准。

为了处理大多数影像，您需要将您的矢量形状从其原始的CRS（通常是地理经纬度的EPSG:4326）转换为数据文件指定的CRS（通常是UTM投影，但有很多选择）。然后，您需要应用仿射变换将世界坐标转换为像素坐标。由于栅格影像通常指定逆变换矩阵（即pix_to_geo变换），您首先需要将其求逆以获得geo_to_pix变换，然后才能将其应用于坐标。现在您已经获得了适合您图像数据的像素坐标！

[BinaryRasterizer]和[Rasterizer]可以通过处理仿射变换来简化这一繁琐的过程。请确保将[Transform]对象传递给[BinaryBuilder]或[LabelBuilder]。在两种情况下，该变换都是一个geo_to_pix变换，这意味着您需要

• 从您的图像中提取CRS并将您的形状转换为该CRS（可能使用proj crate及其与[geo类型][geo]的集成），
• 从您的影像元数据中提取pix_to_geo变换
• 从这些数据创建一个[Transform]实例（GDAL将这些表示为一个[f64; 6]数组）
• 调用transform.inverse以获取相应的geo_to_pix变换（由于并非所有变换都是可逆的，inverse返回一个Option）
• 将生成的[Transform]传递给[BinaryBuilder]或[LabelBuilder]。