利用gdal将渔网矢量转换为指定分辨率的掩膜

最新推荐文章于 2023-03-14 13:44:24 发布

放眼望地球

最新推荐文章于 2023-03-14 13:44:24 发布

阅读量226

点赞数

文章标签： python Powered by 金山文档

本文链接：https://blog.csdn.net/Cyj_0703/article/details/128714289

版权

该代码段展示了如何从ESRIShapefile读取矢量数据，遍历其特征，然后利用GDAL创建新的栅格文件。它涉及到了空间参考、分辨率设定、栅格化过程以及压缩选项。最终，每个矢量要素被转化为一个单波段的8位TIFF文件。

摘要由CSDN通过智能技术生成

from osgeo import gdal
from osgeo import ogr
from osgeo.gdalconst import GDT_Byte, GDT_Float32, GDT_UInt16
import os
from tqdm import tqdm
import random

shp_path = 'G:/shp/China_fishnet/China_fishnet.shp'
save_path ='G:/shp/mask/'
os.chdir(r"G:/shp/") #目标路径

#创建驱动器
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.Open(shp_path, 0)
layers = ds.GetLayer() #获取图层文件
spatial = layers.GetSpatialRef() #获取矢量文件的空间信息
numFeatures = layers.GetFeatureCount() #统计图层中包含的要素
pixel_size=0.00001 #分辨率(经纬度)

for i in tqdm(range(numFeatures)):  
    feature = layers.GetNextFeature()  #遍历矢量文件中的每一个要素
    #创建一个新的图层，将每一个要素考入到新的图层中
    extfile='gd.shp'
    newds = driver.CreateDataSource(extfile)
    layernew = newds.CreateLayer('gd',spatial,ogr.wkbPolygon) #创建图层
    
    #必须要创建相同的字段，才能将要素拷入到该图层中
    fieldDefn = ogr.FieldDefn('Name', ogr.OFTString)
    fieldDefn.SetWidth(10)
    layernew.CreateField(fieldDefn)
                
    layernew.CreateFeature(feature) #拷入当前要素
    x_min, x_max, y_min, y_max = layernew.GetExtent()#获取目标图层中矢量范围
    
    x_res = int((round((x_max - x_min),1) / pixel_size)) #计算生成tif的行列数，round是为了取整
    y_res = int((round((y_max - y_min),1) / pixel_size)) 
    
    spatialRef = layernew.GetSpatialRef() #获取目标图层中矢量的空间信息

    #设置生成栅格的环境，命名等
    id = feature.GetField('Name')    
    save_name=save_path+id+'.tif'  
    select= "Name=" + "'"+ str(id)+ "'"
    
    #创建目标栅格 
    target_ds = gdal.GetDriverByName('GTiff').Create(save_name, x_res, y_res, bands=1,
                                                 eType=gdal.GDT_Byte,options=["TILED=YES", "COMPRESS=LZW"])
    target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
    target_ds.SetProjection(str(spatialRef))     
    gdal.RasterizeLayer(target_ds, [1],layernew,burn_values=[random.randint(0,255)]) #矢量栅格化 这一步必须要用layer裁剪，要素不行

    del layernew, target_ds