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 8d298a1193..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 @@ -12,10 +12,16 @@ 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 +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 @@ -51,32 +57,38 @@ def build_mesh( name='cellWidth') cw_filename = 'cellWidthVsLatLon.nc' da.to_netcdf(cw_filename) - plot_cellWidth=True + 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() + 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='jet') - ax.coastlines() + im = ax.imshow(cellWidth, origin='lower', + transform=ccrs.PlateCarree(), + 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.title( + 'Grid cell size, km, min: {:.1f} max: {:.1f}'.format( + cellWidth.min(),cellWidth.max())) plt.colorbar(im, shrink=.60) - plt.savefig('cellWidthGlobal.png') + fig.canvas.draw() + plt.tight_layout() + plt.savefig('cellWidthGlobal.png', bbox_inches='tight') + plt.close() else: cellWidth, x, y, geom_points, geom_edges = define_base_mesh.cellWidthVsXY() @@ -130,6 +142,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') 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