编写在Python中如何在DSM影像上划分格网的代码

以下是在Python中在DSM影像上划分格网的示例代码：

import gdal
import numpy as np

def create_grid(raster_file, grid_size):
    """
    Create a grid on top of a DSM raster file.
    :param raster_file: Path to DSM raster file.
    :param grid_size: Size of grid cells in meters.
    :return: A numpy array representing the grid, with each cell value corresponding to a grid cell ID.
    """
    # Open raster file
    ds = gdal.Open(raster_file)
    band = ds.GetRasterBand(1)
    nodata_val = band.GetNoDataValue()
    geo_transform = ds.GetGeoTransform()
    x_origin = geo_transform[0]
    y_origin = geo_transform[3]
    cell_width = geo_transform[1]
    cell_height = geo_transform[5]
    cols = ds.RasterXSize
    rows = ds.RasterYSize
    data = band.ReadAsArray(0, 0, cols, rows)
    
    # Create grid
    grid_cols = int(np.ceil(float(cols) * cell_width / grid_size))
    grid_rows = int(np.ceil(float(rows) * abs(cell_height) / grid_size))
    grid = np.zeros((grid_rows, grid_cols), dtype=np.int32)

    for i in range(grid_rows):
        for j in range(grid_cols):
            x = x_origin + j * grid_size
            y = y_origin + i * grid_size
            col = int((x - geo_transform[0]) / geo_transform[1])
            row = int((y - geo_transform[3]) / geo_transform[5])
            
            # Check if cell is valid and assign grid value
            if col >= 0 and col < cols and row >= 0 and row < rows:
                if data[row, col] != nodata_val:
                    grid[i, j] = int(i * grid_cols + j + 1)
    return grid

请注意，这个代码仅仅是一个示例，实际应用可能需要根据特定的数据集和需求进行调整和修改。此外，对于较大的数据集，该代码可能需要进行优化以提高执行效率。