Flesh out literate prog. comments for recipe 20

Author: sadielbartholomew

docs/source/recipes/plot_20_recipe.py

@@ -1,82 +1,90 @@
 """
-Calculating and plotting the divergence of sea current vectors
-==============================================================
+Calculating and plotting the divergence of sea currents
+=======================================================
 
 In this recipe, we will calculate the divergence of depth-averaged
 currents in the Irish Sea, then plot the divergence as a contour
-fill plot underneath the vectors themselves in a vector plot.
+fill plot underneath the vectors themselves in the form of a vector plot.
 """
 
+# %%
+# 1. Import cf-python and cf-plot:
 import cfplot as cfp
 import cf
 
-f = cf.read("~/recipes_break/new/POLCOMS_WAM_ZUV_01_16012006.nc")
-print(f)
+# %%
+# 2. Read the fields in:
+f = cf.read(
+    "~/summerstudents/final-recipes/new-required-datasets/POLCOMS_WAM_ZUV_01_16012006.nc")
 
-# Get separate vector components
+# %%
+# 3. Get the separate vector components, which are stored as separate fields.
+# The first, 'u', corresponds to the eastward component and the second, 'v',
+# the northward component:
 u = f[0]
 v = f[1]
-print(u)
-print(v)
 
-# First get rid of the ocean sigma coord size 1 axis
+# %%
+# 4. Squeeze the fields to remove the size 1 axes in each case:
 u = u.squeeze()
 v = v.squeeze()
 
-# Now we need to use some means to condense the u and v fields in the same way into
-# having 1 time point, not 720 - for example we can just pick a time value out:
-print("Times are", v.construct("T").data.datetime_as_string)
-chosen_time = "2006-01-15 23:30:00"  # 720 choices to choose from!
-v_1 = v.subspace(T=cf.dt(chosen_time))
+# %%
+# 5. Consider the currents at a set point in time. To do this we
+# select one of the 720 datetime sample points in the fields to
+# investigate, in this case by subspacing to pick out a particular
+# datetime value we saw within the time coordinate data of the field (but
+# you could also use indexing or filtering to select a specific value).
+# Once we subspace to one datetime, we squeeze out the size 1 time axis
+# in each case:
+chosen_time = "2006-01-15 23:30:00"  # 720 choices to pick from, try this one!
 u_1 = u.subspace(T=cf.dt(chosen_time))
-v_1 = v_1.squeeze()
+v_1 = v.subspace(T=cf.dt(chosen_time))
 u_1 = u_1.squeeze()
-print(u_1)
-print(v_1)
-
-# Are now in a plottable form! Let's give it a go:
-### cfp.vect(u=u_1, v=v_1)
-# Need to play around with the relative length and spacing of the vectors, using these paramters:
-###cfp.vect(u=u_1, v=v_1, key_length=10, scale=50, stride=2)
-
-# Note that there appear to be some really large vectors all pointing in the
-# same direction which are spamming the plot. We need to remove these. By
-# looking at the data we can see what these are and work out how to remove them:
-
-# ... shows more of the array
-
-# Can see there are lots of -9999 values, seemingly used as a fill/placeholder value
-# so we need to remove those so we can plot the menaingful vectors
-# Apply steps to mask the -9999 fill values, which spam the plot, to x_1
-u_2 = u_1.where(cf.lt(-9e+03), cf.masked)
-v_2 = v_1.where(cf.lt(-9e+03), cf.masked)
-print(u_2)
-print(v_2)
-
-# We can even plot the final field, effective wave height, as the
-# background contour!
-w = f[2]
-w_1 = w.subspace(T=cf.dt(chosen_time))
-# This field also needs masking for those data points.
-w_2 = w_1.where(cf.lt(-9e+03), cf.masked)
-print(w_2)
+v_1 = v_1.squeeze()
 
+# %%
+# 6.
+# When inspecting the u and v fields using cf inspection methods such as
+# from print(u_1.data) and u_1.data.dump(), for example, we can see that there are
+# lots of -9999 values in their data array, apparently used as a
+# fill/placeholder value, including to indicate undefined data over the land.
+# In order for these to not skew the data and dominate the plot, we need
+# to mask values matching this, so that only meaningful values remain.
+u_2 = u_1.where(cf.lt(-9000), cf.masked)
+v_2 = v_1.where(cf.lt(-9000), cf.masked)
 
-# Plot divergence in the background
+# %%
+# 7. Calculate the divergence using the 'div_xy' function operating on the
+# vector eastward and northward components as the first and second argument
+# respectively. We need to calculate this for the latitude-longitude plane
+# of the Earth, defined in spherical polar coordinates, so we must specify
+# the Earth's radius for the appropriate calculation:
 div = cf.div_xy(u_2, v_2, radius="earth")
 
-
-# Our final basic plot:
-cfp.mapset(resolution="10m")  # makes UK coastline more high-res
-cfp.gopen(file=f"irish-sea-currents-with-divergence-{chosen_time}.png")
+# %%
+# 8. Generate the final plot. First we configure the overall plot by
+# making the map higher resolution, to show the coastlines of the UK and
+# Ireland in greater detail, and changing the colourmap to better reflect
+# the data which can be positive or negative, i.e. has 0 as the 'middle'
+# value of significance, so should use a diverging colour map. Then we
+# plot the current vectors, noting we had to play around with the 'stride'
+# and 'scale' parameter values to adjust the vector spacing and size so that
+# the vector field is best represented and visible without over-cluttering
+# the plot. Finally we plot the divergence as a contour plot without any
+# lines showing. This compound plot is saved on one canvas using 'gopen'
+# and 'gclose' to wrap the two plotting calls:
+cfp.mapset(resolution="10m")
 cfp.cscale("ncl_default")
+cfp.gopen(
+    file=f"irish-sea-currents-divergence-{chosen_time.replace(' ', '-')}.png")
 cfp.vect(u=u_2, v=v_2, stride=6, scale=3, key_length=1)
 cfp.con(
     div,
     lines=False,
     title=(
-        f"Depth-averaged Irish Sea currents at {chosen_time} "
-        "with their divergence shown"
+        f"Depth-averaged Irish Sea currents at {chosen_time} with "
+        "their divergence"
     )
 )
 cfp.gclose()