请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 14328|回复: 15

[混合编程] IDL/ENVI+python+shp之获取行政区边界案例

[复制链接]

新浪微博达人勋

发表于 2020-8-25 16:07:50 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
昨天同学需要画个有市县边界的海南岛地图,他从气象家园下载了一组shp文件(见附件),是中国行政分区图。我原本想用IDL编程读取shp文件,写起来有点费事,于是我想起了ENVI应该可以直接加载。
我打开ENVI,把shp文件拖进去,果然直接就出图了,海南的分区十分粗糙,只给了2个市的边界。
既然shp文件里没有海南岛更多分区的边界数据,我计划从高德地图接口获取数据,见下一楼。

ENVI加载shp.png 海南岛.png



diquJie_polyline.dbf

736.65 KB, 下载次数: 2, 下载积分: 金钱 -5

diquJie_polyline.prj

151 Bytes, 下载次数: 3, 下载积分: 金钱 -5

diquJie_polyline.shp

12.96 MB, 下载次数: 23, 下载积分: 金钱 -5

diquJie_polyline.shx

33.76 KB, 下载次数: 3, 下载积分: 金钱 -5

评分

参与人数 1金钱 +20 收起 理由
wet510 + 20 很给力!

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2020-8-25 16:20:38 | 显示全部楼层
获取行政区边界(python+高德地图)


海南行政分区图.png



源代码:


import json
import requests
import matplotlib.pyplot as plt
#设置字体,避免乱码
plt.rcParams['font.family'] = 'Microsoft YaHei'


key='********************'#高德web开发K码,可免费申请
url='https://restapi.amap.com/v3/config/district?key='+key+\
    '&keywords=海南省&subdistrict=3&extensions=all&output=JSON'
#subdistrict=3表示往下检索3级行政区,
#extensions=all表示需要边界,extensions=base不返回边界数据


re=requests.get(url)
x=json.loads(requests.get(url).text)


#检索结果是一个序列,因为检索的地名可能有重复的,
#我第一次检索关键词用“海南”,结果全国有3个“海南”,
#另2个在青海和内蒙,蒙古称湖为海,海南即湖之南


x=x['districts']
n=len(x)#地区个数
x=x[0]#取序列第一个


#撑满画布隐藏坐标需要提前设置,图画好了再设置没用
plt.axes([0,0,1,1])#撑满画布
plt.axis('off')#隐藏坐标轴


#提取海南岛边界
bound=x['polyline']
cut=bound.split('|')#边界以“|”分段,例如该区有飞地或岛屿
partNum=len(cut)#分段数
for i in range(partNum):
    lng=[]
    lat=[]
    cut2=cut.split(';')#分离坐标
    m=len(cut2)
    for j in range(m):
        cut3=cut2[j].split(',')#分离经纬度
        lng.append(float(cut3[0]))
        lat.append(float(cut3[1]))


    plt.plot(lng,lat)#把这段画出来




sub=x['districts']#下一级分区
n=len(sub)#分区个数
for i in range(n):
    # print(sub['name'])#分区名
    #其实按地区码检索更准确,以下偷懒用地名检索:
    url='https://restapi.amap.com/v3/config/district?key='+key+\
        '&keywords='+sub['name']+'&subdistrict=1&extensions=all&output=JSON'
    re=requests.get(url)
    x=json.loads(requests.get(url).text)
    if x['count']!='1':
        print('出现多地同名')
    x=x['districts']
    x=x[0]
    #获取子区中心坐标,用于标记地名:
    centerlnglat=x['center'].split(',')
    centerlng=float(centerlnglat[0])
    centerlat=float(centerlnglat[1])
    #标注地名,由于坐标是文字最左端,所以文字以五角星开头正好能代表标记点位置:
    plt.text(centerlng,centerlat,'☆'+sub['name'],color='r')
    #提取市区县边界:
    bound=x['polyline']
    cut=bound.split('|')
    partNum=len(cut)
    print(x['name']+'的边界有',partNum,'段')
    for i in range(partNum):
        lng=[]
        lat=[]
        cut2=cut.split(';')
        m=len(cut2)
        for j in range(m):
            cut3=cut2[j].split(',')
            lng.append(float(cut3[0]))
            lat.append(float(cut3[1]))   
        plt.plot(lng,lat)
   
