Exercise 3 - Exploring map projections

Exercise 3 - Exploring map projections#

Aim: To create a map of data from reanalysis heat fluxes.

Data: You will need to download files from ICDC. We will be using NCEP Reanalysis heat fluxes at the ocean surface, and you will need one snapshot of sensible heat flux, latex heat flux, net shortwave and net longwave radiation flux.

Directions: Create an *.ipynb and 4 figures: once of heat flux components using matplotlib and one with the net heat flux, and the next two figures using pyGMT to plot net heat flux with two different projections.

Create a notebook#

Create an *.ipynb containing the commands for this assignment, or copy this file.

File naming convention

Name your python notebook something useful ex<X>-<Lastname>-<slug>-seaocn.ipynb where you replace <X> with the exercise number and <slug> with the short slug to name the topic, and <Lastname> with your last name.

Figures should be named something like ex<X>fig<Y>-<Lastname>-<slug>-seaocn.png where you replace <X> with the exercise number, <Y> with the figure number, and <Lastname> with your last name.
Import necessary packages.

For example, matplotlib and pandas and numpy and xarray. You may also need
```
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import xarray as xr
from datetime import datetime
```
If you are missing any of these packages, please refer to Resources: Python.

# Your code here
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import xarray as xr
from datetime import datetime

# Some extra colormaps
import cmocean

# Fancy but non-intuitive
import pygmt

Downloading data#

Download some data. First decide what year you’d like to work with; anything from 1948 to 2023 is fine. Then navigate to each of the 4 components of surface heat fluxes on the ICDC page https://www.cen.uni-hamburg.de/en/icdc/data/atmosphere/reanalysis-atmosphere/ncep.html. And choose “Access NCEP data via OPeNDAP”. Download one of each heat flux file for your chosen year. Note that these files are about 30 mb each.

Tip

Double check when you commit that these files are not getting committed to the gitlab repository.
Load one file. Make a basic exploration. How big are the data? What are the coordinates?

# Your code here

Using loops#

Concatenate to build file names To make the code easier to change, we’ll recommend building the filename out of components. This gives you the possibility later of downloading a different year and only updating part of the code.
```
myyear = 2001
filepath = 'data/' # If you have a file folder containing data
fpre = ['lhtfl', 'nswrs']  # Add to this all the names of the variables
fpost = '.sfc.gauss.' + str(myyear) + '.nc'
```
This will give us a filename once we contenate together the various pieces:
```
fname = filepath + fpre[0] + fpost
```
Learn about loops. Now you could load these one after another on separate lines. Or you can use a loop.

In python, one of the basic loop types is a for-do loop.

See also

About for-do loops: https://www.w3schools.com/python/python_for_loops.asp

Here, we’ll use a loop that goes through a simple vector which has the length the same as the length of fpre. First check what is the length of fpre?
```
print(len(fpre))
```
It should be 4. If it’s not, then go back and edit where you define the list fpre. Verify that it is indeed of type list by doing
```
print(type(fpre))
```
Write a loop to build filenames Now construct your for-do loop which builds the different variable names

