<jupyter_start><jupyter_text>Example Map Plotting<jupyter_code># By line: RRB 2020-07-20
# Script aims to:
# - Load a netCDF file
# - Extract one variable: CO
# - Create contour plot of variable as world map with coastlines
# - Add cyclic point
# - Customize contours and colorbar
# - Add axes labels
# - Add grid lines<jupyter_output><empty_output><jupyter_text>At the start of a Jupyter notebook you need to import all modules that you will use.<jupyter_code>import matplotlib.pyplot as plt
import cartopy.crs as ccrs                 # For plotting maps
import cartopy.feature as cfeature         # For plotting maps
from cartopy.util import add_cyclic_point  # For plotting maps
from pathlib import Path                   # System agnostic paths
import xarray as xr                        # For loading the data arrays
import numpy as np                         # For array creation and calculations
from matplotlib.colors import BoundaryNorm<jupyter_output><empty_output><jupyter_text>Define the directories and file of interest for your results.<jupyter_code>result_dir = Path("/home/buchholz/Documents/code_database/untracked/my-notebook/CAM_Chem_examples")
file = "CAM_chem_merra2_FCSD_1deg_QFED_monthoutput_CO_201801.nc"
file_to_open = result_dir / file
#the netcdf file is now held in an xarray dataset named 'nc_load' and can be referenced later in the notebook
nc_load = xr.open_dataset(file_to_open)
#to see what the netCDF file contains, uncomment below
#nc_load<jupyter_output><empty_output><jupyter_text>Extract the variable of choice at the time and level of choice<jupyter_code>#extract variable
var_sel = nc_load['CO']
#print(var_sel)

#select the surface level at a specific time and convert to ppbv from vmr
#select the surface level for an average over three times and convert to ppbv from vmr
var_srf = var_sel.isel(time=0, lev=55)
var_srf = var_srf*1e09 # 10-9 to ppb
print(var_srf.shape)

#extract grid variables
lat = var_sel.coords['lat']
lon = var_sel.coords['lon']<jupyter_output>(192, 288)<jupyter_text>Add cyclic point to avoid white stripe at lon=0.<jupyter_code>var_srf_cyc, lon_cyc = add_cyclic_point(var_srf, coord=lon)<jupyter_output><empty_output><jupyter_text>Plot the value over the globe.<jupyter_code>plt.figure(figsize=(20,8))

#Define projection
ax = plt.axes(projection=ccrs.PlateCarree())

#define contour levels
cmap = plt.get_cmap('Spectral_r')
clev = np.arange(30, 500, 10)

#plot the data
#plt.contourf(lon_cyc,lat,var_srf_cyc,clev,cmap='Spectral_r',extend='both')
norm = BoundaryNorm(clev, ncolors=cmap.N, clip=True)
plt.pcolormesh(lon_cyc,lat,var_srf_cyc,cmap=cmap, norm=norm)


# add coastlines
ax.add_feature(cfeature.COASTLINE)

#add lat lon grids
gl = ax.gridlines(draw_labels=True, color='grey', alpha=0.8, linestyle='--')
gl.xlabels_top = False
gl.ylabels_right = False

# Titles5
# Main
plt.title("Global map of CAM-chem CO, January 2018",fontsize=18)

# y-axis
ax.text(-0.04, 0.5, 'Latitude', va='bottom', ha='center',
        rotation='vertical', rotation_mode='anchor',
        transform=ax.transAxes)
# x-axis
ax.text(0.5, -0.08, 'Longitude', va='bottom', ha='center',
        rotation='horizontal', rotation_mode='anchor',
        transform=ax.transAxes)
# legend
ax.text(1.15, 0.5, 'CO (ppb)', va='bottom', ha='center',
        rotation='vertical', rotation_mode='anchor',
        transform=ax.transAxes)

plt.colorbar()
plt.show()<jupyter_output><empty_output><jupyter_text>Area weighted regridding<jupyter_code>re=6.37122e06                       # Earth radius (in metres)
rad=4.0 * np.arctan(1.0) / 180.0    # Convert degrees to radians (pi radians per 180 deg)
con  = re * rad                     # constant for determining arc length<jupyter_output><empty_output><jupyter_text>Original grid<jupyter_code>dlon = lon[2].values-lon[1].values
dlat = lat[2].values-lat[1].values
clat = np.cos(lat * rad)            # cosine of latitude
dx = con * dlon * clat              #dx (in metres) at each latitude
dy = con * dlat                     #dy (in metres) is constant
dydx = dy * dx                      #dydx(nlat)
area_wgt = new((/nlat, mlon/), float)
var_srf
     wgt  = new((/nlat, mlon/), float)
     wgt  = conform(wgt, dydx, 0)
      wgt!0 = "lat"
      wgt&lat = newlat
      wgt!1 = "lon"
      wgt&lon = lon
oldlat = xr.DataArray(np.linspace(89.875,-89.875,720).astype('float32'))
oldlat.attrs ['long_name'] = 'latitude'
oldlat.attrs ['units'] = 'degrees_north'
oldlon = xr.DataArray(np.linspace(-179.875,179.875,1440).astype('float32'))
oldlon.attrs ['long_name'] = 'longitude'
oldlon.attrs ['units'] = 'degrees_east'
dlon_old = oldlon[2].values-oldlon[1].values
from math import gcd
gcd(0.25,1.25)
a = np.array([[0,1,2],[2,2,3]])
print(a.shape)
weights = np.array([16,4,2])
print(weights.shape)
np.dot(a,weights)
a = np.repeat(3, 4)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup