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

15195775117

ANDYKYLE 发表于 2020-8-27 15:16
给大佬点赞      开发的一把好手

多谢，希望能给有需要的人一点启发

15195775117

后续：制作无海南岛的shp

由于一楼的全国shp图中，海岸线与行政区边界并存，他想把其中的海南岛部分剔除，用我制作的海南岛行政地图shp替换那一块区域。
读取的数据在附件cnmap.rar；
以下代码是读写单个文件的，可以自己改下读写的文件名，读完3个文件后，最后输出的结果在“无海南岛.rar”

代码：

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 27 15:15:34 2020
制作没有海南岛的shp
@author: Administrator
"""

import fiona
import numpy as np
from shapely.geometry import LineString, mapping

#输出的shp文件：
outshpfile=r"C:\Users\Administrator\Desktop\NoHaiNan\diquJie_polyline无海南岛.shp"
schema = {'geometry': 'LineString',#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#对连续的几何体进行计数

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

# 输出文件的属性字段信息
fields = c.schema['properties']
print('文件的属性字段信息：')
for k, v in fields.items():
    print(f'{k} -> {v}')

for f in c.items():
    # 输入要素的详细信息
    # 要素是以GeoJSON表示的
    x=f[1]
    #字典x的键:'type', 'id', 'properties', 'geometry'
    geom=x['geometry']
    crd=geom['coordinates']
    n=len(crd)
    # print('几何体类型与点个数:',geom['type'],n)
    if geom['type']!='LineString':
        print('非LineString')
    coordinates=[0.0,0.0]
    dotNum=0
    for i in range(n):
        y=crd
                #海南岛范围:
        if y[0]>108.490 and y[0]<111.064 and y[1]>18.149 and y[1]<20.172:
            continue
        coordinates=np.vstack((coordinates,[y[0],y[1]]))
        dotNum=dotNum+1
    if dotNum==0:
        continue
    coordinates=coordinates[1:,:]
    lineString = LineString(coordinates)  # 使用地理坐标定义LineString对象
    lineString = mapping(lineString)  # 将LineString对象转为GeoJSON格式
    feature = {'geometry': lineString,
                'properties': {'id':geomNum, 'name': '边界'}}
    ObjShp.write(feature)
    geomNum=geomNum+1

ObjShp.close()

15195775117

突然想到的瑕疵
高德地图提供的行政区边界坐标应该是高德坐标系的，而一楼的shp按理是GPS坐标；
高德和GPS坐标经常差个几百米，但是在该案例的空间尺度下几乎是看不出来的

yusw

看到你的程序，真是叹为观止！请问可以将中国地图的shp放在idl里，然后用mapcontinents调用吗

15195775117

yusw 发表于 2020-9-9 07:48
看到你的程序，真是叹为观止！请问可以将中国地图的shp放在idl里，然后用mapcontinents调用吗

过誉了，哪有那么夸张！都是些基础的
IDL调用当然可以，你打开智能化工具窗口，就可以直接打开shp文件：

yusw

15195775117 发表于 2020-9-9 09:16
过誉了，哪有那么夸张！都是些基础的
IDL调用当然可以，你打开智能化工具窗口，就可以直接打开shp文件： ...

多谢楼主！