python shp文件_使用python生成shapefile

在卫星遥感图像处理中，经常要使用shpfile文件(文件后缀为.shp)作为label或者辅助数据。同时，在后期分析中，也会利用shapefile，在QGIS中分析tif图像的特征和属性。所以生成符合自己需求的shapefile十分重要。这篇文章将会介绍如何利用python+gdal+QGIS生成shapefile,并且可视化。

原始数据Tif 图像 (10m/pixel)与图像坐标相对应的信息(本文中使用一个大的shapefile，里面提供了荷兰不同区域的各种统计调查数据，包括人口，房屋，车辆等)

通过tif图像可以得到该图像的坐标信息，我们的目标就是生成一个与这个图像坐标相对应的shapefile，其中包含该区域的城市名，人口密度等(这些数据可以根据需求自己设置)。

获取坐标信息

可以利用rasterio来读取tif文件，获取该图像的坐标信息。

import rasterio

img = rasterio.open('your_path_to_tiff/Nerland_60.tif')

print(img.bounds)

print(img.crs)

同时在QGIS中打开tif图像和shapefile，可以发现该区域的相关信息：区域名称：Gooise Meern

人均汽车保有量：1.2

区域总人口：17410

区域房屋总量：7675

创建shapefile

import osgeo.ogr as ogr

import osgeo.osr as osr

driver = ogr.GetDriverByName("ESRI Shapefile")

data_source = driver.CreateDataSource("Gooise.shp")

srs = osr.SpatialReference()

srs.ImportFromEPSG(32631)

layer = data_source.CreateLayer("Gooise", srs, ogr.wkbMultiPolygon)

field_name = ogr.FieldDefn("Name", ogr.OFTString)

field_name.SetWidth(24)

layer.CreateField(field_name)

field_name = ogr.FieldDefn("Car_per_house", ogr.OFTString)

field_name.SetWidth(24)

layer.CreateField(field_name)

feature = ogr.Feature(layer.GetLayerDefn())

feature.SetField("Name", 'Gooise')

feature = ogr.Feature(layer.GetLayerDefn())

feature.SetField("Car_per_house", '1.2')

wkt = 'polygon((646080 5797460,648640 5797460,648640 5794900,646080 5794900,646080 5797460))'

polygon = ogr.CreateGeometryFromWkt(wkt)

feature.SetGeometry(polygon)

layer.CreateFeature(feature)

feature = None

data_source = None

代码解释：

首先要创建一个操作shapefile的driver, 然后设置好对应的crs。shapefile的形状可以自己设置，可以是polygon, 可以是point,可以是linestring, 具体可以看维基上对应的介绍。需要格外注意的一点是，上面步骤完成以后一定不要忘了最后两行释放资源，否则不会成功生成shapefile。

结果生成的shapefile与原tif图像完全重合，右上角可以看到新添加的属性。