#由于只要海南主岛,限定下范围
plt.ylim(18,20.5)
plt.xlim(108,111.5)


plt.show()




番外:


由于使用IDL经验较多,其实高德地图数据用IDL获取也是一样的,都是两三句代码的事,参考我以前的帖子:
IDL之调用高德地图API计算地图坐标距离
IDL与python之获取网页图片
IDL之GPS坐标转高德地图坐标
IDL之高德地图逆地理编码,坐标转地址
IDL之调用高德静态地图

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-8-25 16:28:48 | 显示全部楼层
本帖最后由 15195775117 于 2020-8-25 17:28 编辑

区分海岸线和行政区边界

我画好图给同学看,他说地图的海边有问题,不细致,我比较了下,原来行政边界和陆地边界不一样,下图比较典型,行政区是可以包括水域的。这样的话,之前的shp文件并不是纯行政区边界线,而是包括了海岸线了。把行政区边界和海岸线结合在一起,才是陆地的行政分区地图。
行政区边界与大陆边界.png

接下来,为了便于同学制图,我得把行政区边界数据写进shp里;
平时只是读nc或shp,却没写过标准格式的文件,毕竟自己不生产数据,也很少分享数据,处于只进不出的状态;
shp文件怎么写,还得研究下,后续会更新在楼下。

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-8-26 14:18:02 | 显示全部楼层
寻找读写shp的python包



今天找包来读写shp文件,第一次用的是shapefile包,命令行安装:pip install pyshp,

import shapefile
shpfile=r"C:\Users\Administrator\Desktop\shp_example\diquJie_polyline.shp"
sf = shapefile.Reader(shpfile)

刚一读就报错,说是解码错误:
File "C:\ProgramData\Anaconda3\lib\site-packages\shapefile.py", line 104, in u
    return v.decode(encoding, encodingErrors)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xca in position 4: invalid continuation byte

以下是参考来源:
对python 读取线的shp文件实例详解


此路不通,我另找了个包Fiona,它不能用pip单独安装,因为存在好几个依赖包,所以用anaconda安装。用conda list命令查看了下,我的fiona是1.8.4,对应文档如下:
Fiona 1.8.4 文档


安装好后,查看一楼的shp文件:
import fiona
shpfile=r"C:\Users\Administrator\Desktop\shp_example\diquJie_polyline.shp"
shps = fiona.open(shpfile)
print('文件编码:',shps.encoding)
print(shps.schema)

输出:
文件编码: iso-8859-1
{'properties': OrderedDict([('FID_êD', 'int:9'), ('GB', 'int:9'), ('NAME', 'str:34'), ('NAME_2', 'str:34'), ('COLORID', 'str:50'), ('AREA', 'float:19.11'), ('LEN', 'float:19.11')]), 'geometry': 'LineString'}

看样子有眉目了


密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-8-26 17:06:36 | 显示全部楼层
本帖最后由 15195775117 于 2020-10-13 17:27 编辑

欲写先读---读取美国州界案例


3个附件为美国州界shp文件,是从IDL的软件数据里找的,画的图和代码如下。
观察数据特点,如果一个州的边界是一笔画的,其边界几何类型就是Polygon;如果是多段,就是MultiPolygon。
边界坐标数据有的是形如(经度、纬度)的元组组成的list,有的是单个的形如(经度、纬度)的元组,个人感觉有些不统一。
美国州界.png
import fiona
import matplotlib.pyplot as plt


shpfile=r"C:\Users\Administrator\Desktop\shp_example\states.shp"
c = fiona.open(shpfile)
print('文件编码:',c.encoding)
print(c.schema['properties'])


