From a95a22d9141410d1c9b03eaaf17e7c2a27295fd4 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 12 Apr 2020 13:30:09 +0200 Subject: [PATCH 1/4] Generalize EC cell widths function Many supposed parameters are actually hard-coded. Documentation has been updated to adhere to numpydoc formatting. Unnecessary looping has been removed. --- .../jigsaw_to_MPAS/mesh_definition_tools.py | 114 +++++++++--------- 1 file changed, 56 insertions(+), 58 deletions(-) diff --git a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/mesh_definition_tools.py b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/mesh_definition_tools.py index 395aefcb0a..cead245037 100755 --- a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/mesh_definition_tools.py +++ b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/mesh_definition_tools.py @@ -65,74 +65,72 @@ def mergeCellWidthVsLat( return cellWidthOut -def EC_CellWidthVsLat(lat): - ''' - EC_CellWidthVsLat - Create Eddy Closure spacing as a function of lat. - This is inted as part of the workflow to make an MPAS global mesh. +def EC_CellWidthVsLat(lat, cellWidthEq=30.0, cellWidthMidLat=60.0, + cellWidthPole=35.0, latPosEq=15.0, latPosPole=73.0, + latTransition=40.0, latWidthEq=6.0, latWidthPole=9.0): + """ + Create Eddy Closure spacing as a function of lat. This is intended as part + of the workflow to make an MPAS global mesh. - Syntax: cellWidthOut = EC_CellWidthVsLat(lat, cellWidthEq, cellWidthMidLat, cellWidthPole, - latPosEq, latPosPole, latTransition, - latWidthEq, latWidthPole) - Inputs: - lat - vector of length n, with entries between -90 and 90, degrees + Parameters + ---------- + lat : numpy.ndarray + vector of length n, with entries between -90 and 90, degrees - Optional inputs: - Default values for Cell width, km - cellWidthEq = 30.0 # Eq is equatorial lat - cellWidthMidLat = 60.0 # MidLat is mid lat - cellWidthPole = 35.0 # Pole is polar lat - - Default values for lat positions in degrees - latPosEq = 15.0 # position of center of transition region - latPosPole = 73.0 # position of center of transition region - latTransition = 40 # lat to change from Eq to Pole function - latWidthEq = 6.0 # width of transition region - latWidthPole = 9.0 # width of transition region + cellWidthEq : float, optional + Cell width in km at the equator - Outputs: - cellWidthOut - vector of length n, entrie are cell width as a function of lat + cellWidthMidLat : float, optional + Cell width in km at mid latitudes - Example: - EC60to30 = EC_CellWidthVsLat(lat) - EC120to60 = EC_CellWidthVsLat(lat,60,120,70) - ''' + cellWidthPole : float, optional + Cell width in km at the poles + + latPosEq : float, optional + Latitude in degrees of center of the equatorial transition region + + latPosPole : float, optional + Latitude in degrees of center of the polar transition region + + latTransition : float, optional + Latitude in degrees of the change from equatorial to polar function + + latWidthEq : float, optional + Width in degrees latitude of the equatorial transition region + + latWidthPole : float, optional + Width in degrees latitude of the polar transition region + + Returns + ------- + + cellWidthOut : numpy.ndarray + 1D array of same length as ``lat`` with entries that are cell width as + a function of lat + + Examples + -------- + Default + + >>> EC60to30 = EC_CellWidthVsLat(lat) + + Half the default resolution: - # Default values for Cell width, km - cellWidthEq = 30.0 # Eq is equatorial lat - cellWidthMidLat = 60.0 # MidLat is mid lat - cellWidthPole = 35.0 # Pole is polar lat - - # Default values for lat positions in degrees - latPosEq = 15.0 # position of center of transition region - latPosPole = 73.0 # position of center of transition region - latTransition = 40 # lat to change from Eq to Pole function - latWidthEq = 6.0 # width of transition region - latWidthPole = 9.0 # width of transition region - - # try - # cellWidthEq = varargin{1} # - # cellWidthMidLat = varargin{2} # - # cellWidthPole = varargin{3} # - # latPosEq = varargin{4} # - # latPosPole = varargin{5} # - # latTransition = varargin{6} # - # latWidthEq = varargin{7} # - # latWidthPole = varargin{8} # + >>> EC120to60 = EC_CellWidthVsLat(lat, cellWidthEq=60., cellWidthMidLat=120., cellWidthPole=70.) + """ minCellWidth = min(cellWidthEq, min(cellWidthMidLat, cellWidthPole)) densityEq = (minCellWidth / cellWidthEq)**4 densityMidLat = (minCellWidth / cellWidthMidLat)**4 densityPole = (minCellWidth / cellWidthPole)**4 - densityEC = np.zeros(lat.shape) - cellWidthOut = np.zeros(lat.shape) - for j in range(lat.size): - if np.abs(lat[j]) < latTransition: - densityEC[j] = ((densityEq - densityMidLat) * (1.0 + np.tanh( - (latPosEq - np.abs(lat[j])) / latWidthEq)) / 2.0) + densityMidLat - else: - densityEC[j] = ((densityMidLat - densityPole) * (1.0 + np.tanh( - (latPosPole - np.abs(lat[j])) / latWidthPole)) / 2.0) + densityPole - cellWidthOut[j] = minCellWidth / densityEC[j]**0.25 + densityEqToMid = ((densityEq - densityMidLat) * (1.0 + np.tanh( + (latPosEq - np.abs(lat)) / latWidthEq)) / 2.0) + densityMidLat + densityMidToPole = ((densityMidLat - densityPole) * (1.0 + np.tanh( + (latPosPole - np.abs(lat)) / latWidthPole)) / 2.0) + densityPole + mask = np.abs(lat) < latTransition + densityEC = np.array(densityMidToPole) + densityEC[mask] = densityEqToMid[mask] + cellWidthOut = minCellWidth / densityEC**0.25 return cellWidthOut From b3cc8a43414cda0eea58f699015abe2dad74cbb0 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 12 Apr 2020 13:32:35 +0200 Subject: [PATCH 2/4] Improve plotting of cell widths in define_base_mesh Imports should be at the top of the file. Colormap is now perceputally uniform (not "jet"!) Plot layout is now tighter (though there's still considerable whitespace. --- .../compass/ocean/jigsaw_to_MPAS/build_mesh.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py index 8d298a1193..8b7a627c87 100755 --- a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py +++ b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py @@ -12,10 +12,12 @@ from __future__ import absolute_import, division, print_function, \ unicode_literals -import os import xarray import argparse +import matplotlib +matplotlib.use('Agg') import matplotlib.pyplot as plt +import cartopy.crs as ccrs from mpas_tools.conversion import convert from mpas_tools.io import write_netcdf @@ -30,6 +32,7 @@ import define_base_mesh + def build_mesh( preserve_floodplain=False, floodplain_elevation=20.0, @@ -53,17 +56,14 @@ def build_mesh( da.to_netcdf(cw_filename) plot_cellWidth=True if plot_cellWidth: - import matplotlib - from cartopy import config - import cartopy.crs as ccrs - matplotlib.use('Agg') fig = plt.figure() fig.set_size_inches(16.0, 8.0) plt.clf() ax = plt.axes(projection=ccrs.PlateCarree()) ax.set_global() - im = ax.imshow(cellWidth, origin='lower', transform=ccrs.PlateCarree( - ), extent=[-180, 180, -90, 90], cmap='jet') + im = ax.imshow(cellWidth, origin='lower', + transform=ccrs.PlateCarree(), + extent=[-180, 180, -90, 90], cmap='viridis_r') ax.coastlines() gl = ax.gridlines( crs=ccrs.PlateCarree(), @@ -76,6 +76,8 @@ def build_mesh( gl.ylabels_right = False plt.title('Grid cell size, km') plt.colorbar(im, shrink=.60) + fig.canvas.draw() + plt.tight_layout() plt.savefig('cellWidthGlobal.png') else: From 964e31ae76f252e98989589df57ef351b004bd51 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 20 Apr 2020 19:00:39 +0200 Subject: [PATCH 3/4] Tighten plotting and add more colormap support Add support for a colormap form SciVisColor (3Wbgy5). --- .../compass/ocean/jigsaw_to_MPAS/3Wbgy5.xml | 55 +++++++++++++++++++ .../ocean/jigsaw_to_MPAS/build_mesh.py | 51 ++++++++++++++--- 2 files changed, 98 insertions(+), 8 deletions(-) create mode 100644 testing_and_setup/compass/ocean/jigsaw_to_MPAS/3Wbgy5.xml diff --git a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/3Wbgy5.xml b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/3Wbgy5.xml new file mode 100644 index 0000000000..ea7bde2fe3 --- /dev/null +++ b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/3Wbgy5.xml @@ -0,0 +1,55 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py index 8b7a627c87..effad7536c 100755 --- a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py +++ b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py @@ -18,6 +18,10 @@ matplotlib.use('Agg') import matplotlib.pyplot as plt import cartopy.crs as ccrs +import cartopy +import xml.etree.ElementTree as ET +import pkg_resources +from matplotlib.colors import LinearSegmentedColormap from mpas_tools.conversion import convert from mpas_tools.io import write_netcdf @@ -32,7 +36,6 @@ import define_base_mesh - def build_mesh( preserve_floodplain=False, floodplain_elevation=20.0, @@ -56,29 +59,34 @@ def build_mesh( da.to_netcdf(cw_filename) plot_cellWidth=True if plot_cellWidth: - fig = plt.figure() - fig.set_size_inches(16.0, 8.0) - plt.clf() + map_name = '3Wbgy5' + xmlFile = pkg_resources.resource_filename( + __name__, '{}.xml'.format(map_name)) + _read_xml_colormap(xmlFile, map_name) + + fig = plt.figure(figsize=[16.0, 8.0]) ax = plt.axes(projection=ccrs.PlateCarree()) ax.set_global() im = ax.imshow(cellWidth, origin='lower', transform=ccrs.PlateCarree(), - extent=[-180, 180, -90, 90], cmap='viridis_r') - ax.coastlines() + extent=[-180, 180, -90, 90], cmap=map_name, + zorder=0) + ax.add_feature(cartopy.feature.LAND, edgecolor='black', zorder=1) gl = ax.gridlines( crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color='gray', alpha=0.5, - linestyle='-') + linestyle='-', zorder=2) gl.xlabels_top = False gl.ylabels_right = False plt.title('Grid cell size, km') plt.colorbar(im, shrink=.60) fig.canvas.draw() plt.tight_layout() - plt.savefig('cellWidthGlobal.png') + plt.savefig('cellWidthGlobal.png', bbox_inches='tight') + plt.close() else: cellWidth, x, y, geom_points, geom_edges = define_base_mesh.cellWidthVsXY() @@ -132,6 +140,33 @@ def build_mesh( print("***********************************************") +def _read_xml_colormap(xmlFile, mapName): + """Read in an XML colormap""" + + xml = ET.parse(xmlFile) + + root = xml.getroot() + colormap = root.findall('ColorMap') + if len(colormap) > 0: + colormap = colormap[0] + colorDict = {'red': [], 'green': [], 'blue': []} + for point in colormap.findall('Point'): + x = float(point.get('x')) + color = [float(point.get('r')), float(point.get('g')), + float(point.get('b'))] + colorDict['red'].append((x, color[0], color[0])) + colorDict['green'].append((x, color[1], color[1])) + colorDict['blue'].append((x, color[2], color[2])) + cmap = LinearSegmentedColormap(mapName, colorDict, 256) + + _register_colormap_and_reverse(mapName, cmap) + + +def _register_colormap_and_reverse(mapName, cmap): + plt.register_cmap(mapName, cmap) + plt.register_cmap('{}_r'.format(mapName), cmap.reversed()) + + if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('--preserve_floodplain', action='store_true') From 0307149ca4310b7ac721e46cffb840bfa06e33a8 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Mon, 20 Apr 2020 14:24:08 -0600 Subject: [PATCH 4/4] Add min and max to plot title --- .../compass/ocean/jigsaw_to_MPAS/build_mesh.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py index effad7536c..23857fea8e 100755 --- a/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py +++ b/testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py @@ -57,7 +57,7 @@ def build_mesh( name='cellWidth') cw_filename = 'cellWidthVsLatLon.nc' da.to_netcdf(cw_filename) - plot_cellWidth=True + plot_cellWidth = True if plot_cellWidth: map_name = '3Wbgy5' xmlFile = pkg_resources.resource_filename( @@ -81,7 +81,9 @@ def build_mesh( linestyle='-', zorder=2) gl.xlabels_top = False gl.ylabels_right = False - plt.title('Grid cell size, km') + plt.title( + 'Grid cell size, km, min: {:.1f} max: {:.1f}'.format( + cellWidth.min(),cellWidth.max())) plt.colorbar(im, shrink=.60) fig.canvas.draw() plt.tight_layout()