分库查询 二代征信 silverlight background postman 网站后台管理模板 pmp学习视频 jq延时 excel动态图表制作 java运行软件 less的比较级 matlab中不等于怎么表示 js数组截取前5个 python计算器 javalabel java使用mysql java中的正则表达式 java入门课程 如何查看java版本 java终止线程 java数组追加 linuxshell linux内核编程 pr滤镜插件 淘宝自动发货软件 疯狂java 小洛快跑 dnf瞎子传说套选择 ps3d字体 快手规则 深入解析windows操作系统 文件压缩工具 脚本模板 pip安装教程 迅雷共享会员 伏魔战记隐藏英雄 mywi iphone常去地点怎么查看 电视应用安装器 图片轮播
当前位置: 首页 > 学习教程  > python

python读取shp文件,

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

python读取shp地图文件 为了小白画图初体验四,得先做一些铺垫,也是自己踩的坑,然后学的习。其中一个基础就是怎么读取shp文件。 shp文件一般有三个,.shp.dbf.shx,存放的就是地图信息,更直接一点的就是点的…

python读取shp地图文件

为了小白画图初体验四,得先做一些铺垫,也是自己踩的坑,然后学的习。其中一个基础就是怎么读取shp文件。
shp文件一般有三个,.shp.dbf.shx,存放的就是地图信息,更直接一点的就是点的信息,由点构成线,由线构成多边形。而读取出来的点信息是以经纬度表示的。
气象绘图能涉及的一般是shapefile库和Basemap库,两个都可以读取shp文件。
这里主要介绍shapefile库读取的方法,和Basemap库出现“utf-8”错误的解决方法。

shapefile库读取shp文件

概述

import shapefile
sf = shapefile.Reader('Zhejiang_province')
shapes = sf.shapes()

这是读取shp文件的基本操作,shapes里面存放的就是地图信息类的数组。,里面一个或者多个shapefile类

sf = shapefile.Reader('Zhejiang_province')
shapes = sf.shapes()
pts = shapes[0].points
prt =shapes[0].parts

上图的文件是浙江省省界,就一个文件

sf = shapefile.Reader("country1")
for shape_rec in sf.shapeRecords():
    pts = shape_rec.shape.points#边界点
    prt = shape_rec.shape.parts

上图的文件是国家界,就多个类,每一个存放了每个国界的信息
下面以单一的文件给大家作介绍,方便理解,多文件无非就是写个循环,加个判断。

简单介绍

浙江省省界文件

sf = shapefile.Reader('Zhejiang_province')
shapes = sf.shapes()
pts = shapes[0].points
prt =shapes[0].parts

.points里面存放所有点的信息,每个点用(x,y)表示,其实就是所表示区域的边界点,代码里读出的就是浙江省的所有边界点,用经纬度来表示位置。这里存在一个问题,那就是以浙江省为例,有舟山群岛,所有点的坐标不连续,其实相当于许多小区域的组合,但是.points里面不管这些,只会表示点,咱们画图来看一下

import shapefile
import matplotlib.pyplot as plt

sf = shapefile.Reader('Zhejiang_province')
shapes = sf.shapes()
codes = []
pts = shapes[0].points

x,y = zip(*pts)#把经纬度分别给到x,y
fig = plt.figure(figsize=[12,18])
ax = fig.add_subplot(111)
ax.plot(x, y, '-', lw=1, color='k')
plt.show()

在这里插入图片描述
可以看出,画完陆地去画岛的时候很生硬,直接扯过去了(因为用的是plot,如果是散点图不会画成这样,但我们不需要散点)。
这里涉及了另一个知识点

shapefile.Reader类属性

包含有4个属性:数据类型(shapeType),代表该"几何数据"对象的数据类型(点,shapeType=1,线,shapeType=3,多边形,shapeType=5);数据范围(bbox),只针对多点数据,代表该"几何数据"对象的边界范围;数据块(parts),只针对线或者多边形,代表该"几何数据"对象各个块的第一个点的索引;点集(points),代表该"几何数据"对象的所有点坐标。

点集(points):这个已经介绍过了

数据类型(shapeType):是shp文件内部按什么状态存放,以上面浙江省界例子来看是多边形shapeType=5,即区块

sf = shapefile.Reader('Zhejiang_province')
shapes = sf.shapes()
print(shapes[0].shapeType)

数据范围(bbox):就是经纬度的最大最小值,四个,上下左右
重点说一下
数据块(parts):既然上文提到了shp文件所描绘的边界会出现断开的现象,那从哪开始断开呢,就在这个属性里面记录。parts属性是一个列表,记录了每个区块第一个点的索引值,即下标。
气象绘图目前精确一点的mask方法就是利用parts,points结合matplotlib中的路径功能来做mask。

气象绘图利用shp文件做掩膜

先放代码,还是以浙江省界为例

fig = plt.figure(figsize=[12,18])
ax = fig.add_subplot(111)

sf = shapefile.Reader('Zhejiang_province')
shapes = sf.shapes()
codes = []
pts = shapes[0].points#边界点
prt = list(shapes[0].parts) + [len(pts)]#区块起始索引
for i in range(len(prt) - 1):
    codes += [Path.MOVETO]#点移动
    codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)#画线
    codes += [Path.CLOSEPOLY]#这块画完,循环结束,下一块
clip = Path(pts, codes)#利用数据和路径生成一个画图动作
clip = PathPatch(clip, transform=ax.transData)#再加入ax的变换

prt的含义:上文已经说过了parts的意义,设想一下,两块区域的的起始坐标中间夹的不就是一块区域的所有坐标吗,但是最后一块区域没有结束值,所有加了[len(pts)],这就是最后一个点的索引。

for循环的作用:建立一个绘图路径,利用区块起始点的索引生成一个画图动作

最后生成的clip就是我们需要的掩膜mask,在合适的位置配合下列代码完成掩膜白化效果:

cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')
for contour in cs.collections:
        contour.set_clip_path(clip)

大家可以不用理解原理,只需要明白,怎么读取shp文件,找到你需要的边界坐标点(points)和(parts),把他们带入其中,执行这个过程就可以了。

下一个问题:Basemap库出现“utf-8”错误的解决方法

Basemap库还是很好用的,但是有些朋友肯定出现过这个问题

m.readshapefile('Zhejiang_province','Zhejiang Map',color='k',linewidth=1.2)

当我们执行这条命令时报错

  File "C:\ProgramData\Miniconda3\lib\site-packages\shapefile.py", line 104, in u
    return v.decode(encoding, encodingErrors)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xe8 in position 2: invalid continuation byte

这是shapefile需要读取文件得“utf-8”编码的缘故,解决方法,进入一下路径(作者的电脑,大家相应寻找)C:\ProgramData\Miniconda3\Lib\site-packages
找到shapefile.py文件,复制一份(防止操作失误),然后在shapefile.py文件的第104行做修改

#原来的
return v.decode(encoding, encodingErrors)
#修改为
return v.decode('latin-1', encodingErrors)

这样编译就可以通过了,如果以后碰到utf-8的出错了记得改回来。

最后还是提醒大家文章是以单个shapes[0].points为例的,有一部分shp内部有很多个,记得做循环,筛选你需要的
以我自己画的图结束吧
在这里插入图片描述


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

附件下载

上一篇:第四扩展FS Writeup By K龙

下一篇:opencv02

相关教程

    暂无相关的数据...

共有条评论 网友评论

验证码: 看不清楚?