Python中使用GDAL和OGR库来实现将平面划分为小的平面元的代码示例

以下是示例代码：

from osgeo import ogr, osr

def grid_extent(geotransform, x_size, y_size):
    ulx, xres, xskew, uly, yskew, yres = geotransform
    lrx = ulx + (x_size * xres)
    lry = uly + (y_size * yres)
    return ulx, lrx, lry, uly

def create_grid(shp_path, raster_path, grid_size):
    # Open raster file
    raster = ogr.Open(raster_path)
    layer = raster.GetLayer()

    # Compute extent of raster
    geotransform = layer.GetExtent()
    x_min, x_max, y_min, y_max = grid_extent(geotransform, layer.RasterXSize, layer.RasterYSize)

    # Create new shapefile for grid features
    driver = ogr.GetDriverByName("ESRI Shapefile")
    grid_shp = driver.CreateDataSource(shp_path)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    layer = grid_shp.CreateLayer("grid", srs, ogr.wkbPolygon)

    # Create grid cells
    x_range = range(int(x_min), int(x_max), grid_size)
    y_range = range(int(y_min), int(y_max), grid_size)
    for x in x_range:
        for y in y_range:
            # Create new polygon feature
            feature = ogr.Feature(layer.GetLayerDefn())
            feature.SetField("id", f"{x}-{y}")

            # Add geometry to feature
            ring = ogr.Geometry(ogr.wkbLinearRing)
            ring.AddPoint(x, y)
            ring.AddPoint(x + grid_size, y)
            ring.AddPoint(x + grid_size, y + grid_size)
            ring.AddPoint(x, y + grid_size)
            ring.AddPoint(x, y)
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)
            feature.SetGeometry(poly)

            # Add feature to layer
            layer.CreateFeature(feature)

    # Cleanup
    raster = None
    grid_shp = None

这段代码使用GDAL和OGR库来创建一个矩形的网格。输入参数包括一个栅格文件的路径、一个输出的矢量文件的路径、以及每个网格的大小。输出的矢量文件包括了每个网格的四个角点，同时包含了其ID。