matplotlib绘制基于经纬度的任意shp文件

非对称

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. 
    
18. from matplotlib import dedent
    
19. from matplotlib.collections import LineCollection
    
20. 
    
21. __author__ = 'Xy Liu'
    
22. 
    
23. def readshapefile(shapefile, name, drawbounds=True, zorder=None,
    
24.                   linewidth=0.5, color='k', antialiased=1, ax=None,
    
25.                   default_encoding='utf-8'):
    
26.     """
    
27.     Read in shape file, optionally draw boundaries on map.
    
28. 
    
29.     .. note::
    
30.       - Assumes shapes are 2D
    
31.       - only works for Point, MultiPoint, Polyline and Polygon shapes.
    
32.       - vertices/points must be in geographic (lat/lon) coordinates.
    
33. 
    
34.     Mandatory Arguments:
    
35. 
    
36.     .. tabularcolumns:: |l|L|
    
37. 
    
38.     ==============   ====================================================
    
39.     Argument         Description
    
40.     ==============   ====================================================
    
41.     shapefile        path to shapefile components.  Example:
    
42.                      shapefile='/home/jeff/esri/world_borders' assumes
    
43.                      that world_borders.shp, world_borders.shx and
    
44.                      world_borders.dbf live in /home/jeff/esri.
    
45.     name             name for Basemap attribute to hold the shapefile
    
46.                      vertices or points in map projection
    
47.                      coordinates. Class attribute name+'_info' is a list
    
48.                      of dictionaries, one for each shape, containing
    
49.                      attributes of each shape from dbf file, For
    
50.                      example, if name='counties', self.counties
    
51.                      will be a list of x,y vertices for each shape in
    
52.                      map projection  coordinates and self.counties_info
    
53.                      will be a list of dictionaries with shape
    
54.                      attributes.  Rings in individual Polygon
    
55.                      shapes are split out into separate polygons, and
    
56.                      additional keys 'RINGNUM' and 'SHAPENUM' are added
    
57.                      to the shape attribute dictionary.
    
58.     ==============   ====================================================
    
59. 
    
60.     The following optional keyword arguments are only relevant for Polyline
    
61.     and Polygon shape types, for Point and MultiPoint shapes they are
    
62.     ignored.
    
63. 
    
64.     .. tabularcolumns:: |l|L|
    
65. 
    
66.     ==============   ====================================================
    
67.     Keyword          Description
    
68.     ==============   ====================================================
    
69.     drawbounds       draw boundaries of shapes (default True).
    
70.     zorder           shape boundary zorder (if not specified,
    
71.                      default for mathplotlib.lines.LineCollection
    
72.                      is used).
    
73.     linewidth        shape boundary line width (default 0.5)
    
74.     color            shape boundary line color (default black)
    
75.     antialiased      antialiasing switch for shape boundaries
    
76.                      (default True).
    
77.     ax               axes instance (overrides default axes instance)
    
78.     ==============   ====================================================
    
79. 
    
80.     A tuple (num_shapes, type, min, max) containing shape file info
    
81.     is returned.
    
82.     num_shapes is the number of shapes, type is the type code (one of
    
83.     the SHPT* constants defined in the shapelib module, see
    
84.     http://shapelib.maptools.org/shp_api.html) and min and
    
85.     max are 4-element lists with the minimum and maximum values of the
    
86.     vertices. If ``drawbounds=True`` a
    
87.     matplotlib.patches.LineCollection object is appended to the tuple.
    
88.     """
    
89.     import shapefile as shp
    
90.     from shapefile import Reader
    
91.     shp.default_encoding = default_encoding
    
92.     if not os.path.exists('%s.shp' % shapefile):
    
93.         raise IOError('cannot locate %s.shp' % shapefile)
    
94.     if not os.path.exists('%s.shx' % shapefile):
    
95.         raise IOError('cannot locate %s.shx' % shapefile)
    
96.     if not os.path.exists('%s.dbf' % shapefile):
    
97.         raise IOError('cannot locate %s.dbf' % shapefile)
    
98.     # open shapefile, read vertices for each object, convert
    
99.     # to map projection coordinates (only works for 2D shape types).
    
100.     try:
     
101.         shf = Reader(shapefile)
     
102.     except:
     
103.         raise IOError('error reading shapefile %s.shp' % shapefile)
     
104.     fields = shf.fields
     
105.     coords = []
     
106.     attributes = []
     
107.     msg = dedent("""
     
108.         shapefile must have lat/lon vertices  - it looks like this one has vertices
     
109.         in map projection coordinates. You can convert the shapefile to geographic
     
110.         coordinates using the shpproj utility from the shapelib tools
     
111.         (http://shapelib.maptools.org/shapelib-tools.html)""")
     
