请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

新浪微博登陆

只需一步, 快速开始

QQ登录

只需一步,快速开始

搜索
查看: 202|回复: 5

[源代码] matplotlib绘制基于经纬度的任意shp文件

[复制链接] |关注本帖

新浪微博达人勋

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 非对称 于 2018-2-13 18:14 编辑
  1. # -*- coding: utf-8 -*-
  2. """
  3. 特此申明,以下源代码改写自basemap的readshapefile
  4. 适合在无投影下的经纬度坐标中绘制。
  5. 需要在投影坐标系绘制,请自行用basemap的readshapefile绘制
  6. -------------------------------------------------
  7.    File Name:     readshp
  8.    Description :
  9.    Author :        Xy Liu
  10.    date:          2018/2/13
  11. -------------------------------------------------
  12.    Change Activity:
  13.                    2018/2/13
  14. -------------------------------------------------
  15. """
  16. import os

  17. from matplotlib import dedent
  18. from matplotlib.collections import LineCollection

  19. __author__ = 'Xy Liu'

  20. def readshapefile(shapefile, name, drawbounds=True, zorder=None,
  21.                   linewidth=0.5, color='k', antialiased=1, ax=None,
  22.                   default_encoding='utf-8'):
  23.     """
  24.     Read in shape file, optionally draw boundaries on map.

  25.     .. note::
  26.       - Assumes shapes are 2D
  27.       - only works for Point, MultiPoint, Polyline and Polygon shapes.
  28.       - vertices/points must be in geographic (lat/lon) coordinates.

  29.     Mandatory Arguments:

  30.     .. tabularcolumns:: |l|L|

  31.     ==============   ====================================================
  32.     Argument         Description
  33.     ==============   ====================================================
  34.     shapefile        path to shapefile components.  Example:
  35.                      shapefile='/home/jeff/esri/world_borders' assumes
  36.                      that world_borders.shp, world_borders.shx and
  37.                      world_borders.dbf live in /home/jeff/esri.
  38.     name             name for Basemap attribute to hold the shapefile
  39.                      vertices or points in map projection
  40.                      coordinates. Class attribute name+'_info' is a list
  41.                      of dictionaries, one for each shape, containing
  42.                      attributes of each shape from dbf file, For
  43.                      example, if name='counties', self.counties
  44.                      will be a list of x,y vertices for each shape in
  45.                      map projection  coordinates and self.counties_info
  46.                      will be a list of dictionaries with shape
  47.                      attributes.  Rings in individual Polygon
  48.                      shapes are split out into separate polygons, and
  49.                      additional keys 'RINGNUM' and 'SHAPENUM' are added
  50.                      to the shape attribute dictionary.
  51.     ==============   ====================================================

  52.     The following optional keyword arguments are only relevant for Polyline
  53.     and Polygon shape types, for Point and MultiPoint shapes they are
  54.     ignored.

  55.     .. tabularcolumns:: |l|L|

  56.     ==============   ====================================================
  57.     Keyword          Description
  58.     ==============   ====================================================
  59.     drawbounds       draw boundaries of shapes (default True).
  60.     zorder           shape boundary zorder (if not specified,
  61.                      default for mathplotlib.lines.LineCollection
  62.                      is used).
  63.     linewidth        shape boundary line width (default 0.5)
  64.     color            shape boundary line color (default black)
  65.     antialiased      antialiasing switch for shape boundaries
  66.                      (default True).
  67.     ax               axes instance (overrides default axes instance)
  68.     ==============   ====================================================

  69.     A tuple (num_shapes, type, min, max) containing shape file info
  70.     is returned.
  71.     num_shapes is the number of shapes, type is the type code (one of
  72.     the SHPT* constants defined in the shapelib module, see
  73.     http://shapelib.maptools.org/shp_api.html) and min and
  74.     max are 4-element lists with the minimum and maximum values of the
  75.     vertices. If ``drawbounds=True`` a
  76.     matplotlib.patches.LineCollection object is appended to the tuple.
  77.     """
  78.     import shapefile as shp
  79.     from shapefile import Reader
  80.     shp.default_encoding = default_encoding
  81.     if not os.path.exists('%s.shp' % shapefile):
  82.         raise IOError('cannot locate %s.shp' % shapefile)
  83.     if not os.path.exists('%s.shx' % shapefile):
  84.         raise IOError('cannot locate %s.shx' % shapefile)
  85.     if not os.path.exists('%s.dbf' % shapefile):
  86.         raise IOError('cannot locate %s.dbf' % shapefile)
  87.     # open shapefile, read vertices for each object, convert
  88.     # to map projection coordinates (only works for 2D shape types).
  89.     try:
  90.         shf = Reader(shapefile)
  91.     except:
  92.         raise IOError('error reading shapefile %s.shp' % shapefile)
  93.     fields = shf.fields
  94.     coords = []
  95.     attributes = []
  96.     msg = dedent("""
  97.         shapefile must have lat/lon vertices  - it looks like this one has vertices
  98.         in map projection coordinates. You can convert the shapefile to geographic
  99.         coordinates using the shpproj utility from the shapelib tools
  100.         (http://shapelib.maptools.org/shapelib-tools.html)""")
  101.     shptype = shf.shapes()[0].shapeType
  102.     bbox = shf.bbox.tolist()
  103.     info = dict()
  104.     info['info'] = (shf.numRecords, shptype, bbox[0:2] + [0., 0.], bbox[2:] + [0., 0.])
  105.     npoly = 0
  106.     for shprec in shf.shapeRecords():
  107.         shp = shprec.shape
  108.         rec = shprec.record
  109.         npoly = npoly + 1
  110.         if shptype != shp.shapeType:
  111.             raise ValueError('readshapefile can only handle a single shape type per file')
  112.         if shptype not in [1, 3, 5, 8]:
  113.             raise ValueError('readshapefile can only handle 2D shape types')
  114.         verts = shp.points
  115.         if shptype in [1, 8]:  # a Point or MultiPoint shape.
  116.             lons, lats = list(zip(*verts))
  117.             if max(lons) > 721. or min(lons) < -721. or max(lats) > 90.01 or min(lats) < -90.01:
  118.                 raise ValueError(msg)
  119.             # if latitude is slightly greater than 90, truncate to 90
  120.             lats = [max(min(lat, 90.0), -90.0) for lat in lats]
  121.             if len(verts) > 1:  # MultiPoint
  122.                 x, y = lons, lats
  123.                 coords.append(list(zip(x, y)))
  124.             else:  # single Point
  125.                 x, y = lons[0], lats[0]
  126.                 coords.append((x, y))
  127.             attdict = {}
  128.             for r, key in zip(rec, fields[1:]):
  129.                 attdict[key[0]] = r
  130.             attributes.append(attdict)
  131.         else:  # a Polyline or Polygon shape.
  132.             parts = shp.parts.tolist()
  133.             ringnum = 0
  134.             for indx1, indx2 in zip(parts, parts[1:] + [len(verts)]):
  135.                 ringnum = ringnum + 1
  136.                 lons, lats = list(zip(*verts[indx1:indx2]))
  137.                 if max(lons) > 721. or min(lons) < -721. or max(lats) > 90.01 or min(lats) < -90.01:
  138.                     raise ValueError(msg)
  139.                 # if latitude is slightly greater than 90, truncate to 90
  140.                 lats = [max(min(lat, 90.0), -90.0) for lat in lats]
  141.                 x, y = lons, lats
  142.                 coords.append(list(zip(x, y)))
  143.                 attdict = {}
  144.                 for r, key in zip(rec, fields[1:]):
  145.                     attdict[key[0]] = r
  146.                 # add information about ring number to dictionary.
  147.                 attdict['RINGNUM'] = ringnum
  148.                 attdict['SHAPENUM'] = npoly
  149.                 attributes.append(attdict)
  150.     # draw shape boundaries for polylines, polygons  using LineCollection.
  151.     if shptype not in [1, 8] and drawbounds:
  152.         # get current axes instance (if none specified).
  153.         import matplotlib.pyplot as plt
  154.         ax = ax or plt.gca()
  155.         # make LineCollections for each polygon.
  156.         lines = LineCollection(coords, antialiaseds=(antialiased,))
  157.         lines.set_color(color)
  158.         lines.set_linewidth(linewidth)
  159.         lines.set_label('_nolabel_')
  160.         if zorder is not None:
  161.             lines.set_zorder(zorder)
  162.         ax.add_collection(lines)
  163.     # set axes limits to fit map region.
  164.     # self.set_axes_limits(ax=ax)
  165.     # # clip boundaries to map limbs
  166.     # lines,c = self._cliplimb(ax,lines)
  167.     # info = info + (lines,)
  168.     info[name] = coords
  169.     info[name + '_info'] = attributes
  170.     return info


  171. if __name__ == '__main__':
  172.     file = 'test.shp'
  173.     readshapefile(file.replace('.shp', ''), os.path.basename(file), color='k', linewidth=1)
