登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
python读取shp地图文件
为了小白画图初体验四，得先做一些铺垫，也是自己踩的坑，然后学的习。其中一个基础就是怎么读取shp文件。
shp文件一般有三个，.shp.dbf.shx，存放的就是地图信息，更直接一点的就是点的信息，由点构成线，由线构成多边形。而读取出来的点信息是以经纬度表示的。
气象绘图能涉及的一般是shapefile库和Basemap库，两个都可以读取shp文件。其实Basemap还是依赖shapefile
这里主要介绍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 -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内部有很多个，记得做循环，筛选你需要的
以我自己画的图结束吧 还是建议大家前往CSDN论坛看代码，能舒服很多：htps://blog.csdn.net/weixin_42372313/article/details/113728014
| 
[复制链接]
窥视卡
雷达卡
发表于 2021-2-6 22:09:03
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
{:eb502:}{:eb502:}{:eb502:}