Newer
Older
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
#
# Input data
#
# We are going to use a netCDF file for testing the following code.
# The test file islocated in a dedicated online repository:
# https://nextcloud.cism.ucl.ac.be/s/H7XFSi8J4E8dnRK
#
# Download the netCDF prmsl.2000.nc :
# a classic (regular grid) pressure file
#
from netCDF4 import Dataset
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.util import add_cyclic_point
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.path as mpath
my_example_nc_file = 'prmsl.2000.nc'
fh = Dataset(my_example_nc_file, mode='r')
lons = fh.variables['lon'][:]
lats = fh.variables['lat'][:]
times = fh.variables['time'][:]
times_unit = fh.variables['time'].units
var = fh.variables['prmsl']
var_units = fh.variables['prmsl'].units
#------------------------------------------------------------
# Compute a circle in axes coordinates,
# which we can use as a boundary for the map.
# https://scitools.org.uk/cartopy/docs/latest/gallery/lines_and_polygons/always_circular_stereo.html
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
fig = plt.subplots(figsize=(6, 6), dpi=102)
plt.axis('off')
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=0.0, central_latitude=90.0))
xlims = [-180,180]
ylims = [50,90]
ax.set_extent(xlims+ylims, crs=ccrs.PlateCarree())
ax.set_boundary(circle, transform=ax.transAxes)
ax.stock_img()
ax.coastlines("110m", linewidth=0.5, color="black")
# Plot Data
data = np.squeeze(var[1,:,:])
#cs = ax.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree(), shading='auto')
lon = np.arange(0, 360, 2)
data, lon = add_cyclic_point(data, coord=lon)
cs = ax.contourf(lon, lats, data, transform=ccrs.PlateCarree(), transform_first=False)
# Doing the gridlines we want
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.5, color="grey")
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(-180,181,18))
gl.ylocator = mticker.FixedLocator(np.arange(-90,91,10))
# Add Colorbar
cbar = plt.colorbar(cs,orientation="horizontal")
cbar.set_label(var_units)
# Add Title
plt.title('Daily Pressure at Mean Sea Level')
plt.show()