遥感

python实现使用GDAL实现矢量转栅格

2019-11-28  本文已影响0人  GIS无情老博士

需求

现在有一个shp文件和栅格数据,需要将shp转换成和栅格数据空间位置一致且像元大小一致的栅格数据。ArcGIS速度比较慢,使用GDAL将shp文件转为和目标栅格同样大小的栅格


转换前
目标大小栅格

代码

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()

结果

最终栅格化的结果如下图所示


后.png
上一篇下一篇

猜你喜欢

热点阅读