在做遥感图像处理时,矢量数据与栅格数据的转换是常有的事,本节l将利用GDAL、ogr来完成矢量与栅格之间的转换.
1.矢量转栅格
import gdal import ogr import gdalconst #from osgeo import ogr, osr, gdal,gdalconst def shpRaster(shp,templatePic,output,field,nodata): ndsm = templatePic data = gdal.Open(ndsm, gdalconst.GA_ReadOnly) geo_transform = data.GetGeoTransform() proj=data.GetProjection() x_res = data.RasterXSize y_res = data.RasterYSize mb_v = ogr.Open(shp) mb_l = mb_v.GetLayer() target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Int16) target_ds.SetGeoTransform(geo_transform) target_ds.SetProjection(proj) band = target_ds.GetRasterBand(1) NoData_value = nodata band.SetNoDataValue(NoData_value) band.FlushCache() gdal.RasterizeLayer(target_ds, [1], mb_l, options=["ATTRIBUTE=%s"%field,'ALL_TOUCHED=TRUE']) target_ds = None if __name__ == '__main__': shp=r"D:\test\1.shp" #矢量数据路径 templatePic=r"D:\test\1.tif" #模板栅格,一个tif,地理变换信息从这里读,栅格大小与该栅格一致 output=r"D:\test" #保存路径 field="Value" #转换字段 nodata=0 #矢量空白区转换后的值 shpRaster(shp, templatePic, output, field, nodata)
2.栅格转矢量
from osgeo import ogr, osr, gdal def raster_to_shape(rasterfile, shapefile,name): data = gdal.Open(rasterfile, gdal.GA_ReadOnly) inband = data.GetRasterBand(1) drv = ogr.GetDriverByName('ESRI Shapefile') Polygon = drv.CreateDataSource(shapefile) prj = osr.SpatialReference() prj.ImportFromWkt(data.GetProjection()) ## 使用栅格的投影信息 Polygon_layer = Polygon.CreateLayer(name, srs=prj, geom_type=ogr.wkbPolygon) newField = ogr.FieldDefn('Value', ogr.OFTInteger) Polygon_layer.CreateField(newField) gdal.FPolygonize(inband, None, Polygon_layer, 0) if __name__ == "__main__": rasterfile = r"F:\algorithm\1.tif" shapefile= r"F:\algorithm" name="test" raster_to_shape(rasterfile, shapefile,name)