112.     shptype = shf.shapes()[0].shapeType
     
113.     bbox = shf.bbox.tolist()
     
114.     info = dict()
     
115.     info['info'] = (shf.numRecords, shptype, bbox[0:2] + [0., 0.], bbox[2:] + [0., 0.])
     
116.     npoly = 0
     
117.     for shprec in shf.shapeRecords():
     
118.         shp = shprec.shape
     
119.         rec = shprec.record
     
120.         npoly = npoly + 1
     
121.         if shptype != shp.shapeType:
     
122.             raise ValueError('readshapefile can only handle a single shape type per file')
     
123.         if shptype not in [1, 3, 5, 8]:
     
124.             raise ValueError('readshapefile can only handle 2D shape types')
     
125.         verts = shp.points
     
126.         if shptype in [1, 8]:  # a Point or MultiPoint shape.
     
127.             lons, lats = list(zip(*verts))
     
128.             if max(lons) > 721. or min(lons) < -721. or max(lats) > 90.01 or min(lats) < -90.01:
     
129.                 raise ValueError(msg)
     
130.             # if latitude is slightly greater than 90, truncate to 90
     
131.             lats = [max(min(lat, 90.0), -90.0) for lat in lats]
     
132.             if len(verts) > 1:  # MultiPoint
     
133.                 x, y = lons, lats
     
134.                 coords.append(list(zip(x, y)))
     
135.             else:  # single Point
     
136.                 x, y = lons[0], lats[0]
     
137.                 coords.append((x, y))
     
138.             attdict = {}
     
139.             for r, key in zip(rec, fields[1:]):
     
140.                 attdict[key[0]] = r
     
141.             attributes.append(attdict)
     
142.         else:  # a Polyline or Polygon shape.
     
143.             parts = shp.parts.tolist()
     
144.             ringnum = 0
     
145.             for indx1, indx2 in zip(parts, parts[1:] + [len(verts)]):
     
146.                 ringnum = ringnum + 1
     
147.                 lons, lats = list(zip(*verts[indx1:indx2]))
     
148.                 if max(lons) > 721. or min(lons) < -721. or max(lats) > 90.01 or min(lats) < -90.01:
     
149.                     raise ValueError(msg)
     
150.                 # if latitude is slightly greater than 90, truncate to 90
     
151.                 lats = [max(min(lat, 90.0), -90.0) for lat in lats]
     
152.                 x, y = lons, lats
     
153.                 coords.append(list(zip(x, y)))
     
154.                 attdict = {}
     
155.                 for r, key in zip(rec, fields[1:]):
     
156.                     attdict[key[0]] = r
     
157.                 # add information about ring number to dictionary.
     
158.                 attdict['RINGNUM'] = ringnum
     
159.                 attdict['SHAPENUM'] = npoly
     
160.                 attributes.append(attdict)
     
161.     # draw shape boundaries for polylines, polygons  using LineCollection.
     
162.     if shptype not in [1, 8] and drawbounds:
     
163.         # get current axes instance (if none specified).
     
164.         import matplotlib.pyplot as plt
     
165.         ax = ax or plt.gca()
     
166.         # make LineCollections for each polygon.
     
167.         lines = LineCollection(coords, antialiaseds=(antialiased,))
     
168.         lines.set_color(color)
     
169.         lines.set_linewidth(linewidth)
     
170.         lines.set_label('_nolabel_')
     
171.         if zorder is not None:
     
172.             lines.set_zorder(zorder)
     
173.         ax.add_collection(lines)
     
174.     # set axes limits to fit map region.
     
175.     # self.set_axes_limits(ax=ax)
     
176.     # # clip boundaries to map limbs
     
177.     # lines,c = self._cliplimb(ax,lines)
     
178.     # info = info + (lines,)
     
179.     info[name] = coords
     
180.     info[name + '_info'] = attributes
     
181.     return info
     
182. 
     
183. 
     
184. if __name__ == '__main__':
     
185.     file = 'test.shp'
     
186.     readshapefile(file.replace('.shp', ''), os.path.basename(file), color='k', linewidth=1)
     

复制代码

平流层的萝卜

无投影和有投影，可否上个图或者举个例子

非对称

平流层的萝卜 发表于 2018-2-13 09:59
无投影和有投影，可否上个图或者举个例子

这段代码 适合在自定义坐标下绘制shp文件，取决于你的shp文件的坐标系统和自定义坐标系统匹配

尽头的尽头

感谢分享~~快过年了还没闲着，新年快乐~~

weinihou

大神出来分享东西了啊~

顶一个

13top

新年快乐~~

iam428jll

顶一个！

空空大蕃薯

顶一个

kmdqlb

老师，你好，请问这个图中 天气现象 图例是怎么画出来的?谢谢