python实现使用GDAL实现矢量转栅格
2019-11-28 本文已影响0人
GIS无情老博士
需求
现在有一个shp文件和栅格数据,需要将shp转换成和栅格数据空间位置一致且像元大小一致的栅格数据。ArcGIS速度比较慢,使用GDAL将shp文件转为和目标栅格同样大小的栅格
![](https://img.haomeiwen.com/i19853754/ab9f7b15824ba7f7.png)
![](https://img.haomeiwen.com/i19853754/80f8728480834a0d.png)
代码
from osgeo import gdal, gdalconst
from osgeo import ogr
rasterFile = 'F:/**0416.dat' # 原影像
shpFile = 'F:/**小麦.shp' # 裁剪矩形
dataset = gdal.Open(rasterFile, gdalconst.GA_ReadOnly)
geo_transform = dataset.GetGeoTransform()
cols = dataset.RasterXSize # 列数
rows = dataset.RasterYSize # 行数
x_min = geo_transform[0]
y_min = geo_transform[3]
pixel_width = geo_transform[1]
shp = ogr.Open(shpFile, 0)
m_layer = shp.GetLayerByIndex(0)
target_ds = gdal.GetDriverByName('GTiff').Create("F:/小麦test.dat", xsize=cols, ysize=rows, bands=1,
eType=gdal.GDT_Byte)
target_ds.SetGeoTransform(geo_transform)
target_ds.SetProjection(dataset.GetProjection())
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(0)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], m_layer, options=["ATTRIBUTE=val"]) # 跟shp字段给栅格像元赋值
# gdal.RasterizeLayer(target_ds, [1], m_layer) # 多边形内像元值的全是255
del dataset
del target_ds
shp.Release()
结果
最终栅格化的结果如下图所示
![](https://img.haomeiwen.com/i19853754/d65c0db0c90e2c6c.png)