# 输出数据的基本信息
print(f'数据范围:{c.bounds}')
print(f'投影定义:{c.crs}')
print(f'数据格式:{c.driver}')
print(f'数据编码:{c.encoding}')
# 输出文件的属性字段信息
fields = c.schema['properties']
print('文件的属性字段信息:')
for k, v in fields.items():
    print(f'{k} -> {v}')
plt.axes([0,0,1,1])#撑满画布
plt.axis('off')#隐藏坐标轴
# 遍历集合中的要素
# f是一个tuple,第一个元素是要素编号,第二个是dict格式的要素
for f in c.items():
    # 输入要素的详细信息
    # 要素是以GeoJSON表示的
    x=f[1]
    # dict_keys(['type', 'id', 'properties', 'geometry'])
    print('州名:',x['properties']['STATE_NAME'])
    geom=x['geometry']
    crd=geom['coordinates']
    n=len(crd)
    print('该州几何体类型与个数:',geom['type'],n)
    for i in range(n):
        y=crd
        m=len(y)
        if m==1:
            y=y[0]
            if isinstance(y,list)==True:
                lng=[]
                lat=[]
                dotnum=len(y)
                for (lng0,lat0) in y:
                    lng.append(lng0)
                    lat.append(lat0)
                plt.plot(lng,lat,lw=0.5)            
        if m>1:
            lng=[]
            lat=[]
            for j in range(m):
                if isinstance(y[j],tuple)==True:
                    lng.append(y[j][0])
                    lat.append(y[j][1])
            plt.plot(lng,lat,lw=0.5)
plt.show()



输出结果:


文件编码: iso-8859-1
OrderedDict([('AREA', 'float:12.3'), ('STATE_NAME', 'str:25'), ('STATE_FIPS', 'str:2'), ('SUB_REGION', 'str:7'), ('STATE_ABBR', 'str:2'), ('POP1990', 'int:10'), ('POP1996', 'int:10')])
数据范围:(-178.21502685546875, 18.924781799316406, -66.9698486328125, 71.40664672851562)
投影定义:{}
数据格式:ESRI Shapefile
数据编码:iso-8859-1
文件的属性字段信息:
AREA -> float:12.3
STATE_NAME -> str:25
STATE_FIPS -> str:2
SUB_REGION -> str:7
STATE_ABBR -> str:2
POP1990 -> int:10
POP1996 -> int:10


州名: Washington
该州几何体类型与个数: MultiPolygon 3
州名: Montana
该州几何体类型与个数: Polygon 1
州名: Maine
该州几何体类型与个数: MultiPolygon 2
... ...
州名: Alaska
该州几何体类型与个数: MultiPolygon 32
州名: Michigan
该州几何体类型与个数: MultiPolygon 5


states.dbf

3.69 KB, 下载次数: 1, 下载积分: 金钱 -5

states.shp

217.13 KB, 下载次数: 1, 下载积分: 金钱 -5

states.shx

508 Bytes, 下载次数: 1, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-8-26 17:39:57 | 显示全部楼层
本帖最后由 15195775117 于 2020-8-26 17:40 编辑

fiona读取中国行政区边界
中国省市界.png
import fiona
import matplotlib.pyplot as plt

shpfile=r"C:\Users\Administrator\Desktop\shp_example\diquJie_polyline.shp"
c = fiona.open(shpfile)
print('文件编码:',c.encoding)
print(c.schema['properties'])

# 输出数据的基本信息
print(f'数据范围:{c.bounds}')
print(f'投影定义:{c.crs}')
print(f'数据格式:{c.driver}')
print(f'数据编码:{c.encoding}')

# 输出文件的属性字段信息
fields = c.schema['properties']
print('文件的属性字段信息:')
for k, v in fields.items():
    print(f'{k} -> {v}')
