本文摘自php中文网，作者coldplay.xixi，侵删。
免费学习推荐： python视频教程
1、叠加分析
??叠加分析操作：plot颜色： ‘r’ 红色， ‘g’ 绿色， ‘b’ 蓝色， ‘c’ 青色， ‘y’ 黄色， ‘m’ 品红， ‘k’ 黑色， ‘w’ 白色。
??新奥尔良城市边界、水体和湿地的简单地图：
1.新奥尔良城市沼泽区域分析：
1
2
3
import osfrom osgeo import ogrfrom ospybook.vectorplotter import VectorPlotter
data_dir = r'E:\Google chrome\Download\gis with python\osgeopy data'# 得到新奥尔良附近的一个特定的沼泽特征vp = VectorPlotter(True)water_ds = ogr.Open(os.path.join(data_dir, 'US', 'wtrbdyp010.shp'))water_lyr = water_ds.GetLayer(0)water_lyr.SetAttributeFilter('WaterbdyID = 1011327')marsh_feat = water_lyr.GetNextFeature()marsh_geom = marsh_feat.geometry().Clone()vp.plot(marsh_geom, 'c')# 获得新奥尔良边城市边界nola_ds = ogr.Open(os.path.join(data_dir, 'Louisiana', 'NOLA.shp'))nola_lyr = nola_ds.GetLayer(0)nola_feat = nola_lyr.GetNextFeature()nola_geom = nola_feat.geometry().Clone()vp.plot(nola_geom, fill=False, ec='red', ls='dashed', lw=3)# 相交沼泽和边界多边形得到沼泽的部分# 位于新奥尔良城市边界内intersection = marsh_geom.Intersection(nola_geom)vp.plot(intersection, 'yellow', hatch='x')vp.draw()
2.计算城市的湿地面积：
1
2
3
# 获得城市内的湿地多边形# 将多边形的面积进行累加# 除以城市面积water_lyr.SetAttributeFilter("Feature != 'Lake'") # 限定对象water_lyr.SetSpatialFilter(nola_geom)wetlands_area = 0# 累加多边形面积for feat in water_lyr:
intersect = feat.geometry().Intersection(nola_geom)
wetlands_area += intersect.GetArea()pcnt = wetlands_area / nola_geom.GetArea()print('{:.1%} of New Orleans is wetland'.format(pcnt))
1
28.7% of New Orleans is wetland
??注 ：通过空间过滤和属性过滤，将不必要的要素过滤，这样可以显著减少处理时间。
??3.两图层求交：
1
2
3
# 将湖泊数据排除# 在内存中创建一个临时图层# 将图层相交，将结果储存在临时图层中water_lyr.SetAttributeFilter("Feature != 'Lake'")water_lyr.SetSpatialFilter(nola_geom)wetlands_area = 0for feat in water_lyr:
intersect = feat.geometry().Intersection(nola_geom) # 求交
wetlands_area += intersect.GetArea()pcnt = wetlands_area / nola_geom.GetArea()print('{:.1%} of New Orleans is wetland'.format(pcnt))water_lyr.SetSpatialFilter(None)water_lyr.SetAttributeFilter("Feature != 'Lake'")memory_driver = ogr.GetDriverByName('Memory')temp_ds = memory_driver.CreateDataSource('temp')temp_lyr = temp_ds.CreateLayer('temp')nola_lyr.Intersection(water_lyr, temp_lyr)sql = 'SELECT SUM(OGR_GEOM_AREA) AS area FROM temp'lyr = temp_ds.ExecuteSQL(sql)pcnt = lyr.GetFeature(0).GetField('area') / nola_geom.GetArea()print('{:.1%} of New Orleans is wetland'.format(pcnt))
1
28.7% of New Orleans is wetland
2、邻近分析（确定要素间的距离）
??OGR包含两个邻近分析工具：测量几何要素的距离；创建缓冲区。
??1.确定美国有多少城市位于火山10英里（1英里=1609.3米）的范围之内。确定火山附近城市数量的存在问题的方法：
1
2
3
4
5
6
from osgeo import ogr
shp_ds = ogr.Open(r'E:\Google chrome\Download\gis with python\osgeopy data\US')volcano_lyr = shp_ds.GetLayer('us_volcanos_albers')cities_lyr = shp_ds.GetLayer('cities_albers')# 在内存中创建一个临时层来存储缓冲区memory_driver = ogr.GetDriverByName('memory')memory_ds = memory_driver.CreateDataSource('temp')buff_lyr = memory_ds.CreateLayer('buffer')buff_feat = ogr.Feature(buff_lyr.GetLayerDefn())# 缓缓冲每一个火山点，将结果添加到缓冲图层中for volcano_feat in volcano_lyr:
buff_geom = volcano_feat.geometry().Buffer(16000)
tmp = buff_feat.SetGeometry(buff_geom)
tmp = buff_lyr.CreateFeature(buff_feat)# 将城市图层与火山缓冲区图层相交result_lyr = memory_ds.CreateLayer('result')buff_lyr.Intersection(cities_lyr, result_lyr)print('Cities: {}'.format(result_lyr.GetFeatureCount()))
??2.一个更好地确定火山附近城市数量方法：
1
2
3
4
5
from osgeo import ogr
shp_ds = ogr.Open(r'E:\Google chrome\Download\gis with python\osgeopy data\US')volcano_lyr = shp_ds.GetLayer('us_volcanos_albers')cities_lyr = shp_ds.GetLayer('cities_albers')# 将缓冲区添加到一个复合多边形，而不是一个临时图层multipoly = ogr.Geometry(ogr.wkbMultiPolygon)for volcano_feat in volcano_lyr:
buff_geom = volcano_feat.geometry().Buffer(16000)
multipoly.AddGeometry(buff_geom)# 将所有的缓冲区联合在一起得到一个可以使用的多边形作为空间过滤器cities_lyr.SetSpatialFilter(multipoly.UnionCascaded())print('Cities: {}'.format(cities_lyr.GetFeatureCount()))
注 ：UnionCascaded()：有效地将所有的多边形合并成一个复合多边形
??3.计算特定的城市与火山的距离：
1
2
3
import osfrom osgeo import ogrfrom ospybook.vectorplotter import VectorPlotter
data_dir = r'E:\Google chrome\Download\gis with python\osgeopy data'shp_ds = ogr.Open(os.path.join(data_dir, 'US'))volcano_lyr = shp_ds.GetLayer('us_volcanos_albers')cities_lyr = shp_ds.GetLayer('cities_albers')# 西雅图到雷尼尔山的距离volcano_lyr.SetAttributeFilter("NAME = 'Rainier'")feat = volcano_lyr.GetNextFeature()rainier = feat.geometry().Clone()cities_lyr.SetSpatialFilter(None)cities_lyr.SetAttributeFilter("NAME = 'Seattle'")feat = cities_lyr.GetNextFeature()seattle = feat.geometry().Clone()meters = round(rainier.Distance(seattle))miles = meters / 1600print('{} meters ({} miles)'.format(meters, miles))
1
92656 meters (57.91 miles)
??3. 用2.5D几何对象，表示两点之间的距离：
1
# 2Dpt1_2d = ogr.Geometry(ogr.wkbPoint)pt1_2d.AddPoint(15, 15)pt2_2d = ogr.Geometry(ogr.wkbPoint)pt2_2d.AddPoint(15, 19)print(pt1_2d.Distance(pt2_2d))
1
# 2.5Dpt1_25d = ogr.Geometry(ogr.wkbPoint25D)pt1_25d.AddPoint(15, 15, 0)pt2_25d = ogr.Geometry(ogr.wkbPoint25D)pt2_25d.AddPoint(15, 19, 3)print(pt1_25d.Distance(pt2_25d))
??将高程Z值考虑进去，真正的距离是5。
1
# 用2D计算面积ring = ogr.Geometry(ogr.wkbLinearRing)ring.AddPoint(10, 10)ring.AddPoint(10, 20)ring.AddPoint(20, 20)ring.AddPoint(20, 10)poly_2d = ogr.Geometry(ogr.wkbPolygon)poly_2d.AddGeometry(ring)poly_2d.CloseRings()print(poly_2d.GetArea())
1
# 用2.5D计算面积ring = ogr.Geometry(ogr.wkbLinearRing)ring.AddPoint(10, 10, 0)ring.AddPoint(10, 20, 0)ring.AddPoint(20, 20, 10)ring.AddPoint(20, 10, 10)poly_25d = ogr.Geometry(ogr.wkbPolygon25D)poly_25d.AddGeometry(ring)poly_25d.CloseRings()print(poly_25d.GetArea())
??2.5D的面积实际上是141。
1
# 叠加操作同样忽略了高程值Zprint(poly_2d.Contains(pt1_2d))print(poly_25d.Contains(pt1_2d))
相关免费学习推荐： python教程 （视频）
以上就是Python地理数据处理之分析使用GR进行矢量的详细内容，更多文章请关注木庄网络博客 ！！
相关阅读 >>
web自动化测试(三)selenium+beatuifulsoup
Python 数据分析师要学什么
Python 采集--数据的储存
Python 是哪个国家的人开发的语言
学习Python 的几个不同阶段
Python 实现在idle中输入多行的方法
Python 可以开发软件吗
Python 如何实现堆栈与队列的实例详解
Python 如何求列表平均值？
Python 中py2exe打包工具的用法详解
更多相关阅读请进入《Python 》频道 >>
￥69.8元
python入门书籍，非常畅销，超高好评，python官方公认好书。
转载请注明出处：木庄网络博客 » Python地理数据处理之分析使用GR进行矢量
