爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
12
返回列表 发新帖
楼主: 15195775117

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

[复制链接]

新浪微博达人勋

 楼主| 发表于 2020-8-27 16:54:42 | 显示全部楼层
ANDYKYLE 发表于 2020-8-27 15:16
给大佬点赞      开发的一把好手

多谢,希望能给有需要的人一点启发
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-8-27 17:05:36 | 显示全部楼层
后续:制作无海南岛的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()

cnmap.rar

6.4 MB, 下载次数: 0, 下载积分: 金钱 -5

无海南岛.rar

6.34 MB, 下载次数: 0, 下载积分: 金钱 -5

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

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-8-28 08:50:58 | 显示全部楼层
突然想到的瑕疵
高德地图提供的行政区边界坐标应该是高德坐标系的,而一楼的shp按理是GPS坐标;
高德和GPS坐标经常差个几百米,但是在该案例的空间尺度下几乎是看不出来的
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-9-9 07:48:03 | 显示全部楼层
看到你的程序,真是叹为观止!请问可以将中国地图的shp放在idl里,然后用mapcontinents调用吗
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-9-9 09:16:43 | 显示全部楼层
yusw 发表于 2020-9-9 07:48
看到你的程序,真是叹为观止!请问可以将中国地图的shp放在idl里,然后用mapcontinents调用吗

过誉了,哪有那么夸张!都是些基础的
IDL调用当然可以,你打开智能化工具窗口,就可以直接打开shp文件:
360截图20200909090836176.png 360截图20200909091510362.png
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-9-9 09:26:30 | 显示全部楼层
15195775117 发表于 2020-9-9 09:16
过誉了,哪有那么夸张!都是些基础的
IDL调用当然可以,你打开智能化工具窗口,就可以直接打开shp文件: ...

多谢楼主!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

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