print('')
plt.axes([0,0,1,1])#撑满画布
plt.axis('off')#隐藏坐标轴
# 遍历集合中的要素
# f是一个tuple,第一个元素是要素编号,第二个是dict格式的要素
for f in c.items():
    # 输入要素的详细信息
    # 要素是以GeoJSON表示的
    x=f[1]
    #字典x的键:'type', 'id', 'properties', 'geometry'
    # print(x['properties'].values())#乱码不管
    geom=x['geometry']
    crd=geom['coordinates']
    n=len(crd)
    # print('几何体类型与个数:',geom['type'],n)
    if geom['type']!='LineString':
        print('非LineString')
    lng=[]
    lat=[]
    for i in range(n):
        y=crd
        lng.append(y[0])
        lat.append(y[1])
    plt.plot(lng,lat,lw=0.5)


plt.show()

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-8-27 10:06:42 | 显示全部楼层
本帖最后由 15195775117 于 2020-8-27 10:08 编辑

小试牛刀,先写一个简单的三角形shp
写shp的代码:

import fiona
from shapely.geometry import Polygon, mapping
outshpfile=r"C:\Users\Administrator\Desktop\writeSHP\out.shp"
# schema是一个字典结构,指定了geometry及其它属性结构
schema = {'geometry': 'Polygon',
          'properties': {'id': 'int', 'name': 'str'}}
ObjShp=fiona.open(outshpfile, mode='w', driver='ESRI Shapefile',
                schema=schema, crs='EPSG:4326', encoding='utf-8')
coordinates = [(117.4219, 40.21),
                (116.8066, 39.5947),
                (117.2461, 40.5176),
                (117.4219, 40.21)]#这里坐标用[xx,xx]也行,不限定为元祖格式
polygon = Polygon(coordinates)  # 使用地理坐标定义Polygon对象
polygon = mapping(polygon)  # 将Polygon对象转为GeoJSON格式
feature = {'geometry': polygon,
            'properties': {'id': 1, 'name': '三角形'}}
ObjShp.write(feature)
ObjShp.close()
结果:
生成了5个文件:cpg,dbf,prj,shp,shx
把shp文件拖入ENVI,顺利显示:
写shp.png
小技巧:
shp文件在ENVI加载后,ENVI会保留它,以至于删不掉shp文件,可以通过关闭所有文件来释放内存:
ENVI清理内存.png




密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-8-27 10:36:10 | 显示全部楼层
shp写入多个几何体
2个三角形.png
import fiona
from shapely.geometry import Polygon, mapping
outshpfile=r"C:\Users\Administrator\Desktop\writeSHP\out.shp"
# schema是一个字典结构,指定了geometry及其它属性结构
schema = {'geometry': 'Polygon',
          'properties': {'id': 'int', 'name': 'str'}}
ObjShp=fiona.open(outshpfile, mode='w', driver='ESRI Shapefile',
                schema=schema, crs='EPSG:4326', encoding='utf-8')

#第一个几何体:
coordinates = [(117.4219, 40.21),
                (116.8066, 39.5947),
                (117.2461, 40.5176),
                (117.4219, 40.21)]
polygon = Polygon(coordinates)  # 使用地理坐标定义Polygon对象
polygon = mapping(polygon)  # 将Polygon对象转为GeoJSON格式
feature = {'geometry': polygon,
            'properties': {'id': 1, 'name': '三角形'}}
ObjShp.write(feature)
#第二个几何体:

coordinates = [(117.7319, 40.92),
                (116.0166, 39.0047),
                (117.9561, 40.0276),
                (117.7319, 40.92)]
polygon = Polygon(coordinates)  # 使用地理坐标定义Polygon对象
polygon = mapping(polygon)  # 将Polygon对象转为GeoJSON格式
feature = {'geometry': polygon,
            'properties': {'id': 2, 'name': '三角形2'}}
ObjShp.write(feature)

ObjShp.close()



密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-8-27 11:37:38 | 显示全部楼层
本帖最后由 15195775117 于 2020-8-27 11:42 编辑

收官:获取行政区边界并写入shp

