使用GDAL从tif文件创建shapefile

问题描述

我正在使用gdal库加载tiff文件并创建shapefile.

I am using library gdal to load a tiff file and create a shapefile.

当我使用QGIS GUI加载shapefile时，标高上没有任何信息.

When I load my shapefile with QGIS GUI, There are no informations on the elevation.

我想在转换时保持海拔高度.

I would like to keep the elevation while the transformation.

import os
from osgeo import gdal,ogr,osr,gdalnumeric
import numpy as np

# this allows GDAL to throw Python Exceptions
gdal.UseExceptions()


print "reading tif file..."
try:
    ds = gdal.Open( "file.tif" )
except RuntimeError, e:
    print 'Unable to open file'
    print e
    sys.exit(1)

try:
    srcband = ds.GetRasterBand(1)
except RuntimeError, e:
    # for example, try GetRasterBand(10)
    print 'Band ( %i ) not found' % band_num
    print e
    sys.exit(1)


band = ds.GetRasterBand(1)
# elevation 2D numpy array
elevation = band.ReadAsArray()

# create shapefile datasource from geotiff file
#
print "creating shapefile..."
dst_layername = "Shape"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )
gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None 

致谢

推荐答案

以下是您的脚本，该脚本已更新为创建一个名为海拔"的字段，并将波段1中的栅格值提取到该字段中.

Here is your script which is updated to create a field called "elevation" and extract the raster value in band 1 into that field.

import os
from osgeo import gdal,ogr
import sys 

# mapping between gdal type and ogr field type
type_mapping = { gdal.GDT_Byte: ogr.OFTInteger,
                 gdal.GDT_UInt16: ogr.OFTInteger,   
                 gdal.GDT_Int16: ogr.OFTInteger,    
                 gdal.GDT_UInt32: ogr.OFTInteger,
                 gdal.GDT_Int32: ogr.OFTInteger,
                 gdal.GDT_Float32: ogr.OFTReal,
                 gdal.GDT_Float64: ogr.OFTReal,
                 gdal.GDT_CInt16: ogr.OFTInteger,
                 gdal.GDT_CInt32: ogr.OFTInteger,
                 gdal.GDT_CFloat32: ogr.OFTReal,
                 gdal.GDT_CFloat64: ogr.OFTReal}

# this allows GDAL to throw Python Exceptions
gdal.UseExceptions()


print "reading tif file..."
try:
    ds = gdal.Open(sys.argv[1])
except RuntimeError, e:
    print 'Unable to open file'
    print e
    sys.exit(1)

try:
    srcband = ds.GetRasterBand(1)
except RuntimeError, e:
    # for example, try GetRasterBand(10)
    print 'Band ( %i ) not found' % 1
    print e
    sys.exit(1)

# create shapefile datasource from geotiff file
print "creating shapefile..."
dst_layername = "Shape"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( "poly_out.shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )
raster_field = ogr.FieldDefn('elevation', type_mapping[srcband.DataType])
dst_layer.CreateField(raster_field)
gdal.Polygonize( srcband, None, dst_layer, 0, [], callback=None)

这篇关于使用GDAL从tif文件创建shapefile的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！