Python geographic data processing analysis uses GR for vectors-Python Tutorial-php.cn

Python geographic data processing analysis uses GR for vectors

Free learning recommendations: python video tutorial

1. Overlay analysis

Overlay analysis operation:
Python geographic data processing analysis uses GR for vectors
Plot color: 'r' red, 'g' green, 'b' Blue, 'c' cyan, 'y' yellow, 'm' magenta, 'k' black, 'w' white.

# syntax

Simple map of New Orleans city boundaries, water bodies and wetlands: ## sync

synopsis Python geographic data processing analysis uses GR for vectors 1. Analysis of urban swamp areas in New Orleans:

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

Copy after login

Python geographic data processing analysis uses GR for vectors 2. Calculate the wetland area of the city: Filtering out unnecessary features can significantly reduce processing time.
　3. Intersection of two layers:

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

Copy after login

28.7% of New Orleans is wetland

Copy after login

2. Proximity analysis (determining the distance between elements)

OGR includes two proximity analysis tools: measuring distances between geometric features and creating buffers.

1. Determine how many cities in the United States are within 10 miles of a volcano.

A problematic way to determine the number of cities near a volcano:

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

Copy after login

28.7% of New Orleans is wetland

Copy after login

2. A better way to determine the number of cities near a volcano:

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

Copy after login

Cities: 83

Copy after login

Note: UnionCascaded(): Effectively combines all polygons into one composite polygon

In the first example, whenever a city is within the volcano buffer, it will be copied to the output result. Note that a city located within multiple 16,000-meter buffer zones will be included more than once. # synchronization

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

Copy after login

Cities: 78

Copy after login

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

Copy after login

92656 meters (57.91 miles)

Copy after login

　 Taking the elevation Z value into account, the true distance is 5.

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

Copy after login

4.0
Copy after login
Copy after login

# 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))
Copy after login

4.0
Copy after login
Copy after login
The area of 2.5D is actually 141.

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

Copy after login

rrree

Related free learning recommendations:

python tutorial

(Video)

The above is the detailed content of Python geographic data processing analysis uses GR for vectors. For more information, please follow other related articles on the PHP Chinese website!