经过以上几楼的摸索铺垫,写出同学需要的shp已经水到渠成了。
红色是海岸线,是一楼的shp里的,绿色的是行政区边界,从高德地图扒的,行政区边界比海岸线大一圈


1977908067.jpg




代码:
import json
import fiona
import requests
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Polygon, mapping
#设置字体,避免乱码
plt.rcParams['font.family'] = 'Microsoft YaHei'
#需要写入的文件:
outshpfile=r"C:\Users\Administrator\Desktop\writeHaiNanShp\out.shp"
schema = {'geometry': 'Polygon',
          'properties': {'id': 'int', 'name': 'str'}}
ObjShp=fiona.open(outshpfile, mode='w', driver='ESRI Shapefile',
                schema=schema, crs='EPSG:4326', encoding='utf-8')
geomNum=1#记录连续性几何体的个数
key='*****************************'#高德web开发K码
url='https://restapi.amap.com/v3/config/district?key='+key+\
    '&keywords=海南省&subdistrict=3&extensions=all&output=JSON'
#subdistrict=3表示往下检索3个行政区,
#extensions=all表示需要边界,extensions=base不返回边界数据
re=requests.get(url)
x=json.loads(requests.get(url).text)
#检索结果是一个序列,因为检索的地名可能有重复的,
#我第一次检索关键词用“海南”,结果全国有3个“海南”,
#另2个在青海和内蒙,蒙古称湖为海,海南即湖之南
x=x['districts']
n=len(x)#地区个数
x=x[0]#取序列第一个
#提取海南岛边界
bound=x['polyline']
cut=bound.split('|')#边界以“|”分段,例如该区有飞地或岛屿
partNum=len(cut)#分段数
print('分段数=',partNum)
for i in range(partNum):
    coordinates=[0.0,0.0]
    cut2=cut.split(';')#分离坐标
    m=len(cut2)
    for j in range(m):
        cut3=cut2[j].split(',')#分离经纬度
        coordinates=np.vstack((coordinates,[float(cut3[0]),float(cut3[1])]))

    coordinates=coordinates[1:,:]
    polygon = Polygon(coordinates)  # 使用地理坐标定义Polygon对象
    polygon = mapping(polygon)  # 将Polygon对象转为GeoJSON格式
    feature = {'geometry': polygon,
                'properties': {'id':geomNum, 'name': '边界'}}
    ObjShp.write(feature)
    geomNum=geomNum+1

sub=x['districts']#下一级分区
n=len(sub)#分区个数
for i in range(n):
    # print(sub['name'])#分区名
    #其实按地区码检索更准确,以下偷懒用地名检索:
    url='https://restapi.amap.com/v3/config/district?key='+key+\
        '&keywords='+sub['name']+'&subdistrict=1&extensions=all&output=JSON'
    re=requests.get(url)
    x=json.loads(requests.get(url).text)
    if x['count']!='1':
        print('出现多地同名')
    x=x['districts']
    x=x[0]
    #提取市区县边界:
    bound=x['polyline']
    cut=bound.split('|')
    partNum=len(cut)
    print(x['name']+'的边界有',partNum,'段')
    for i in range(partNum):
        coordinates=[0.0,0.0]
        cut2=cut.split(';')
        m=len(cut2)
        for j in range(m):
            cut3=cut2[j].split(',')
            coordinates=np.vstack((coordinates,[float(cut3[0]),float(cut3[1])]))
        coordinates=coordinates[1:,:]
        polygon = Polygon(coordinates)  # 使用地理坐标定义Polygon对象
        polygon = mapping(polygon)  # 将Polygon对象转为GeoJSON格式
        feature = {'geometry': polygon,
                    'properties': {'id':geomNum, 'name': '边界'}}
        ObjShp.write(feature)
        geomNum=geomNum+1        
ObjShp.close()
print('输出结束')
生成5个文件见附件

海南岛行政区边界.rar

265.43 KB, 下载次数: 1, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-8-27 15:16:36 | 显示全部楼层
给大佬点赞      开发的一把好手
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表