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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
| import salem
from salem import geogrid_simulator
from salem import wgs84
fpath = r'./namelist.wps'
g, maps,rect= geogrid_simulator(fpath)
x0,y0 = pyproj.transform(g[0].proj,wgs84,g[0].extent[0], g[0].extent[2])
x1,y1 = pyproj.transform(g[0].proj,wgs84,g[0].extent[1], g[0].extent[3])
def background(ax):
fname = './bou2_4p.shp'
shape_feature = ShapelyFeature(Reader(fname).geometries(),
ccrs.PlateCarree(),facecolor='none',edgecolor ='k',linewidth = 0.25)
ax.add_feature(shape_feature)
fname = './中国轮廓_含南海.shp'
shape_feature = ShapelyFeature(Reader(fname).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor ='k',linewidth = 0.75)
ax.add_feature(shape_feature)
fname = './country1.shp'
shape_feature = ShapelyFeature(Reader(fname).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor ='k',linewidth = 0.75)
ax.add_feature(shape_feature)
ax.set_extent([x0,125,y0+3,55], crs=ccrs.PlateCarree())
xticks = [80,90,100,110,120,130,140]
yticks = [10,20,30,40,50,60]
fig.canvas.draw()
ax.gridlines(xlocs=xticks, ylocs=yticks, color='gray', alpha=0.5, linestyle='--',linewidth = 0.75)
ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lambert_xticks(ax, xticks)
lambert_yticks(ax, yticks)
def find_side(ls, side):
"""
Given a shapely LineString which is assumed to be rectangular, return the
line corresponding to a given side of the rectangle.
"""
minx, miny, maxx, maxy = ls.bounds
points = {'left': [(minx, miny), (minx, maxy)],
'right': [(maxx, miny), (maxx, maxy)],
'bottom': [(minx, miny), (maxx, miny)],
'top': [(minx, maxy), (maxx, maxy)],}
return sgeom.LineString(points[side])
def lambert_xticks(ax, ticks):
"""Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
te = lambda xy: xy[0]
lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T
xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te)
ax.xaxis.tick_bottom()
ax.set_xticks(xticks)
ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels])
def lambert_yticks(ax, ticks):
"""Draw ricks on the left y-axis of a Lamber Conformal projection."""
te = lambda xy: xy[1]
lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
yticks, yticklabels = _lambert_ticks(ax, ticks, 'left', lc, te)
ax.yaxis.tick_left()
ax.set_yticks(yticks)
ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels])
def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):
"""Get the tick locations and labels for an axis of a Lambert Conformal projection."""
outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist())
axis = find_side(outline_patch, tick_location)
n_steps = 30
extent = ax.get_extent(ccrs.PlateCarree())
_ticks = []
for t in ticks:
xy = line_constructor(t, n_steps, extent)
proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1])
xyt = proj_xyz[..., :2]
ls = sgeom.LineString(xyt.tolist())
locs = axis.intersection(ls)
if not locs:
tick = [None]
else:
tick = tick_extractor(locs.xy)
_ticks.append(tick[0])
ticklabels = copy(ticks)
while True:
try:
index = _ticks.index(None)
except ValueError:
break
_ticks.pop(index)
ticklabels.pop(index)
return _ticks, ticklabels
proj = salem.proj_to_cartopy(g[0].proj)
fig = plt.figure(figsize=(6.,5), frameon=True)
ax1 = plt.subplot(111, projection=proj)
cs = cn_back_plot(ax1)
ax1.background_img(name='ne_shaded3',resolution='low')
ax1.add_geometries([rect[0][0]], crs = proj,facecolor='none', edgecolor='k', lw=2.5)
plt.savefig("./模拟嵌套网格.png",dpi=400)
plt.show()
|
Kommentare