From 07f593c0c101f12ef5e84019264c3709ab47ef46 Mon Sep 17 00:00:00 2001 From: Jiaqi Li Date: Mon, 12 Jul 2021 10:44:17 -0600 Subject: [PATCH 1/8] initial commit --- Plots/Contours/NCL_polar_8.py | 117 ++++++++++++++++++++++++++++++++++ 1 file changed, 117 insertions(+) create mode 100755 Plots/Contours/NCL_polar_8.py diff --git a/Plots/Contours/NCL_polar_8.py b/Plots/Contours/NCL_polar_8.py new file mode 100755 index 000000000..cf4b78e29 --- /dev/null +++ b/Plots/Contours/NCL_polar_8.py @@ -0,0 +1,117 @@ +""" +NCL_polar_8.py +============== +This script illustrates the following concepts: + - Drawing filled contours and streamlines over a polar stereographic map + - Drawing the northern hemisphere of a polar stereographic map + +See following URLs to see the reproduced NCL plot & script: + - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/polar_8_lg.ncl + - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/polar_8_lg.png +""" +############################################################################### +# Import packages: + +import numpy as np +import xarray as xr +import cartopy.feature as cfeature +import cartopy.crs as ccrs +import matplotlib.pyplot as plt +import matplotlib.ticker as mticker + +import geocat.datafiles as gdf +from geocat.viz import util as gvutil + +############################################################################### +# Read in data: + +# Open a netCDF data file using xarray default engine and load the data into xarrays +ds = xr.open_dataset(gdf.get("netcdf_files/atmos.nc")) +u = ds.U[0, 1, :, :] +v = ds.V[0, 1, :, :] +t = ds.TS[0, :, :] + +# Extract slices of data + +############################################################################### +# Fix the artifact of not-shown-data around 0 and 360-degree longitudes +wrap_U = gvutil.xr_add_cyclic_longitudes(u, "lon") +wrap_V = gvutil.xr_add_cyclic_longitudes(v, "lon") +wrap_T = gvutil.xr_add_cyclic_longitudes(t, "lon") + +############################################################################### +# Plot: + +# Generate axes, using Cartopy, drawing coastlines, and adding features +fig = plt.figure(figsize=(10, 12)) +projection = ccrs.NorthPolarStereo() +ax = plt.axes(projection=projection) +ax.add_feature(cfeature.LAND, facecolor='lightgray') + +# Set map boundary to include latitudes between 0 and 40 and longitudes +# between -180 and 180 only +gvutil.set_map_boundary(ax, [-180, 180], [0, 40], south_pad=1) + +# Set draw_labels to False so that you can manually manipulate it later +gl = ax.gridlines(ccrs.PlateCarree(), + draw_labels=False, + linestyle="--", + color='black') + +# Manipulate latitude and longitude gridline numbers and spacing +gl.ylocator = mticker.FixedLocator(np.arange(0, 90, 15)) +gl.xlocator = mticker.FixedLocator(np.arange(-180, 180, 30)) + +# Manipulate longitude labels (0, 30 E, 60 E, ..., 30 W, etc.) +ticks = np.arange(0, 210, 30) +etick = ['0'] + [ + r'%dE' % tick for tick in ticks if (tick != 0) & (tick != 180) +] + ['180'] +wtick = [r'%dW' % tick for tick in ticks if (tick != 0) & (tick != 180)] +labels = etick + wtick +xticks = np.arange(0, 360, 30) +yticks = np.full_like(xticks, -5) # Latitude where the labels will be drawn + +for xtick, ytick, label in zip(xticks, yticks, labels): + if label == '180': + ax.text(xtick, + ytick, + label, + fontsize=14, + horizontalalignment='center', + verticalalignment='top', + transform=ccrs.Geodetic()) + elif label == '0': + ax.text(xtick, + ytick, + label, + fontsize=14, + horizontalalignment='center', + verticalalignment='bottom', + transform=ccrs.Geodetic()) + else: + ax.text(xtick, + ytick, + label, + fontsize=14, + horizontalalignment='center', + verticalalignment='center', + transform=ccrs.Geodetic()) + +# Contour-plot U-data +p = wrap_U.plot.contour(ax=ax, + vmin=-8, + vmax=16, + transform=ccrs.PlateCarree(), + levels=np.arange(-12, 44, 4), + linewidths=0.5, + cmap='black', + add_labels=False) + +# Use geocat.viz.util convenience function to add left and right titles +gvutil.set_titles_and_labels(ax, + lefttitle="Surface temperature", + righttitle="K") + +# Show the plot +plt.show() From f791334c35d9bb4d24677cd07c5f7f13bacf3656 Mon Sep 17 00:00:00 2001 From: Jiaqi Li Date: Mon, 12 Jul 2021 13:34:56 -0600 Subject: [PATCH 2/8] without vectors --- Plots/Contours/NCL_polar_8.py | 91 ++++++++++++++++++++++++----------- 1 file changed, 63 insertions(+), 28 deletions(-) diff --git a/Plots/Contours/NCL_polar_8.py b/Plots/Contours/NCL_polar_8.py index cf4b78e29..77ae11a7f 100755 --- a/Plots/Contours/NCL_polar_8.py +++ b/Plots/Contours/NCL_polar_8.py @@ -6,7 +6,7 @@ - Drawing the northern hemisphere of a polar stereographic map See following URLs to see the reproduced NCL plot & script: - - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/polar_8_lg.ncl + - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/polar_8.ncl - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/polar_8_lg.png """ ############################################################################### @@ -26,51 +26,55 @@ # Read in data: # Open a netCDF data file using xarray default engine and load the data into xarrays -ds = xr.open_dataset(gdf.get("netcdf_files/atmos.nc")) +ds = xr.open_dataset(gdf.get("netcdf_files/atmos.nc"), decode_times=False) u = ds.U[0, 1, :, :] v = ds.V[0, 1, :, :] t = ds.TS[0, :, :] # Extract slices of data +U = u.sel(lat=slice(60, 90)) +V = u.sel(lat=slice(60, 90)) +T = t.sel(lat=slice(60, 90)) ############################################################################### # Fix the artifact of not-shown-data around 0 and 360-degree longitudes -wrap_U = gvutil.xr_add_cyclic_longitudes(u, "lon") -wrap_V = gvutil.xr_add_cyclic_longitudes(v, "lon") -wrap_T = gvutil.xr_add_cyclic_longitudes(t, "lon") +wrap_U = gvutil.xr_add_cyclic_longitudes(U, "lon") +wrap_V = gvutil.xr_add_cyclic_longitudes(V, "lon") +wrap_T = gvutil.xr_add_cyclic_longitudes(T, "lon") ############################################################################### # Plot: -# Generate axes, using Cartopy, drawing coastlines, and adding features +# Generate axes with Cartopy projections fig = plt.figure(figsize=(10, 12)) projection = ccrs.NorthPolarStereo() ax = plt.axes(projection=projection) -ax.add_feature(cfeature.LAND, facecolor='lightgray') + +# Use Cartopy to add land feature +land_110m = cfeature.NaturalEarthFeature('physical', 'land', '110m') +ax.add_feature(land_110m, facecolor='white', edgecolor='black') # Set map boundary to include latitudes between 0 and 40 and longitudes # between -180 and 180 only -gvutil.set_map_boundary(ax, [-180, 180], [0, 40], south_pad=1) +gvutil.set_map_boundary(ax, [-180, 180], [60, 90], south_pad=1) -# Set draw_labels to False so that you can manually manipulate it later +# Set draw_labels to False to manually set labels later gl = ax.gridlines(ccrs.PlateCarree(), draw_labels=False, - linestyle="--", + linestyle=(0, (4, 10)), color='black') # Manipulate latitude and longitude gridline numbers and spacing -gl.ylocator = mticker.FixedLocator(np.arange(0, 90, 15)) +gl.ylocator = mticker.FixedLocator(np.arange(60, 90, 15)) gl.xlocator = mticker.FixedLocator(np.arange(-180, 180, 30)) # Manipulate longitude labels (0, 30 E, 60 E, ..., 30 W, etc.) -ticks = np.arange(0, 210, 30) -etick = ['0'] + [ - r'%dE' % tick for tick in ticks if (tick != 0) & (tick != 180) -] + ['180'] -wtick = [r'%dW' % tick for tick in ticks if (tick != 0) & (tick != 180)] +ticks = np.arange(30, 151, 30) +etick = ['0'] + [r'%dE' % tick for tick in ticks] + ['180'] +wtick = [r'%dW' % tick for tick in ticks] labels = etick + wtick -xticks = np.arange(0, 360, 30) -yticks = np.full_like(xticks, -5) # Latitude where the labels will be drawn +xticks = np.arange(0, 360, 30) # Longitude of the labels +yticks = np.full_like(xticks, 58) # Latitude of the labels for xtick, ytick, label in zip(xticks, yticks, labels): if label == '180': @@ -98,20 +102,51 @@ verticalalignment='center', transform=ccrs.Geodetic()) -# Contour-plot U-data -p = wrap_U.plot.contour(ax=ax, - vmin=-8, - vmax=16, - transform=ccrs.PlateCarree(), - levels=np.arange(-12, 44, 4), - linewidths=0.5, - cmap='black', - add_labels=False) +# Set contour levels +levels = np.arange(249, 283, 3) + +# Contourf-plot T-data +p = wrap_T.plot.contourf(ax=ax, + alpha=0.85, + transform=ccrs.PlateCarree(), + levels=levels, + cmap='viridis', + add_labels=False, + add_colorbar=False, + zorder=3) + +# Draw vector plot +# (there is no matplotlib equivalent to "CurlyVector" yet) +Q = ax.quiver(wrap_T['lon'], + wrap_T['lat'], + wrap_U.data, + wrap_V.data, + zorder=4, + pivot="middle", + width=0.001, + transform=ccrs.PlateCarree()) + +# Add colorbar +clb = plt.colorbar(p, + ax=ax, + pad=0.08, + shrink=0.85, + aspect=12, + ticks=levels, + extendrect=True, + extendfrac='auto', + orientation='horizontal') + +# Set colorbar ticks +clb.ax.xaxis.set_tick_params(length=0, labelsize=14, pad=9) # Use geocat.viz.util convenience function to add left and right titles gvutil.set_titles_and_labels(ax, lefttitle="Surface temperature", - righttitle="K") + righttitle="K", + lefttitlefontsize=22, + righttitlefontsize=22) # Show the plot +plt.tight_layout() plt.show() From da562acb3a20525a45828001391d88c603313fb3 Mon Sep 17 00:00:00 2001 From: Jiaqi Li Date: Mon, 12 Jul 2021 13:38:16 -0600 Subject: [PATCH 3/8] without vectors --- Plots/Contours/NCL_polar_8.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/Plots/Contours/NCL_polar_8.py b/Plots/Contours/NCL_polar_8.py index 77ae11a7f..f82aef8ef 100755 --- a/Plots/Contours/NCL_polar_8.py +++ b/Plots/Contours/NCL_polar_8.py @@ -117,14 +117,15 @@ # Draw vector plot # (there is no matplotlib equivalent to "CurlyVector" yet) -Q = ax.quiver(wrap_T['lon'], - wrap_T['lat'], - wrap_U.data, - wrap_V.data, - zorder=4, - pivot="middle", - width=0.001, - transform=ccrs.PlateCarree()) +# Q = ax.quiver(wrap_T['lon'], +# wrap_T['lat'], +# wrap_U.data, +# wrap_V.data, +# zorder=4, +# pivot="middle", +# width=0.001, +# scale=1000, +# transform=ccrs.PlateCarree()) # Add colorbar clb = plt.colorbar(p, From 1ab007cef509b56e0b07f0981ab880b3dde2fa17 Mon Sep 17 00:00:00 2001 From: Jiaqi Li Date: Thu, 15 Jul 2021 08:07:31 -0600 Subject: [PATCH 4/8] minor changes --- Plots/Contours/NCL_polar_8.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Plots/Contours/NCL_polar_8.py b/Plots/Contours/NCL_polar_8.py index f82aef8ef..35089db73 100755 --- a/Plots/Contours/NCL_polar_8.py +++ b/Plots/Contours/NCL_polar_8.py @@ -46,13 +46,13 @@ # Plot: # Generate axes with Cartopy projections -fig = plt.figure(figsize=(10, 12)) +fig = plt.figure(figsize=(8, 10)) projection = ccrs.NorthPolarStereo() ax = plt.axes(projection=projection) # Use Cartopy to add land feature land_110m = cfeature.NaturalEarthFeature('physical', 'land', '110m') -ax.add_feature(land_110m, facecolor='white', edgecolor='black') +ax.add_feature(land_110m, facecolor='none', edgecolor='gray') # Set map boundary to include latitudes between 0 and 40 and longitudes # between -180 and 180 only From fd149bb2f0bdb8f363ec9678201be40f5f6a58b8 Mon Sep 17 00:00:00 2001 From: Jiaqi Li Date: Thu, 15 Jul 2021 08:18:43 -0600 Subject: [PATCH 5/8] minor changes --- Plots/Contours/NCL_polar_8.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Plots/Contours/NCL_polar_8.py b/Plots/Contours/NCL_polar_8.py index 35089db73..3f04778d5 100755 --- a/Plots/Contours/NCL_polar_8.py +++ b/Plots/Contours/NCL_polar_8.py @@ -71,7 +71,7 @@ # Manipulate longitude labels (0, 30 E, 60 E, ..., 30 W, etc.) ticks = np.arange(30, 151, 30) etick = ['0'] + [r'%dE' % tick for tick in ticks] + ['180'] -wtick = [r'%dW' % tick for tick in ticks] +wtick = [r'%dW' % tick for tick in ticks[::-1]] labels = etick + wtick xticks = np.arange(0, 360, 30) # Longitude of the labels yticks = np.full_like(xticks, 58) # Latitude of the labels @@ -81,7 +81,7 @@ ax.text(xtick, ytick, label, - fontsize=14, + fontsize=13, horizontalalignment='center', verticalalignment='top', transform=ccrs.Geodetic()) @@ -89,7 +89,7 @@ ax.text(xtick, ytick, label, - fontsize=14, + fontsize=13, horizontalalignment='center', verticalalignment='bottom', transform=ccrs.Geodetic()) @@ -97,7 +97,7 @@ ax.text(xtick, ytick, label, - fontsize=14, + fontsize=13, horizontalalignment='center', verticalalignment='center', transform=ccrs.Geodetic()) @@ -130,23 +130,23 @@ # Add colorbar clb = plt.colorbar(p, ax=ax, - pad=0.08, + pad=0.12, shrink=0.85, - aspect=12, + aspect=9, ticks=levels, extendrect=True, extendfrac='auto', orientation='horizontal') # Set colorbar ticks -clb.ax.xaxis.set_tick_params(length=0, labelsize=14, pad=9) +clb.ax.xaxis.set_tick_params(length=0, labelsize=13, pad=9) # Use geocat.viz.util convenience function to add left and right titles gvutil.set_titles_and_labels(ax, lefttitle="Surface temperature", righttitle="K", - lefttitlefontsize=22, - righttitlefontsize=22) + lefttitlefontsize=16, + righttitlefontsize=16) # Show the plot plt.tight_layout() From b846bb11a92b73127c1640486efa680c7a125b6c Mon Sep 17 00:00:00 2001 From: Jiaqi Li Date: Mon, 2 Aug 2021 21:08:38 -0600 Subject: [PATCH 6/8] in progress --- Plots/Contours/NCL_polar_8.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/Plots/Contours/NCL_polar_8.py b/Plots/Contours/NCL_polar_8.py index 3f04778d5..fb13ea067 100755 --- a/Plots/Contours/NCL_polar_8.py +++ b/Plots/Contours/NCL_polar_8.py @@ -84,7 +84,7 @@ fontsize=13, horizontalalignment='center', verticalalignment='top', - transform=ccrs.Geodetic()) + transform=ccrs.PlateCarree()) elif label == '0': ax.text(xtick, ytick, @@ -92,7 +92,7 @@ fontsize=13, horizontalalignment='center', verticalalignment='bottom', - transform=ccrs.Geodetic()) + transform=ccrs.PlateCarree()) else: ax.text(xtick, ytick, @@ -100,7 +100,7 @@ fontsize=13, horizontalalignment='center', verticalalignment='center', - transform=ccrs.Geodetic()) + transform=ccrs.PlateCarree()) # Set contour levels levels = np.arange(249, 283, 3) @@ -117,15 +117,16 @@ # Draw vector plot # (there is no matplotlib equivalent to "CurlyVector" yet) -# Q = ax.quiver(wrap_T['lon'], -# wrap_T['lat'], -# wrap_U.data, -# wrap_V.data, -# zorder=4, -# pivot="middle", -# width=0.001, -# scale=1000, -# transform=ccrs.PlateCarree()) +Q = ax.quiver(wrap_U['lon'], + wrap_U['lat'], + wrap_U.data, + wrap_V.data, + zorder=4, + pivot="middle", + width=0.001, + color='white', + transform=ccrs.PlateCarree(), + regrid_shape=20) # Add colorbar clb = plt.colorbar(p, From 974239cad5d0b4ff3b6a1ed875dc1643a2c66e3e Mon Sep 17 00:00:00 2001 From: Katelyn FitzGerald <7872563+kafitzgerald@users.noreply.github.com> Date: Thu, 1 Feb 2024 12:22:04 -0700 Subject: [PATCH 7/8] move polar_8 example to new location and update imports --- {Plots => Gallery}/Contours/NCL_polar_8.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) rename {Plots => Gallery}/Contours/NCL_polar_8.py (90%) diff --git a/Plots/Contours/NCL_polar_8.py b/Gallery/Contours/NCL_polar_8.py similarity index 90% rename from Plots/Contours/NCL_polar_8.py rename to Gallery/Contours/NCL_polar_8.py index fb13ea067..aa9c138dd 100755 --- a/Plots/Contours/NCL_polar_8.py +++ b/Gallery/Contours/NCL_polar_8.py @@ -20,7 +20,7 @@ import matplotlib.ticker as mticker import geocat.datafiles as gdf -from geocat.viz import util as gvutil +import geocat.viz as gv ############################################################################### # Read in data: @@ -38,9 +38,9 @@ ############################################################################### # Fix the artifact of not-shown-data around 0 and 360-degree longitudes -wrap_U = gvutil.xr_add_cyclic_longitudes(U, "lon") -wrap_V = gvutil.xr_add_cyclic_longitudes(V, "lon") -wrap_T = gvutil.xr_add_cyclic_longitudes(T, "lon") +wrap_U = gv.xr_add_cyclic_longitudes(U, "lon") +wrap_V = gv.xr_add_cyclic_longitudes(V, "lon") +wrap_T = gv.xr_add_cyclic_longitudes(T, "lon") ############################################################################### # Plot: @@ -56,7 +56,7 @@ # Set map boundary to include latitudes between 0 and 40 and longitudes # between -180 and 180 only -gvutil.set_map_boundary(ax, [-180, 180], [60, 90], south_pad=1) +gv.set_map_boundary(ax, [-180, 180], [60, 90], south_pad=1) # Set draw_labels to False to manually set labels later gl = ax.gridlines(ccrs.PlateCarree(), @@ -143,11 +143,11 @@ clb.ax.xaxis.set_tick_params(length=0, labelsize=13, pad=9) # Use geocat.viz.util convenience function to add left and right titles -gvutil.set_titles_and_labels(ax, - lefttitle="Surface temperature", - righttitle="K", - lefttitlefontsize=16, - righttitlefontsize=16) +gv.set_titles_and_labels(ax, + lefttitle="Surface temperature", + righttitle="K", + lefttitlefontsize=16, + righttitlefontsize=16) # Show the plot plt.tight_layout() From d61e70b224a317f8709255a30926b33f4ca66fd9 Mon Sep 17 00:00:00 2001 From: Katelyn FitzGerald <7872563+kafitzgerald@users.noreply.github.com> Date: Thu, 1 Feb 2024 13:41:40 -0700 Subject: [PATCH 8/8] add rotating and rescaling to the polar quiver example --- Gallery/Contours/NCL_polar_8.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/Gallery/Contours/NCL_polar_8.py b/Gallery/Contours/NCL_polar_8.py index aa9c138dd..53bfdc132 100755 --- a/Gallery/Contours/NCL_polar_8.py +++ b/Gallery/Contours/NCL_polar_8.py @@ -25,21 +25,25 @@ ############################################################################### # Read in data: -# Open a netCDF data file using xarray default engine and load the data into xarrays +# Open a netCDF data file using xarray default engine ds = xr.open_dataset(gdf.get("netcdf_files/atmos.nc"), decode_times=False) -u = ds.U[0, 1, :, :] -v = ds.V[0, 1, :, :] -t = ds.TS[0, :, :] +U = ds.U[0, 1, :, :].sel(lat=slice(60, 90)) +V = ds.V[0, 1, :, :].sel(lat=slice(60, 90)) +T = ds.TS[0, :, :].sel(lat=slice(60, 90)) -# Extract slices of data -U = u.sel(lat=slice(60, 90)) -V = u.sel(lat=slice(60, 90)) -T = t.sel(lat=slice(60, 90)) +# Rotate and rescale the wind vectors following recommendations from SciTools/cartopy#1179 +U_source_crs = U / np.cos(U["lat"] / 180. * np.pi) +V_source_crs = V +magnitude = np.sqrt(U**2 + V**2) +magnitude_source_crs = np.sqrt(U_source_crs**2 + V_source_crs**2) + +U_rs = U_source_crs * magnitude / magnitude_source_crs +V_rs = V_source_crs * magnitude / magnitude_source_crs ############################################################################### # Fix the artifact of not-shown-data around 0 and 360-degree longitudes -wrap_U = gv.xr_add_cyclic_longitudes(U, "lon") -wrap_V = gv.xr_add_cyclic_longitudes(V, "lon") +wrap_U = gv.xr_add_cyclic_longitudes(U_rs, "lon") +wrap_V = gv.xr_add_cyclic_longitudes(V_rs, "lon") wrap_T = gv.xr_add_cyclic_longitudes(T, "lon") ###############################################################################