Python地理数据处理 九：发电厂选址问题

  区域的合适性从1到7，一般数值大于等于3的区域是合适的。将这些区域与人口普查数据结合，选出风速合适并且每平方千米人口不到0.5%的地方。

  人口普查数据包含每个人口普查单元的人口数据，但没有人口密度数据。

  1. 要获得人口属性字段：

# 添加一个浮点字段，计算人口密度
# 通过HD01_S001字段获取人口普查数据
census_fn = os.path.join(data_dir, 'California', 'ca_census_albers.shp')
census_ds = ogr.Open(census_fn, True)
census_lyr = census_ds.GetLayer()
density_field = ogr.FieldDefn('popsqkm', ogr.OFTReal)
census_lyr.CreateField(density_field)
for row in census_lyr:
    pop = row.GetField('HD01_S001')
    sqkm = row.geometry().GetArea() / 1000000
    row.SetField('popsqkm', pop / sqkm)
    census_lyr.SetFeature(row)

# 获取因皮里尔县区域几何图像
county_fn = os.path.join(data_dir, 'US', 'countyp010.shp')
county_ds = ogr.Open(county_fn)
county_lyr = county_ds.GetLayer()
county_lyr.SetAttributeFilter("COUNTY ='Imperial County'")
county_row = county_lyr.GetNextFeature()
county_geom = county_row.geometry().Clone()

del county_ds

  2. 因为县数据使用的是经纬度坐标值，普查区和风能数据使用的是米，所以需要转换坐标系：

# 转换坐标系
county_geom.TransformTo(census_lyr.GetSpatialRef())
census_lyr.SetSpatialFilter(county_geom)

# 设置属性过滤器
# 人口密度低于0.5%
census_lyr.SetAttributeFilter('popsqkm < 0.5')

# 对风力数据集进行属性过滤，只显示等级为3或更好的区域：
wind_fn = os.path.join(data_dir, 'California', 'california_50m_wind_albers.shp')
wind_ds = ogr.Open(wind_fn)
wind_lyr = wind_ds.GetLayer()
wind_lyr.SetAttributeFilter('WPC >= 3')

# 创建一个新的shp文件，并将风力等级字段和人口密度字段加入：
out_fn = os.path.join(data_dir, 'California', 'wind_farm.shp')
out_ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource(out_fn)
out_lyr = out_ds.CreateLayer('wind_farm', wind_lyr.GetSpatialRef(), ogr.wkbPolygon)
out_lyr.CreateField(ogr.FieldDefn('wind', ogr.OFTInteger))
out_lyr.CreateField(ogr.FieldDefn('popsqkm', ogr.OFTReal))
out_row = ogr.Feature(out_lyr.GetLayerDefn())

创建结果：

  3. 将人口普查数据与风力数据相交，并储存到新创建的shp文件中：

#将人口普查数据与县界相交
for census_row in census_lyr:
    census_geom = census_row.geometry()
    census_geom = census_geom.Intersection(county_geom)
    wind_lyr.SetSpatialFilter(census_geom)

    print('Intersecting census tract with {0} wind polygons'.format(
        wind_lyr.GetFeatureCount()))

    # 检查是否存在风力多边形
    if wind_lyr.GetFeatureCount() > 0:
        out_row.SetField('popsqkm', census_row.GetField('popsqkm'))
        for wind_row in wind_lyr:
            wind_geom = wind_row.geometry()

            # 检查人口普查数据是否与风力多边形相交
            if census_geom.Intersect(wind_geom):
                new_geom = census_geom.Intersection(wind_geom)
                out_row.SetField('wind', wind_row.GetField('WPC'))
                out_row.SetGeometry(new_geom)
                out_lyr.CreateFeature(out_row)
                
del out_ds

运行结果：

  存在的问题： 虚线区域（黄色）代表县界数据，实线区域（蓝色）代表人口普查边界，两者并没有完全重合。所以，可以使用虚线以内的数据，进行风力发电厂选址。

  分析得到的适宜建造风力发电厂的位置（颜色越深，风力越强）：

  存在的问题： 将图像放大之后，会发现有很多大的多边形，如果将大的多边形合并为小的多边形，效果会更好一些，最快的方法是使用 UnionCascaded() 方法。

  6.将小的多边形合并成大的多边形：

import os
from osgeo import ogr
from ospybook.vectorplotter import VectorPlotter

folder = r'E:\osgeopy data\California'

ds = ogr.Open(folder, True)
in_lyr = ds.GetLayerByName('wind_farm')
out_lyr = ds.CreateLayer(
    'wind_farm2', in_lyr.GetSpatialRef(), ogr.wkbPolygon)
out_row = ogr.Feature(out_lyr.GetLayerDefn())

# 创建一个复合多边形来容纳要合并的小多边形
multipoly = ogr.Geometry(ogr.wkbMultiPolygon)

# 循环遍历原始输出中的行并获得几何图形
for in_row in in_lyr:
    in_geom = in_row.geometry().Clone()
    in_geom_type = in_geom.GetGeometryType()

    # 如果几何图形是一个多边形，继续并将其添加到复合多边形
    if in_geom_type == ogr.wkbPolygon:
        multipoly.AddGeometry(in_geom)

    # 打散复合多边形
    elif in_geom_type == ogr.wkbMultiPolygon:
        for i in range(in_geom.GetGeometryCount()):
            multipoly.AddGeometry(
                in_geom.GetGeometryRef(i))

# 将所有要素合并
multipoly = multipoly.UnionCascaded()

# 只保留大的多边形
for i in range(multipoly.GetGeometryCount()):
    poly = multipoly.GetGeometryRef(i)
    if poly.GetArea() > 1000000:
        out_row.SetGeometry(poly)
        out_lyr.CreateFeature(out_row)
        
del ds

  保留了面积至少为1平方千米的区域，只保留大的多边形，使得处理更加容易。

结果对比：