使用rasterio进行矢量转栅格

网上流行的多是用gdal进行矢量转栅格，而gdal用起来需要考虑地理坐标信息等，容易出错。通过geopandas和rasterio同样可以实现。

需要注意的是对于没有地理信息的矢量和图像，需要特别处理下更新下地理信息，即翻转下y坐标。

    # for no geotransform data
    trans = raster.transform
    if trans.a == 1 and trans.e == 1 and trans.b == 0 and trans.c == 0:
        # origin = 1.0, 0.0, 0.0, 0.0, 1.0, 0.0  c a b f d e =
        geotransform = (0, 1, 0, 0, 0, -1)  # c a b f d e
        geotransform = Affine.from_gdal(*geotransform)

完整的矢量转栅格代码如下：

import os.path
from os.path import join, exists
import numpy as np
import rasterio as rio
from rasterio import features, enums
import geopandas as gpd
from tqdm import tqdm
from affine import Affine


def vector2img(vectorFileName, templateTifFileName, outputFileName=None, field=None):
    # Read in vector
    vector = gpd.read_file(vectorFileName)

    # Get list of geometries for all features in vector file
    geom = [shapes for shapes in vector.geometry]

    # Open example raster
    raster = rio.open(templateTifFileName)
    geotransform = raster.transform

    # for no geotransform data
    trans = raster.transform
    if trans.a == 1 and trans.e == 1 and trans.b == 0 and trans.c == 0:
        # origin = 1.0, 0.0, 0.0, 0.0, 1.0, 0.0  c a b f d e =
        geotransform = (0, 1, 0, 0, 0, -1)  # c a b f d e
        geotransform = Affine.from_gdal(*geotransform)

    if len(geom) > 0:  # only rasterize non-empty vector
        if field is None:
            rasterized = features.rasterize(geom,
                                            out_shape=raster.shape,
                                            fill=0,
                                            out=None,
                                            transform=geotransform, # raster.transform
                                            all_touched=False,
                                            default_value=1,
                                            dtype=None)
        else:
            # create a numeric unique value for each row
            vector[field] = range(0, len(vector))

            # create tuples of geometry, value pairs, where value is the attribute value you want to burn
            geom_value = ((geom, value) for geom, value in zip(vector.geometry, vector[field]))

            # Rasterize vector using the shape and transform of the raster
            rasterized = features.rasterize(geom_value,
                                            out_shape=raster.shape,
                                            transform=geotransform,  # raster.transform
                                            all_touched=True,
                                            fill=0,  # background value
                                            merge_alg=enums.MergeAlg.replace,
                                            dtype=np.int16)
    else:
        rasterized = np.zeros([raster.height, raster.width]).astype(np.uint8)

    if outputFileName:
        with rio.open(
                outputFileName, "w",
                driver="GTiff",
                transform=geotransform,
                dtype=rio.uint8,
                count=1,
                width=raster.width,
                height=raster.height,
                compress='lzw') as dst:
            dst.write(rasterized, indexes=1)
    else:
        return rasterized