IDL/ENVI+python+shp之获取行政区边界案例

15195775117

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

15195775117

获取行政区边界(python+高德地图)

源代码：


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之调用高德静态地图

15195775117

区分海岸线和行政区边界

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

接下来，为了便于同学制图，我得把行政区边界数据写进shp里；

平时只是读nc或shp，却没写过标准格式的文件，毕竟自己不生产数据，也很少分享数据，处于只进不出的状态；

shp文件怎么写，还得研究下，后续会更新在楼下。

15195775117

寻找读写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'}

看样子有眉目了

15195775117

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

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

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

15195775117

fiona读取中国行政区边界

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

15195775117

小试牛刀，先写一个简单的三角形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文件在ENVI加载后，ENVI会保留它，以至于删不掉shp文件，可以通过关闭所有文件来释放内存：

15195775117

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

15195775117

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

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

代码：
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个文件见附件

ANDYKYLE

给大佬点赞      开发的一把好手