python实现栅格转点

栅格转点，即将栅格数据（如DEM高程数据）转换成点数据，可以使用Python中的GDAL库实现。下面是一个简单的示例代码：

from osgeo import gdal, ogr

# 打开栅格数据
raster_ds = gdal.Open('/path/to/raster.tif')

# 获取栅格数据的投影和地理变换信息
proj = raster_ds.GetProjection()
geotransform = raster_ds.GetGeoTransform()

# 获取栅格数据的行列数
cols = raster_ds.RasterXSize
rows = raster_ds.RasterYSize

# 创建点矢量文件
driver = ogr.GetDriverByName('ESRI Shapefile')
point_ds = driver.CreateDataSource('/path/to/point.shp')
layer = point_ds.CreateLayer('point', geom_type=ogr.wkbPoint)
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))

# 遍历栅格数据，将每个像元转换成点矢量并写入矢量文件
for i in range(rows):
    for j in range(cols):
        data = raster_ds.GetRasterBand(1).ReadAsArray(j, i, 1, 1)
        if data[0][0] is not None:
            x = geotransform[0] + j * geotransform[1] + i * geotransform[2]
            y = geotransform[3] + j * geotransform[4] + i * geotransform[5]
            point = ogr.Geometry(ogr.wkbPoint)
            point.AddPoint(x, y)
            feature = ogr.Feature(layer.GetLayerDefn())
            feature.SetGeometry(point)
            feature.SetField('id', i * cols + j)
            layer.CreateFeature(feature)
            feature = None

# 释放资源
raster_ds = None
point_ds = None

这个示例代码将从栅格数据（如DEM高程数据）中遍历每个像元，将其转换成点矢量，并写入一个新的点矢量文件中。需要注意的是，这里假设栅格数据的分辨率是正方形的，即像元的宽度和高度相等。如果不是正方形，需要根据实际情况进行调整。