|
| 1 | +""" |
| 2 | +Calculating and plotting the divergence of sea currents |
| 3 | +======================================================= |
| 4 | +
|
| 5 | +In this recipe, we will calculate the divergence of depth-averaged |
| 6 | +currents in the Irish Sea, then plot the divergence as a contour |
| 7 | +fill plot underneath the vectors themselves in the form of a vector plot. |
| 8 | +""" |
| 9 | + |
| 10 | +# %% |
| 11 | +# 1. Import cf-python and cf-plot: |
| 12 | +import cfplot as cfp |
| 13 | +import cf |
| 14 | + |
| 15 | +# %% |
| 16 | +# 2. Read the fields in. This dataset consists of depth-averaged eastward and |
| 17 | +# northward current components plus the sea surface height above sea level and |
| 18 | +# is a gridded dataset, with grid resolution of 1.85 km, covering the entire |
| 19 | +# Irish Sea area. It was found via the CEDA Archive at the location of: |
| 20 | +# https://catalogue.ceda.ac.uk/uuid/1b89e025eedd49e8976ee0721ec6e9b5, with |
| 21 | +# DOI of https://dx.doi.org/10.5285/031e7ca1-9710-280d-e063-6c86abc014a0: |
| 22 | +f = cf.read("~/recipes/POLCOMS_WAM_ZUV_01_16012006.nc") |
| 23 | + |
| 24 | +# %% |
| 25 | +# 3. Get the separate vector components, which are stored as separate fields. |
| 26 | +# The first, 'u', corresponds to the eastward component and the second, 'v', |
| 27 | +# the northward component: |
| 28 | +u = f[0] |
| 29 | +v = f[1] |
| 30 | + |
| 31 | +# %% |
| 32 | +# 4. Squeeze the fields to remove the size 1 axes in each case: |
| 33 | +u = u.squeeze() |
| 34 | +v = v.squeeze() |
| 35 | + |
| 36 | +# %% |
| 37 | +# 5. Consider the currents at a set point in time. To do this we |
| 38 | +# select one of the 720 datetime sample points in the fields to |
| 39 | +# investigate, in this case by subspacing to pick out a particular |
| 40 | +# datetime value we saw within the time coordinate data of the field (but |
| 41 | +# you could also use indexing or filtering to select a specific value). |
| 42 | +# Once we subspace to one datetime, we squeeze out the size 1 time axis |
| 43 | +# in each case: |
| 44 | +chosen_time = "2006-01-15 23:30:00" # 720 choices to pick from, try this one! |
| 45 | +u_1 = u.subspace(T=cf.dt(chosen_time)) |
| 46 | +v_1 = v.subspace(T=cf.dt(chosen_time)) |
| 47 | +u_1 = u_1.squeeze() |
| 48 | +v_1 = v_1.squeeze() |
| 49 | + |
| 50 | +# %% |
| 51 | +# 6. |
| 52 | +# When inspecting the u and v fields using cf inspection methods such as |
| 53 | +# from print(u_1.data) and u_1.data.dump(), for example, we can see that there are |
| 54 | +# lots of -9999 values in their data array, apparently used as a |
| 55 | +# fill/placeholder value, including to indicate undefined data over the land. |
| 56 | +# In order for these to not skew the data and dominate the plot, we need |
| 57 | +# to mask values matching this, so that only meaningful values remain. |
| 58 | +u_2 = u_1.where(cf.lt(-9000), cf.masked) |
| 59 | +v_2 = v_1.where(cf.lt(-9000), cf.masked) |
| 60 | + |
| 61 | +# %% |
| 62 | +# 7. Calculate the divergence using the 'div_xy' function operating on the |
| 63 | +# vector eastward and northward components as the first and second argument |
| 64 | +# respectively. We need to calculate this for the latitude-longitude plane |
| 65 | +# of the Earth, defined in spherical polar coordinates, so we must specify |
| 66 | +# the Earth's radius for the appropriate calculation: |
| 67 | +div = cf.div_xy(u_2, v_2, radius="earth") |
| 68 | + |
| 69 | +# %% |
| 70 | +# 8. Generate the final plot. First we configure the overall plot by |
| 71 | +# making the map higher resolution, to show the coastlines of the UK and |
| 72 | +# Ireland in greater detail, and changing the colourmap to better reflect |
| 73 | +# the data which can be positive or negative, i.e. has 0 as the 'middle' |
| 74 | +# value of significance, so should use a diverging colour map. Then we |
| 75 | +# plot the current vectors, noting we had to play around with the 'stride' |
| 76 | +# and 'scale' parameter values to adjust the vector spacing and size so that |
| 77 | +# the vector field is best represented and visible without over-cluttering |
| 78 | +# the plot. Finally we plot the divergence as a contour plot without any |
| 79 | +# lines showing. This compound plot is saved on one canvas using 'gopen' |
| 80 | +# and 'gclose' to wrap the two plotting calls: |
| 81 | +cfp.mapset(resolution="10m") |
| 82 | +cfp.cscale("ncl_default") |
| 83 | +cfp.gopen( |
| 84 | + file=f"irish-sea-currents-divergence-{chosen_time.replace(' ', '-')}.png") |
| 85 | +cfp.vect(u=u_2, v=v_2, stride=6, scale=3, key_length=1) |
| 86 | +cfp.con( |
| 87 | + div, |
| 88 | + lines=False, |
| 89 | + title=( |
| 90 | + f"Depth-averaged Irish Sea currents at {chosen_time} with " |
| 91 | + "their divergence" |
| 92 | + ) |
| 93 | +) |
| 94 | +cfp.gclose() |
0 commit comments