Python地理数据处理之分析使用GR进行矢量

wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

免费学习推荐：python视频教程wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

1、叠加分析wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

  叠加分析操作：

  plot颜色：‘r’ 红色， ‘g’ 绿色， ‘b’ 蓝色， ‘c’ 青色， ‘y’ 黄色， ‘m’ 品红， ‘k’ 黑色， ‘w’ 白色。wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

  新奥尔良城市边界、水体和湿地的简单地图：wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

  1.新奥尔良城市沼泽区域分析：wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

import osfrom osgeo import ogrfrom ospybook.vectorplotter import VectorPlotter data_dir = r'E:Google chromeDownloadgis with pythonosgeopy 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.计算城市的湿地面积：wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

# 获得城市内的湿地多边形# 将多边形的面积进行累加# 除以城市面积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))28.7% of New Orleans is wetland

  注：通过空间过滤和属性过滤，将不必要的要素过滤，这样可以显著减少处理时间。wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

  3.两图层求交：wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

# 将湖泊数据排除# 在内存中创建一个临时图层# 将图层相交，将结果储存在临时图层中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))28.7% of New Orleans is wetland

2、邻近分析（确定要素间的距离）wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

  OGR包含两个邻近分析工具：测量几何要素的距离；创建缓冲区。wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

  1.确定美国有多少城市位于火山10英里（1英里=1609.3米）的范围之内。确定火山附近城市数量的存在问题的方法：wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

from osgeo import ogr shp_ds = ogr.Open(r'E:Google chromeDownloadgis with pythonosgeopy dataUS')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())) Cities: 83

  2.一个更好地确定火山附近城市数量方法：wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

from osgeo import ogr shp_ds = ogr.Open(r'E:Google chromeDownloadgis with pythonosgeopy dataUS')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()))Cities: 78

注：UnionCascaded()：有效地将所有的多边形合并成一个复合多边形
  第一个例子中，每当城市位于火山缓冲区内，就会复制到输出结果中。说明一个城市位于多个16000米缓冲区内，将被列入不止一次。wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

  3.计算特定的城市与火山的距离：wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

import osfrom osgeo import ogrfrom ospybook.vectorplotter import VectorPlotter data_dir = r'E:Google chromeDownloadgis with pythonosgeopy 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))92656 meters (57.91 miles)

  3. 用2.5D几何对象，表示两点之间的距离：wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

# 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))4.0# 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))4.0

  将高程Z值考虑进去，真正的距离是5。wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

# 用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())100.0# 用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())100.0

  2.5D的面积实际上是141。wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

# 叠加操作同样忽略了高程值Zprint(poly_2d.Contains(pt1_2d))print(poly_25d.Contains(pt1_2d))True True

相关免费学习推荐：python教程（视频）wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台

以上就是Python地理数据处理之分析使用GR进行矢量的详细内容，更多请关注少儿编程网其它相关文章！wkz少儿编程网-Scratch_Python_教程_免费儿童编程学习平台