geo-rasterize：用于地理空间应用的 2D 光栅化器，用 Rust 编写

geo-rasterize是一个用于rust 库的 Python 包装器，它具有相同的名称，可以栅格化 匀称的 矢量形状，就像 rasterio的 features.rasterize一样。不同之处在于，虽然rasterio 是在 GDAL 上构建的，但这个库没有依赖关系——您可以安装它而不必担心 GDAL（或任何其他 C 库）！Plusgeo-rasterize的光栅化算法是基于 GDAL 的，所以它应该是一个替代品rasterio.features.rasterize ，它提供了一个非常相似的 API。

我们为以下平台发布了轮子到 PyPI Python 3.7+：

例子

例如，让我们将Point位于的单个栅格化(x=1, y=2) 到一个 5 像素宽和 6 像素高的栅格上。默认情况下，光栅像素将从零开始，我们将1在点接触的每个像素中放置一个：

>>> from shapely.geometry import Point
>>> from geo_rasterize import rasterize
>>> print(rasterize([Point(1, 2)], [1], (5, 6)))
[[0 0 0 0 0]
 [0 0 0 0 0]
 [0 1 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]]

结果正是我们所期望的：一个 5x6numpy数组，其中一个像素设置为 1！请注意，我们必须指定一个形状列表，而不仅仅是一个形状。并且前景值列表（[1]在这种情况下）必须具有相同的长度，因为我们需要为每个形状提供一个前景值。

所以让我们看看多种形状！

>>> from shapely.geometry import Point, LineString
>>> from geo_rasterize import rasterize
>>> shapes = [Point(3, 4), LineString([(0, 3), (3, 0)])]
>>> foregrounds = [3, 7]
>>> raster_size = (4, 5)
>>> print(rasterize(shapes, foregrounds, raster_size))
[[0 0 7 0]
 [0 7 7 0]
 [7 7 0 0]
 [7 0 0 0]
 [0 0 0 3]]

当两个形状在同一个像素中燃烧时会发生什么？这取决于您如何设置由algorithm 参数给出的合并算法。默认值'replace'意味着最后一个形状会覆盖像素，但您也可以将其设置'add'为前景值将相加。这使您可以制作热图！

>>> from shapely.geometry import Point, LineString
>>> from geo_rasterize import rasterize
>>> shapes = [LineString([(0, 0), (5, 5)]), LineString([(5, 0), (0, 5)])]
>>> print(rasterize(shapes, [1, 1], (5, 5), algorithm='add'))
[[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]]

我们的两条线在您会发现的中心交叉2。

您可以使用参数更改默认值background。dtype您可以使用参数设置数组dtype。设置 dtype='uint8'将使光栅阵列所需的空间减少 8 倍。如果您只对二进制光栅化感兴趣，这将特别有用。

地理变换

到目前为止，我们所有的例子都假设我们的形状坐标在图像空间中。换句话说，我们假设x 坐标将在范围内0..width，并且y坐标将在范围内0..height。唉，通常情况并非如此！

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

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

geo-raterize可以通过处理仿射变换来缓解这个繁琐的过程。geo_to_pix只需使用参数传递一个仿射变换数组 （.to_gdal()如果有 affine.Affine实例则调用）。

总而言之，您必须：

• 从您的图像中提取 CRS 并将您的形状转换为该 CRS（可能使用pyproj 及其与 [geo types][geo] 的集成，
• pix_to_geo从图像元数据中提取转换
• 从该数据创建一个Affine实例（GDAL 将这些表示为一个[f64; 6]数组，您可以使用Affine.from_gdal）
• 调用transform.inverse以获取相应的geo_to_pix 变换（请记住，并非所有变换都是可逆的！）
• 调用transform.to_gdal()并使用带有 geo_to_pix参数的结果数组