gdal 压缩tif,使用GDAL从tif文件创建shapefile

I am using library gdal to load a tiff file and create a 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

regards,

解决方案

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)