In [1]:
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import netCDF4 as nc
from mpl_toolkits.basemap import Basemap
import cartopy.crs as ccrs
In [2]:
# create a basemap instance
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution='c',lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
m.drawcoastlines() # draw coastlines.
m.drawmapboundary(fill_color='aqua') # draw a boundary around the map, fill the background.
m.fillcontinents(color='coral',lake_color='aqua') # fill continents, set lake color same as ocean color.
plt.show()
In [3]:
def basemapping(subplot, i):
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution='c',lat_1=45.,lat_2=55,lat_0=50,lon_0=-107., ax = subplot)
if i == 0:
m.drawlsmask(land_color='coral',ocean_color='aqua',lakes=True)
if i == 1:
m.bluemarble()
if i == 2:
m.shadedrelief()
if i == 3:
m.etopo()
return m
In [4]:
fig, axes = plt.subplots(2,2, figsize=(15, 15), facecolor='w', edgecolor='k') # establish subplot figure
fig.subplots_adjust(wspace = 0.1, hspace = 0.001)
axes[0,0].set_title('Land-Sea Mask', fontsize = 16)
axes[0,1].set_title('NASA ‘Blue Marble’', fontsize = 16)
axes[1,0].set_title('Shaded Relief', fontsize = 16)
axes[1,1].set_title('Etopo', fontsize = 16)
for i, ax in enumerate(fig.axes):
axes = basemapping(ax, i)
In [5]:
fig, axes = plt.subplots(1, figsize=(30, 15))
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution='c',lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
m.shadedrelief()
# draw parallels and meridians.
# label parallels on right and top
# meridians on bottom and left
parallels = np.arange(0.,81,10.)
meridians = np.arange(10.,351.,20.)
#labels = [left,right,top,bottom]
m.drawparallels(parallels,labels=[False,True,True,False])
m.drawmeridians(meridians,labels=[True,False,False,True])
# plot blue dot on Miami, FL and label it as such.
lon, lat = -80.19, 25.76 # Location of MIA
# convert to map projection coords.
# Note that lon,lat can be scalars, lists or numpy arrays.
xpt,ypt = m(lon,lat)
# convert back to lat/lon
lonpt, latpt = m(xpt,ypt,inverse=True)
m.plot(xpt,ypt,'rx', markersize = 12, markeredgewidth = 5) # plot a red x there
# put some text next to the +, offset a little bit
# (the offset is in map projection coordinates)
plt.text(xpt-1000000,ypt+300000,'Miami (%5.1fW,%3.1fN)' % (lonpt,latpt), fontsize = 16, bbox=dict(facecolor='red', alpha=0.5))
plt.show()
In [6]:
concat_t = xr.open_dataset('temps_9520.nc')
concat_t
Out[6]:
In [7]:
concat_h = xr.open_dataset('heights_9520.nc')
concat_h
Out[7]:
In [8]:
%%time
data = concat_h.update(concat_t)
data
In [9]:
data['air'] = data['air'] - 273.15 #convert to Celsius
data
Out[9]:
In [10]:
seasonal = data.groupby('time.season').mean(dim = 'time') # use the 'powerful' groupby funtion to find seasonal averages
seasonal
Out[10]:
In [11]:
def mapping(subplot):
ll = 180. # lower lon bounds
ul = 330. # upper lon bounds
lllat=20. # lower lat bounds
urlat=80. # upper lat bounds
m = Basemap(projection='cyl',llcrnrlat=lllat,urcrnrlat=urlat,\
llcrnrlon=ll,urcrnrlon=ul,resolution='l' , ax = subplot, lon_0 = 0) # establishing basemap, ax = subplot = key feature to make function/subplot loop work
m.drawcountries(linewidth=0.2) # map design specifications
m.drawcoastlines(linewidth=0.2)
m.drawmapboundary(fill_color='white')
parallels = np.arange(25.,85.,20.) # setting range of where and how many lat/lon grid lines you want
meridians = np.arange(-180.,181.,30.)
m.drawparallels(parallels, labels=[1,0,0,0], fontsize = 20, dashes=[6,900]) # actually drawing them to map
m.drawmeridians(meridians, labels=[1,1,0,1], fontsize=20, dashes=[6,900])
return m
In [12]:
def contours(subplot, temps, heights):
lat_lon1 = False # set to true if crossing meridian or your data might have gaps or not show up at all!
# Criteria for lat_lon in contour or contourf
# extract raw numpy arrays from xarray
lons = temps.lon
lats = temps.lat
temps = temps.values
hgt = heights.values
# call the mapping function to get your precalled mapping instance
m = mapping(subplot)
px,py = np.meshgrid(lons, lats) # mesh lat/lon numpy arrays
# compute native map projection coordinates of lat/lon grid
# does not actually know those are lat/lons!!
# this is a function of basemap and how it converts to given projection from your lat/lon array
x,y = m(px, py) # apply meshed grid to map projection
levels1 = np.arange(-40, 42, 2) # set interval for temperature contours
# contourf = contour fill for our temperatures
# note still calling to mapping instance
tc = m.contourf(x, y, temps, 35, extend='both', cmap = 'seismic', latlon = lat_lon1, levels = levels1, vmin = -40, vmax = 40)
# contour = line contours
levels2 = np.arange(5000, 6050, 100) # set interval for geopotential contours
h_contour = m.contour(x, y, hgt, 15, extend = 'both', latlon = lat_lon1, colors = 'k', alpha = 0.75, levels = levels2) # alpha sets transparency
h_contour.clabel(fontsize=15, colors='k', inline=1, inline_spacing=8, levels = levels2, fmt='%i', rightside_up=True, use_clabeltext=True);
return tc # return contour plot for use in colorbar
In [13]:
fig, axes = plt.subplots(2,2, figsize=(30, 15), facecolor='w', edgecolor='k', sharex =False, sharey = False) # establish subplot figure
labels = seasonal.season.values # set array for plot labels
for i, ax in enumerate(fig.axes): # enumerate axes so that you dont have [0,1] etc., just 0,1,2,3
ax.set_title(labels[i], fontsize = 35)
ax1 = contours(ax, seasonal['air'][i], seasonal['hgt'][i]) # call created function and loop through and plot seasons based on data
cax = fig.add_axes([0.93, 0.17, 0.01, 0.67]) # customizing colorbar for one large one for plot by creating new axis
cb = plt.colorbar(ax1, cax=cax)
cb.ax.tick_params(labelsize=25)
fig.suptitle('Seasonally Averaged Geopotential Heights (m) and Temperatures ($^\circ$C), 1995-2000', fontsize=40, ha = 'center') # don't forget to label :D
fig.text(1, 0.6, "Temperature ($^\circ$C)", ha='center', fontsize=30, rotation = 270);
#plt.savefig('nameyourplothere.png', bbox_inches = "tight") # bbox_inches limits the extra whitespcae once saved
In [14]:
ax = plt.axes(projection=ccrs.PlateCarree(), )
ax.coastlines(resolution='110m', color = 'green')
# Save the plot by calling plt.savefig() BEFORE plt.show()
# plt.savefig('coastlines.pdf')
# plt.savefig('coastlines.png')
plt.show()
In [15]:
fig, axes = plt.subplots(1, figsize=(15, 15))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
ny_lon, ny_lat = -75, 43
delhi_lon, delhi_lat = 77.23, 28.61
plt.plot([ny_lon, delhi_lon], [ny_lat, delhi_lat],
color='blue', linewidth=2, marker='o',
transform=ccrs.Geodetic(),
)
plt.plot([ny_lon, delhi_lon], [ny_lat, delhi_lat],
color='red', linestyle='--',
transform=ccrs.PlateCarree(),
)
plt.text(ny_lon - 3, ny_lat - 12, 'New York',
horizontalalignment='right',
transform=ccrs.Geodetic())
plt.text(delhi_lon + 3, delhi_lat - 12, 'Delhi',
horizontalalignment='left',
transform=ccrs.Geodetic())
plt.show()
In [16]:
lon = np.linspace(-80, 80, 25)
lat = np.linspace(0, 70, 25)
lon2d, lat2d = np.meshgrid(lon, lat)
data = np.cos(np.deg2rad(lat2d) * 4) + np.sin(np.deg2rad(lon2d) * 4)**2
In [17]:
# The projection keyword determines how the plot will look
plt.figure(figsize=(16, 13))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
ax.contourf(lon, lat, data) # didn't use transform, but looks ok...
plt.show()
In [18]:
# The data are defined in lat/lon coordinate system, so PlateCarree()
# is the appropriate choice:
data_crs = ccrs.PlateCarree()
# The projection keyword determines how the plot will look
plt.figure(figsize=(16, 13))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
ax.contourf(lon, lat, data, transform=data_crs)
plt.show()
In [19]:
# Now we plot a rotated pole projection
projection = ccrs.RotatedPole(pole_longitude=-177.5, pole_latitude=37.5)
plt.figure(figsize=(16, 13))
ax = plt.axes(projection=projection)
ax.set_global()
ax.coastlines()
ax.contourf(lon, lat, data) # didn't use transform, uh oh!
plt.show()
In [20]:
# A rotated pole projection again...
projection = ccrs.RotatedPole(pole_longitude=-177.5, pole_latitude=37.5)
plt.figure(figsize=(16, 13))
ax = plt.axes(projection=projection)
ax.set_global()
ax.coastlines()
# ...but now using the transform argument
ax.contourf(lon, lat, data, transform=data_crs)
plt.show()
In [21]:
# We can choose any projection we like...
projection = ccrs.InterruptedGoodeHomolosine()
plt.figure(figsize=(16, 13))
ax1 = plt.subplot(211, projection=projection)
ax1.set_global()
ax1.coastlines()
ax2 = plt.subplot(212, projection=ccrs.NorthPolarStereo())
ax2.set_extent([-180, 180, 20, 90], crs=ccrs.PlateCarree())
ax2.coastlines()
# ...as long as we provide the correct transform, the plot will be correct
ax1.contourf(lon, lat, data, transform=data_crs)
ax2.contourf(lon, lat, data, transform=data_crs)
plt.show()
In [22]:
import cartopy.crs as ccrs
from cartopy.mpl.geoaxes import GeoAxes
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from mpl_toolkits.axes_grid1 import AxesGrid
In [23]:
seasonal
Out[23]:
In [24]:
fig = plt.figure(figsize=(30, 15))
labels = seasonal.season.values
projection = ccrs.PlateCarree()
axes_class = (GeoAxes,
dict(map_projection=projection))
axgr = AxesGrid(fig, 111, axes_class=axes_class, # setting the axes class with previously defined GeoAxes
nrows_ncols=(2, 2),
axes_pad=1,
cbar_location='right',
cbar_mode='single',
cbar_pad=0.2,
cbar_size='3%',
label_mode='') # note the empty label_mode
for i, ax in enumerate(axgr):
ax.set_extent([-180, -30, 20, 80], crs = ccrs.PlateCarree())
ax.set_title(labels[i], fontsize = 35)
ax.coastlines()
ax.set_xticks(np.arange(-180, -30, 30), crs=projection)
ax.set_yticks(np.arange(25, 80, 20), crs=projection)
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.tick_params(labelsize=25)
levels1 = np.arange(-40, 42, 2)
levels2 = np.arange(5000, 6050, 100)
t = ax.contourf(seasonal.lon, seasonal.lat, seasonal['air'][i],
transform=projection,
cmap='seismic', levels = levels1)
h_contour = ax.contour(seasonal.lon, seasonal.lat, seasonal['hgt'][i], 15, extend = 'both', colors = 'k', alpha = 0.75, levels = levels2) # alpha sets transparency
h_contour.clabel(fontsize=15, colors='k', inline=1, inline_spacing=8, levels = levels2, fmt='%i', rightside_up=True, use_clabeltext=True);
cb = axgr.cbar_axes[0].colorbar(t)
cb.ax.tick_params(labelsize=25)
fig.suptitle('Seasonally Averaged Geopotential Heights (m) and Temperatures ($^\circ$C), 1995-2000', fontsize=40, ha = 'center') # don't forget to label :D
fig.text(0.95, 0.6, "Temperature ($^\circ$C)", ha='center', fontsize=30, rotation = 270);
#plt.savefig('nameyourplothere.png', bbox_inches = "tight") # bbox_inches limits the extra whitespcae once saved
plt.show()
In [ ]: