python 矢量数据转格栅 自己实现

矢量数据转格栅是将矢量数据（如点、线、面等）转换为栅格数据的过程。Python提供了很多用于地理空间数据处理的库，其中最常用的是GDAL和Fiona库。

首先，我们需要安装这两个库。可以使用pip命令进行安装：

pip install gdal
pip install fiona

接下来，我们可以按照以下步骤来实现矢量数据转格栅：

1. 导入所需的库：

import fiona
from osgeo import gdal
from osgeo import ogr
from osgeo import osr

2. 打开矢量数据源，并获取投影信息和几何类型：

source = fiona.open('input.shp', 'r')
layer = source[0]
proj = layer.crs
geometry_type = layer['geometry']

3. 创建栅格数据：

driver = gdal.GetDriverByName('GTiff')
output = driver.Create('output.tif', width, height, 1, gdal.GDT_Byte)
output.SetProjection(proj)
output.GetRasterBand(1).SetNoDataValue(0)
output.GetRasterBand(1).Fill(0)

4. 根据矢量数据的类型进行不同的处理：
- 对于点数据，可以根据点的位置在栅格上设置像素值。

for feat in source:
    geom = feat['geometry']
    pixel_x = int((geom['coordinates'][0] - origin_x) / pixel_size)
    pixel_y = int((origin_y - geom['coordinates'][1]) / pixel_size)
    output.GetRasterBand(1).WriteRaster(pixel_x, pixel_y, 1, 1, 1, 1)

- 对于线数据，可以根据线的像素宽度在栅格上设置像素值。

for feat in source:
    geom = feat['geometry']
    for line in geom['coordinates']:
        for i in range(len(line) - 1):
            start = line[i]
            end = line[i + 1]
            pixel_start_x = int((start[0] - origin_x) / pixel_size)
            pixel_start_y = int((origin_y - start[1]) / pixel_size)
            pixel_end_x = int((end[0] - origin_x) / pixel_size)
            pixel_end_y = int((origin_y - end[1]) / pixel_size)
            
            output.GetRasterBand(1).DrawLine(pixel_start_x, pixel_start_y, pixel_end_x, pixel_end_y)

- 对于面数据，可以根据面的边界进行栅格化。

for feat in source:
    geom = feat['geometry']
    for polygon in geom['coordinates']:
        ring = ogr.Geometry(ogr.wkbLinearRing)
        for point in polygon:
            ring.AddPoint(point[0], point[1])
        polygon = ogr.Geometry(ogr.wkbPolygon)
        polygon.AddGeometry(ring)
        
        for pixel_x in range(output.RasterXSize):
            for pixel_y in range(output.RasterYSize):
                raster_geom = ogr.Geometry(ogr.wkbPoint)
                raster_geom.AddPoint(origin_x + pixel_x * pixel_size, origin_y - pixel_y * pixel_size)
                
                if polygon.Contains(raster_geom):
                    output.GetRasterBand(1).WriteRaster(pixel_x, pixel_y, 1, 1, 1, 1)

5. 保存和清理资源：

output.FlushCache()
output = None
source.close()

以上就是用Python实现矢量数据转格栅的简单步骤。根据具体情况，可能还需要进行投影转换、像素大小计算等处理。通过使用GDAL和Fiona库，我们可以很方便地处理地理空间数据，并进行矢量数据到栅格数据的转换。