1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
import netCDF4 as nc
import numpy as np
import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
def basemask(cs, ax, map, shpfile):
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
if shape_rec.record[3] >= 0:
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
vertices.append(map(pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform = ax.transData)
for contour in cs.collections:
contour.set_clip_path(clip)
shpfile = u'E:\MATLAB\shp\中国行政区_包含南海九段线'
data = nc.Dataset('E:\python\scripts\sss\pres.mon.ltm.nc')
pres = data.variables['pres'][0,14:35,28:57]
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
lons, lone = lon[28], lon[57]
lats, late = lat[35], lat[14]
nx = pres.shape[1]
ny = pres.shape[0]
fig = plt.figure()
ax = fig.add_subplot(111)
map = Basemap(llcrnrlon = lon1, llcrnrlat = lat1, urcrnrlon = lon2, urcrnrlat = lat2,
projection = 'cyl', ax = ax)
map.drawparallels(np.arange(lats, late + 1, 5), labels = [1, 0, 0, 0])
map.drawmeridians(np.arange(lons, lone + 1, 10), labels = [0, 0, 0, 1])
x, y = map.makegrid(nx, ny)
# 转换坐标
x, y = map(x, y[::-1])
map.readshapefile(shpfile, 'China')
cs = map.contourf(x, y, pres, vmin = int(pres.min()), vmax = int(pres.max()))
basemask(cs, ax, map, r'E:\MATLAB\shp\bou2_4p')
axins = zoomed_inset_axes(ax, 0.8, loc = 3)
axins.set_xlim(lons + 38, lons + 52)
axins.set_ylim(lons, lons + 20)
map2 = Basemap(llcrnrlon = lons + 38, llcrnrlat = lats, urcrnrlon = lons + 52, urcrnrlat = lats + 20, ax = axins)
map2.readshapefile(shpfile, 'China')
cs2 = map2.contourf(x, y, pres, vmin = int(pres.min()), vmax = int(pres.max()))
basemask(cs2, axins, map2, r'E:\MATLAB\shp\bou2_4p')
mark_inset(ax, axins, loc1=2, loc2=4, fc = "none", ec = "none")
plt.show()
|