复制代码

18021308.216.png

评分

参与人数 1金钱 +22 贡献 +4 收起 理由
又是那隻貓 + 22 + 4

查看全部评分

密码修改失败请联系qq:937062711

新浪微博达人勋

发表于 2018-2-13 09:59:05 | 显示全部楼层 |取消关注该作者的回复
无投影和有投影,可否上个图或者举个例子
密码修改失败请联系qq:937062711

新浪微博达人勋

 楼主| 发表于 2018-2-13 10:03:19 | 显示全部楼层 |取消关注该作者的回复
平流层的萝卜 发表于 2018-2-13 09:59
无投影和有投影,可否上个图或者举个例子

这段代码 适合在自定义坐标下绘制shp文件,取决于你的shp文件的坐标系统和自定义坐标系统匹配
密码修改失败请联系qq:937062711

新浪微博达人勋

发表于 2018-2-13 10:18:39 | 显示全部楼层 |取消关注该作者的回复
感谢分享~~快过年了还没闲着,新年快乐~~
密码修改失败请联系qq:937062711

新浪微博达人勋

发表于 2018-2-13 10:44:35 | 显示全部楼层 |取消关注该作者的回复
大神出来分享东西了啊~

顶一个
密码修改失败请联系qq:937062711

新浪微博达人勋

发表于 2018-2-14 14:48:59 | 显示全部楼层 |取消关注该作者的回复
新年快乐~~
密码修改失败请联系qq:937062711
回复

使用道具 举报

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

本版积分规则

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

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

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