Python地理數據處理之使用GR進行矢量的示例分析

這篇文章主要介紹Python地理數據處理之使用GR進行矢量的示例分析，文中介紹的非常詳細，具有一定的參考價值，感興趣的小伙伴們一定要看完！

1、疊加分析

??疊加分析操作：

??plot顏色：‘r’ 紅色， ‘g’ 綠色， ‘b’ 藍色， ‘c’ 青色， ‘y’ 黃色， ‘m’ 品紅， ‘k’ 黑色， ‘w’ 白色。

??新奧爾良城市邊界、水體和濕地的簡單地圖：

??1.新奧爾良城市沼澤區(qū)域分析：

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.計算城市的濕地面積：

# 獲得城市內的濕地多邊形# 將多邊形的面積進行累加# 除以城市面積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

??注：通過空間過濾和屬性過濾，將不必要的要素過濾，這樣可以顯著減少處理時間。

??3.兩圖層求交：

# 將湖泊數據排除# 在內存中創(chuàng)建一個臨時圖層# 將圖層相交，將結果儲存在臨時圖層中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、鄰近分析（確定要素間的距離）

??OGR包含兩個鄰近分析工具：測量幾何要素的距離；創(chuàng)建緩沖區(qū)。

??1.確定美國有多少城市位于火山10英里（1英里=1609.3米）的范圍之內。確定火山附近城市數量的存在問題的方法：

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')# 在內存中創(chuàng)建一個臨時層來存儲緩沖區(qū)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)# 將城市圖層與火山緩沖區(qū)圖層相交result_lyr = memory_ds.CreateLayer('result')buff_lyr.Intersection(cities_lyr, result_lyr)print('Cities: {}'.format(result_lyr.GetFeatureCount()))

Cities: 83

??2.一個更好地確定火山附近城市數量方法：

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')# 將緩沖區(qū)添加到一個復合多邊形，而不是一個臨時圖層multipoly = ogr.Geometry(ogr.wkbMultiPolygon)for volcano_feat in volcano_lyr:
    buff_geom = volcano_feat.geometry().Buffer(16000)
    multipoly.AddGeometry(buff_geom)# 將所有的緩沖區(qū)聯(lián)合在一起得到一個可以使用的多邊形作為空間過濾器cities_lyr.SetSpatialFilter(multipoly.UnionCascaded())print('Cities: {}'.format(cities_lyr.GetFeatureCount()))

Cities: 78

注：UnionCascaded()：有效地將所有的多邊形合并成一個復合多邊形
??第一個例子中，每當城市位于火山緩沖區(qū)內，就會復制到輸出結果中。說明一個城市位于多個16000米緩沖區(qū)內，將被列入不止一次。

??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))

92656 meters (57.91 miles)

??3. 用2.5D幾何對象，表示兩點之間的距離：

# 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。

# 用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。

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

True
True

以上是“Python地理數據處理之使用GR進行矢量的示例分析”這篇文章的所有內容，感謝各位的閱讀！希望分享的內容對大家有幫助，更多相關知識，歡迎關注億速云行業(yè)資訊頻道！