本文摘自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进行矢量的详细内容，更多文章请关注木庄网络博客 ！！
相关阅读 >>
Python 的for循环语句怎么写
零基础如何学习Python
Python 如何实现数字炸弹游戏
Python 中count函数的用法详解
对Python 2.7 pandas 中的read_excel详解
Python 调用c# com dll组件的过程详解
厉害了，30行Python 代码爬取英雄联盟全英雄皮肤
Python 列表list
ubuntu怎么安装Python
ipad上可以运行Python 么
更多相关阅读请进入《Python 》频道 >>
￥69.8元
python入门书籍，非常畅销，超高好评，python官方公认好书。
转载请注明出处：木庄网络博客 » Python地理数据处理之分析使用GR进行矢量