Check that it is correctly producing the name of a file.
Add a line to the loop to load the data. If you had a single filename correctly named in the string fname, then you could use xr.open_dataset(fname). For example
```
mydata = xr.open_dataset(fname)
```
This line needs to be inside the loop (i.e., indented to the same indentation as the line fname = filepath....

Note

After doing this, you’ve only actually got one data variable in mydata. Do a print(mydata) within the loop to see how its contents change as you go through the loop.

# Your code here

Loading multiple files: Using a `dictionary`#

Python has a variable type called a “dictionary” which is used to store “key - value” pairs.

Merge data into a single `xarray` dataset#

The dictionary of xarray datasets was kind of useful, but with xarray we don’t need to bother using a dictionary to store the data. Instead, we can use the command xr.merge to combine the similar datatypes (same coordinates, same dimensions, but different variables: latent, sensible heat flux, and shortwave and longwave radiation).

# Merge the datasets using xr.merged
all_flux = xr.merge([
    flux_components['lhtfl'],
    flux_components['shtfl'],
    flux_components['nswrs'],
    flux_components['nlwrs']])
print(all_flux)

# Repeat this, but make an average over Jan - March
ann_flux = all_flux.mean(dim='time', keep_attrs=True)
print(ann_flux)

#winter_flux = 

Making a seasonal average#

Now we’re going to use some of the fancier features of the xarray data construction. We’d like to make an average over 1 January through 31 March for your chosen year. Since we’ve stored the data all in a single xarray dataset, we can calculate the mean with one line of code.

An annual average would be computed as:
```
ann_flux = all_flux.mean(dim='time', keep_attrs=True)
```
What happens if you don’t include the keep_attrs=True option? Try deleting it and see what changes.
How can you tell Python to only average over a certain time range, or specified months?

Fig. 1. Plot with `matplotlib`#

Now we’d like to take a look at the data for a single snapshot (a single time). The example code below will choose the very first frame (where the time index is 0), and plot the latent heat flux. Update the code in order to plot four fields (sensible, latent, shortwave and longwave).

# Plot the fields
# choose the index of the snapshot to show
itime = 0
map1 = winter_flux.lhtfl[itime,:,:]

fig, axs = plt.subplots(2,2)
axs[0,0].contourf(data1.lon, data1.lat, map1, cmap='RdYlBu')
axs[0,0].set_title('Latent heat flux')
axs[0,0].set_ylabel('Latitude')

# Cumbersome date time to string
d = data1.time[itime].dt.strftime('%Y').values
fig.suptitle('NCEP Reanalysis \n' + d + 'winter')

fig.savefig('fig1-Lastname-heatflux.png')

Fig. 2 - Calculate and plot net wintertime heat flux#

Combine your four heat flux term into a single net netflux. This is the total experienced by the ocean at a single lat/lon location.

Warning

Not all signs are the same! Check the individual fields to see what the sign of the fluxes should be.

# Your code here
#net_winter_flux = 

## Fig. 2. Plot with `matplotlib`

#fig.savefig('figf-Lastname-heatflux.png')

Fig 3 & 4. Using `pyGMT`#

PyGMT is especially good for geophysical quantities (and geophysics). However, the formatting language looks a little strange if you’re used to matplotlib.

Update the code below#

The code below should plot fields of heat flux from NCEP.

Update the code to instead plot the single net heat flux. We don’t want to bother with making 4 figures every time, so just use your average from fig 2.
Experiment with a different font color or size for the primary font annotations (FONT_ANNOT_PRIMARY).

See also

Various defaults that can be set: https://docs.generic-mapping-tools.org/latest/gmt.conf.html
Try different figure sizes. Try varying figsize and margins.
Important Can you change the projection? Note that the formatting of projection strings is strange.
For the figure 3 and figure 4, you will want to save these to the shared_figures/ folder for a discussion on map projections.
1. Choose a subregion of the globe. Something spanning a latitude range of no more than 30 degrees, and longitude span of your choice. Plot the data from this region, using a projection of your choice.
2. Choose either (i) a different global projection or (i) a polar projection. Plot the data from this region using a projection of your choice.

# Some sample code.
fig = pygmt.Figure()

with pygmt.config(FONT_ANNOT_PRIMARY="20p,Helvetica,blue", FONT_LABEL="15p,Helvetica,red"):
    with fig.subplot(nrows=2, ncols=2,  figsize=("30c", "22c"), sharex="b", sharey="l", margins="1c"):
        for i in range(len(fpre)):
            fname = filepath + fpre[i] + fpost
            grid = fname + '?' + fpre[i]

            
            with fig.set_panel(panel=i):  # sets the current panel
                fig.basemap(
                    region="g",
                    projection="Cyl_stere/150/-20/?",
                    frame=['WSne+t'+fpre[i], "xa90", "ya30"],
                )
                fig.grdimage(
                    grid=grid,
                    cmap='no_green',
                )
                fig.coast(shorelines="1/0.5p,black")
                #fig.colorbar(frame=['x+l' +fpre[i], "y+lW/m@+2@+"])

fig.show()
fig.savefig('fig4-lastname-pygmt.png')

# Figure 3

# Figure 4

Exercise 3 - Exploring map projections

Contents

Exercise 3 - Exploring map projections#

Create a notebook#

Downloading data#

Using loops#

Loading multiple files: Using a dictionary#

Merge data into a single xarray dataset#