离散点转换点shapefile

1.导入库

使用有python自带的os, sys, csv，以及GDAL库。

GDAL的安装请参考另一篇博文：http://www.cnblogs.com/tianjg/p/8145629.html

import os
import sys
import csv

try:
    from osgeo import gdal, ogr, osr
except ImportError:
    import gdal, ogr, osr

2. 解析txt或csv文件

读取txt点要素的经纬度信息和对应的数字，组织形式如下：

108,34,349.0
108,35,55.8
....

 

3.利用python+gdal将包含离散点的文本转成点矢量

代码如下：

       # 新建输出的点shapefile
 1     shp_driver = ogr.GetDriverByName("ESRI Shapefile")
 2     if os.path.exists(out_shp):
 3         shp_driver.DeleteDataSource(out_shp)
 4     out_dataset = shp_driver.CreateDataSource(out_shp)
 5 　　 # 定义输出矢量的空间参考和属性
 6     out_srs = osr.SpatialReference()
 7     out_srs.ImportFromEPSG(4326)
 8     layer_name = os.path.splitext(os.path.basename(out_shp))[0]
 9     out_layer = out_dataset.CreateLayer(layer_name, 
10     geom_type=ogr.wkbPoint, srs=out_srs)
11 
12     # 创建shaple的表结构
13     field_type = ogr.OFTReal
14     lon_field = ogr.FieldDefn('lon', field_type)
15     lat_field = ogr.FieldDefn('lat', field_type)
16     value = ogr.FieldDefn('aqi', field_type)
17     out_layer.CreateField(lon_field)
18     out_layer.CreateField(lat_field)
19     out_layer.CreateField(value)
20 　　
21     for data in txt_data:
22 
23         data = str(data).split(',')
24 
25         # 创建点Geometry
26         point = ogr.Geometry(ogr.wkbPoint)
27         point.AddPoint(float(data[lon_ind]), float(data[lat_ind]))
28 
29         # 创建要素和数字
30         featureDefn = out_layer.GetLayerDefn()
31         out_feature = ogr.Feature(featureDefn)
32         out_feature.SetGeometry(point)
33         out_feature.SetField('lon', float(data[lon_ind]))
34         out_feature.SetField('lat', float(data[lat_ind]))
35 
36         out_feature.SetField('aqi', float(data[aqi_ind]))
37         out_layer.CreateFeature(out_feature)
38 
39         out_feature = None
40 　　# 关闭
41     out_layer = None
42     out_dataset = None

 

posted @ 2018-01-17 09:29  tianjg  阅读(751)  评论(0)    收藏  举报