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

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

1、疊加分析

??疊加分析操作：

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

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

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

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.兩圖層求交：

# 將湖泊數據排除# 在內存中創建一個臨時圖層# 將圖層相交，將結果儲存在臨時圖層中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包含兩個鄰近分析工具：測量幾何要素的距離；創建緩沖區。

??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')# 在內存中創建一個臨時層來存儲緩沖區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.一個更好地確定火山附近城市數量方法：

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

Cities: 78

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

??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進行矢量的示例分析”這篇文章的所有內容，感謝各位的閱讀！希望分享的內容對大家有幫助，更多相關知識，歡迎關注億速云行業資訊頻道！