定义键盘快捷键 SLAM file delphi optimization drupal ansible coldfusion vue前端 jq选择第一个子元素 sublime分屏快捷键 html好看的字体样式 mac版的matlab好用吗 phpstorm插件 python类 python指数函数 python中文 python编程语言 python如何定义变量 java数组反转 java数据库 安装java环境 java方法重载 java将数据写入文件 变量的类型 java系统学习 pascal教程 flash相册制作 相关软件 eclipse中文版下载 虚拟打印机安装 模拟人生2夜生活 抖音代码 羽毛球拍握法 fireworks下载 quickchm 粉碎文件工具 苹果手机怎么看内存 起义任务线 mysql退出命令
当前位置: 首页 > 学习教程  > python

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

2021/2/7 9:55:34 文章标签: 测试文章如有侵权请发送至邮箱809451989@qq.com投诉后文章立即删除

区域的合适性从1到7,一般数值大于等于3的区域是合适的。将这些区域与人口普查数据结合,选出风速合适并且每平方千米人口不到0.5%的地方。 人口普查数据包含每个人口普查单元的人口数据,但没有人口密度数据。   1. 要获得人口属性字段&…

  区域的合适性从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. 对风力数据集进行属性过滤,只显示等级为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')

  4. 创建一个新的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())

创建结果:

在这里插入图片描述
  5. 将人口普查数据与风力数据相交,并储存到新创建的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:\Google chrome\Download\gis with python\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平方千米的区域,只保留大的多边形,使得处理更加容易。

结果对比:
在这里插入图片描述
在这里插入图片描述


本文链接: http://www.dtmao.cc/news_show_2000207.shtml

附件下载

相关教程

    暂无相关的数据...

共有条评论 网友评论

验证码: 看不清楚?