Skip to content

Generalize EC cell widths function #526

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mark-petersen merged 4 commits into from
Apr 20, 2020
Merged

Generalize EC cell widths function #526

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a95a22d
Generalize EC cell widths function
xylar Apr 12, 2020
b3cc8a4
Improve plotting of cell widths in define_base_mesh
xylar Apr 12, 2020
964e31a
Tighten plotting and add more colormap support
xylar Apr 20, 2020
0307149
Add min and max to plot title
mark-petersen Apr 20, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 55 additions & 0 deletions testing_and_setup/compass/ocean/jigsaw_to_MPAS/3Wbgy5.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
<ColorMaps>
<ColorMap space="CIELAB" indexedLookup="false" name="3Wbgy5">
<Point x="0.000000" o="1" r="0.831373" g="0.909804" b="0.980392"/>
<Point x="0.012500" o="1" r="0.749020" g="0.862745" b="0.960784"/>
<Point x="0.025000" o="1" r="0.694118" g="0.827451" b="0.941176"/>

@mark-petersen mark-petersen Apr 20, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like you downloaded this directly from Paraview? That was smart.

@xylar xylar Apr 20, 2020

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, copied from MPAS-Analysis, which gives proper credit and got it from the original source (LANL's own SciVisColor web page, I believe).

<Point x="0.050000" o="1" r="0.568627" g="0.760784" b="0.921569"/>
<Point x="0.075000" o="1" r="0.450980" g="0.705882" b="0.901961"/>
<Point x="0.100000" o="1" r="0.345098" g="0.643137" b="0.858824"/>
<Point x="0.125000" o="1" r="0.247059" g="0.572549" b="0.819608"/>
<Point x="0.150000" o="1" r="0.180392" g="0.521569" b="0.780392"/>
<Point x="0.160000" o="1" r="0.149020" g="0.490196" b="0.749020"/>
<Point x="0.180000" o="1" r="0.129412" g="0.447059" b="0.709804"/>
<Point x="0.200000" o="1" r="0.101961" g="0.427451" b="0.690196"/>
<Point x="0.210000" o="1" r="0.094118" g="0.403922" b="0.658824"/>
<Point x="0.220000" o="1" r="0.090196" g="0.392157" b="0.639216"/>
<Point x="0.230000" o="1" r="0.082353" g="0.368627" b="0.619608"/>
<Point x="0.240000" o="1" r="0.070588" g="0.352941" b="0.600000"/>
<Point x="0.250000" o="1" r="0.066667" g="0.329412" b="0.568627"/>
<Point x="0.260000" o="1" r="0.074510" g="0.313725" b="0.541176"/>
<Point x="0.270000" o="1" r="0.086275" g="0.305882" b="0.509804"/>
<Point x="0.280000" o="1" r="0.094118" g="0.286275" b="0.478431"/>
<Point x="0.290000" o="1" r="0.101961" g="0.278431" b="0.450980"/>
<Point x="0.300000" o="1" r="0.109804" g="0.266667" b="0.411765"/>
<Point x="0.310000" o="1" r="0.113725" g="0.258824" b="0.380392"/>
<Point x="0.320000" o="1" r="0.113725" g="0.250980" b="0.349020"/>
<Point x="0.330000" o="1" r="0.109804" g="0.266667" b="0.321569"/>
<Point x="0.340000" o="1" r="0.105882" g="0.301961" b="0.262745"/>
<Point x="0.350000" o="1" r="0.094118" g="0.309804" b="0.243137"/>
<Point x="0.360000" o="1" r="0.082353" g="0.321569" b="0.227451"/>
<Point x="0.370000" o="1" r="0.074510" g="0.341176" b="0.219608"/>
<Point x="0.380000" o="1" r="0.070588" g="0.360784" b="0.211765"/>
<Point x="0.390000" o="1" r="0.066667" g="0.380392" b="0.215686"/>
<Point x="0.400000" o="1" r="0.062745" g="0.400000" b="0.176471"/>
<Point x="0.425000" o="1" r="0.074510" g="0.419608" b="0.145098"/>
<Point x="0.450000" o="1" r="0.086275" g="0.439216" b="0.117647"/>
<Point x="0.475000" o="1" r="0.121569" g="0.470588" b="0.117647"/>
<Point x="0.500000" o="1" r="0.184314" g="0.501961" b="0.149020"/>
<Point x="0.525000" o="1" r="0.254902" g="0.541176" b="0.188235"/>
<Point x="0.550000" o="1" r="0.325490" g="0.580392" b="0.231373"/>
<Point x="0.575000" o="1" r="0.403922" g="0.619608" b="0.278431"/>
<Point x="0.600000" o="1" r="0.501961" g="0.670588" b="0.333333"/>
<Point x="0.630000" o="1" r="0.592157" g="0.729412" b="0.400000"/>
<Point x="0.650000" o="1" r="0.741176" g="0.788235" b="0.490196"/>
<Point x="0.670000" o="1" r="0.858824" g="0.858824" b="0.603922"/>
<Point x="0.700000" o="1" r="0.921569" g="0.835294" b="0.580392"/>
<Point x="0.750000" o="1" r="0.901961" g="0.729412" b="0.494118"/>
<Point x="0.800000" o="1" r="0.858824" g="0.584314" b="0.388235"/>
<Point x="0.850000" o="1" r="0.800000" g="0.439216" b="0.321569"/>
<Point x="0.900000" o="1" r="0.678431" g="0.298039" b="0.203922"/>
<Point x="0.950000" o="1" r="0.549020" g="0.168627" b="0.109804"/>
<Point x="0.975000" o="1" r="0.478431" g="0.082353" b="0.047059"/>
<Point x="1.000000" o="1" r="0.450980" g="0.007843" b="0.000000"/>
<NaN r="0.25" g="0" b="0"/>
</ColorMap>
</ColorMaps>
69 changes: 54 additions & 15 deletions testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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()
Comment on lines +88 to +89

@xylar xylar Apr 20, 2020

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This cleans up some whitespace.

@pwolfram pwolfram Apr 20, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍

plt.savefig('cellWidthGlobal.png', bbox_inches='tight')

@xylar xylar Apr 20, 2020

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

bbox_inches='tight' eliminates some of the remaining white space around the image so it looks right.

@mark-petersen mark-petersen Apr 20, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, good to know.

plt.close()

else:
cellWidth, x, y, geom_points, geom_edges = define_base_mesh.cellWidthVsXY()
Expand Down Expand Up @@ -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')
Expand Down
114 changes: 56 additions & 58 deletions testing_and_setup/compass/ocean/jigsaw_to_MPAS/mesh_definition_tools.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down