From 4ead0943849105d2a2576c2c1c017b959719c822 Mon Sep 17 00:00:00 2001 From: George Pulickan <> Date: Thu, 10 Oct 2024 16:24:09 +0100 Subject: [PATCH 1/9] Add George's original recipe for divergenece of sea current vectors --- docs/source/recipes/plot_20_recipe.py | 72 +++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 docs/source/recipes/plot_20_recipe.py diff --git a/docs/source/recipes/plot_20_recipe.py b/docs/source/recipes/plot_20_recipe.py new file mode 100644 index 0000000000..41c87a19bb --- /dev/null +++ b/docs/source/recipes/plot_20_recipe.py @@ -0,0 +1,72 @@ +import cfplot as cfp +import cf + +f = cf.read("~/recipes_break/new/POLCOMS_WAM_ZUV_01_16012006.nc") +print(f) + +# Get separate vector components +u = f[0] +v = f[1] +print(u) +print(v) + +# First get rid of the ocean sigma coord size 1 axis +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: +chosen_time = "2006-01-16 00:00:00" +v_1 = v.subspace(T=cf.dt(chosen_time)) +u_1 = u.subspace(T=cf.dt(chosen_time)) +v_1 = v_1.squeeze() +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: +print(u.data) +print(u[:10].data.array) + +# ... 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) +print(w_2, w_2[:10].data.array) + +# Our final basic plot: +cfp.mapset(resolution="10m") # makes UK coastline more high-res +cfp.gopen(file="irish-sea-currents.png") +# BTW ignore the warnings below - they aren't relevant. +cfp.vect(u=u_2, v=v_2, stride=2, scale=8, key_length=5) +cfp.levs(min=-5, max=5, step=0.5) +cfp.con(w_1, blockfill=True, lines=False) +cfp.gclose() + +# Ideas for TODOs: +# investigate difference days (do this by changing the 'T=cf.dt("2006-01-16 00:00:00")') datetime +# values to different ones in the time coordinate data so you look at different days, or repace it +# with a collapse over some stat e.g. mean to show the mean over all the times, +# calculate divergence, calculate curl / relative voriticity, calculate absolute voriticity, +# explore the other dataset as well (that covers other dates/times) - you could compare the +# two to effectively compare the currents across different dates. From f75a3cb39e2225a7cecd35740c1c7f72561233de Mon Sep 17 00:00:00 2001 From: "Sadie L. Bartholomew" Date: Thu, 10 Oct 2024 16:33:28 +0100 Subject: [PATCH 2/9] Optimise for readability plotted vector arrow size in recipe 20 --- docs/source/recipes/plot_20_recipe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/recipes/plot_20_recipe.py b/docs/source/recipes/plot_20_recipe.py index 41c87a19bb..84ee663fdd 100644 --- a/docs/source/recipes/plot_20_recipe.py +++ b/docs/source/recipes/plot_20_recipe.py @@ -58,7 +58,7 @@ cfp.mapset(resolution="10m") # makes UK coastline more high-res cfp.gopen(file="irish-sea-currents.png") # BTW ignore the warnings below - they aren't relevant. -cfp.vect(u=u_2, v=v_2, stride=2, scale=8, key_length=5) +cfp.vect(u=u_2, v=v_2, stride=5, scale=2, key_length=1) cfp.levs(min=-5, max=5, step=0.5) cfp.con(w_1, blockfill=True, lines=False) cfp.gclose() From a012e57f7eb58f602bfaa814955d0ab4f3fa9838 Mon Sep 17 00:00:00 2001 From: "Sadie L. Bartholomew" Date: Thu, 10 Oct 2024 18:00:44 +0100 Subject: [PATCH 3/9] Choose better colourmap and time snapshot for recipe 20 --- docs/source/recipes/plot_20_recipe.py | 28 +++++++++++---------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/docs/source/recipes/plot_20_recipe.py b/docs/source/recipes/plot_20_recipe.py index 84ee663fdd..4d7dfb53ae 100644 --- a/docs/source/recipes/plot_20_recipe.py +++ b/docs/source/recipes/plot_20_recipe.py @@ -16,7 +16,8 @@ # 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: -chosen_time = "2006-01-16 00:00:00" +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)) u_1 = u.subspace(T=cf.dt(chosen_time)) v_1 = v_1.squeeze() @@ -32,8 +33,6 @@ # 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: -print(u.data) -print(u[:10].data.array) # ... shows more of the array @@ -52,21 +51,16 @@ # This field also needs masking for those data points. w_2 = w_1.where(cf.lt(-9e+03), cf.masked) print(w_2) -print(w_2, w_2[:10].data.array) + + +# Plot divergence in the background +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="irish-sea-currents.png") -# BTW ignore the warnings below - they aren't relevant. -cfp.vect(u=u_2, v=v_2, stride=5, scale=2, key_length=1) -cfp.levs(min=-5, max=5, step=0.5) -cfp.con(w_1, blockfill=True, lines=False) +cfp.gopen(file=f"irish-sea-currents-with-divergence-{chosen_time}.png") +cfp.cscale("ncl_default") +cfp.vect(u=u_2, v=v_2, stride=6, scale=3, key_length=1) +cfp.con(div, lines=False) cfp.gclose() - -# Ideas for TODOs: -# investigate difference days (do this by changing the 'T=cf.dt("2006-01-16 00:00:00")') datetime -# values to different ones in the time coordinate data so you look at different days, or repace it -# with a collapse over some stat e.g. mean to show the mean over all the times, -# calculate divergence, calculate curl / relative voriticity, calculate absolute voriticity, -# explore the other dataset as well (that covers other dates/times) - you could compare the -# two to effectively compare the currents across different dates. From fd167a89bd8e39da9292f0dfd463a15146e8fc11 Mon Sep 17 00:00:00 2001 From: "Sadie L. Bartholomew" Date: Thu, 10 Oct 2024 18:20:23 +0100 Subject: [PATCH 4/9] Add contour plot title to finalise plot output of recipe 20 --- docs/source/recipes/plot_20_recipe.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/docs/source/recipes/plot_20_recipe.py b/docs/source/recipes/plot_20_recipe.py index 4d7dfb53ae..77558f4eed 100644 --- a/docs/source/recipes/plot_20_recipe.py +++ b/docs/source/recipes/plot_20_recipe.py @@ -62,5 +62,12 @@ cfp.gopen(file=f"irish-sea-currents-with-divergence-{chosen_time}.png") cfp.cscale("ncl_default") cfp.vect(u=u_2, v=v_2, stride=6, scale=3, key_length=1) -cfp.con(div, lines=False) +cfp.con( + div, + lines=False, + title=( + f"Depth-averaged Irish Sea currents at {chosen_time} " + "with their divergence shown" + ) +) cfp.gclose() From e71978a26d2c31aae059fa9d6436e087d6f55a6c Mon Sep 17 00:00:00 2001 From: "Sadie L. Bartholomew" Date: Thu, 10 Oct 2024 23:06:46 +0100 Subject: [PATCH 5/9] Add title and summary for recipe 20 --- docs/source/recipes/plot_20_recipe.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/docs/source/recipes/plot_20_recipe.py b/docs/source/recipes/plot_20_recipe.py index 77558f4eed..2da72726ff 100644 --- a/docs/source/recipes/plot_20_recipe.py +++ b/docs/source/recipes/plot_20_recipe.py @@ -1,3 +1,12 @@ +""" +Calculating and plotting the divergence of sea current vectors +============================================================== + +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. +""" + import cfplot as cfp import cf From 20e5c6e069570ef6f8be7ec3afa0415b8ddba0b2 Mon Sep 17 00:00:00 2001 From: "Sadie L. Bartholomew" Date: Fri, 11 Oct 2024 00:41:18 +0100 Subject: [PATCH 6/9] Flesh out literate prog. comments for recipe 20 --- docs/source/recipes/plot_20_recipe.py | 110 ++++++++++++++------------ 1 file changed, 59 insertions(+), 51 deletions(-) diff --git a/docs/source/recipes/plot_20_recipe.py b/docs/source/recipes/plot_20_recipe.py index 2da72726ff..58c87fc1b7 100644 --- a/docs/source/recipes/plot_20_recipe.py +++ b/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() From 43c421a7a57e3da3200a6be05c80ddf38d4d5725 Mon Sep 17 00:00:00 2001 From: "Sadie L. Bartholomew" Date: Fri, 11 Oct 2024 00:46:59 +0100 Subject: [PATCH 7/9] Add details of data source for recipe 20 --- docs/source/recipes/plot_20_recipe.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/source/recipes/plot_20_recipe.py b/docs/source/recipes/plot_20_recipe.py index 58c87fc1b7..b88798b069 100644 --- a/docs/source/recipes/plot_20_recipe.py +++ b/docs/source/recipes/plot_20_recipe.py @@ -13,7 +13,12 @@ import cf # %% -# 2. Read the fields in: +# 2. Read the fields in. This dataset consists of depth-averaged eastward and +# northward current components plus the sea surface height above sea level and +# is a gridded dataset, with grid resolution of 1.85 km, covering the entire +# Irish Sea area. It was found via the CEDA Archive at the location of: +# https://catalogue.ceda.ac.uk/uuid/1b89e025eedd49e8976ee0721ec6e9b5, with +# DOI of https://dx.doi.org/10.5285/031e7ca1-9710-280d-e063-6c86abc014a0: f = cf.read( "~/summerstudents/final-recipes/new-required-datasets/POLCOMS_WAM_ZUV_01_16012006.nc") From 3b15542f5ba17a90e96af629dc4ee0d7d1553f6b Mon Sep 17 00:00:00 2001 From: "Sadie L. Bartholomew" Date: Fri, 11 Oct 2024 00:53:46 +0100 Subject: [PATCH 8/9] Add recipe 20 to recipes listing with filter keywords --- docs/source/recipes/plot_20_recipe.py | 3 +-- docs/source/recipes/recipe_list.txt | 8 ++++++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/docs/source/recipes/plot_20_recipe.py b/docs/source/recipes/plot_20_recipe.py index b88798b069..11c3250842 100644 --- a/docs/source/recipes/plot_20_recipe.py +++ b/docs/source/recipes/plot_20_recipe.py @@ -19,8 +19,7 @@ # Irish Sea area. It was found via the CEDA Archive at the location of: # https://catalogue.ceda.ac.uk/uuid/1b89e025eedd49e8976ee0721ec6e9b5, with # DOI of https://dx.doi.org/10.5285/031e7ca1-9710-280d-e063-6c86abc014a0: -f = cf.read( - "~/summerstudents/final-recipes/new-required-datasets/POLCOMS_WAM_ZUV_01_16012006.nc") +f = cf.read("~/recipes/POLCOMS_WAM_ZUV_01_16012006.nc") # %% # 3. Get the separate vector components, which are stored as separate fields. diff --git a/docs/source/recipes/recipe_list.txt b/docs/source/recipes/recipe_list.txt index f39297f9cb..d8061d0647 100644 --- a/docs/source/recipes/recipe_list.txt +++ b/docs/source/recipes/recipe_list.txt @@ -17,7 +17,7 @@ plot_08_recipe.html#sphx-glr-recipes-plot-08-recipe-py plot_09_recipe.html#sphx-glr-recipes-plot-09-recipe-py
plot_10_recipe.html#sphx-glr-recipes-plot-10-recipe-py -
+
plot_11_recipe.html#sphx-glr-recipes-plot-11-recipe-py
plot_12_recipe.html#sphx-glr-recipes-plot-12-recipe-py @@ -27,4 +27,8 @@ plot_13_recipe.html#sphx-glr-recipes-plot-13-recipe-py plot_14_recipe.html#sphx-glr-recipes-plot-14-recipe-py
plot_15_recipe.html#sphx-glr-recipes-plot-15-recipe-py -
\ No newline at end of file +
+plot_16_recipe.html#sphx-glr-recipes-plot-16-recipe-py +
+plot_20_recipe.html#sphx-glr-recipes-plot-20-recipe-py +
From 33a9debae364ca56b6c8077ac7817f5ced41437c Mon Sep 17 00:00:00 2001 From: "Sadie L. Bartholomew" Date: Fri, 11 Oct 2024 12:34:57 +0100 Subject: [PATCH 9/9] Add new filter keyword 'vectorplot' to appropriate recipes --- docs/source/recipes/recipe_list.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/recipes/recipe_list.txt b/docs/source/recipes/recipe_list.txt index d8061d0647..b03700c115 100644 --- a/docs/source/recipes/recipe_list.txt +++ b/docs/source/recipes/recipe_list.txt @@ -7,7 +7,7 @@ plot_03_recipe.html#sphx-glr-recipes-plot-03-recipe-py plot_04_recipe.html#sphx-glr-recipes-plot-04-recipe-py
plot_05_recipe.html#sphx-glr-recipes-plot-05-recipe-py -
+
plot_06_recipe.html#sphx-glr-recipes-plot-06-recipe-py
plot_07_recipe.html#sphx-glr-recipes-plot-07-recipe-py @@ -31,4 +31,4 @@ plot_15_recipe.html#sphx-glr-recipes-plot-15-recipe-py plot_16_recipe.html#sphx-glr-recipes-plot-16-recipe-py
plot_20_recipe.html#sphx-glr-recipes-plot-20-recipe-py -
+