Source code for compass.ocean.tests.global_ocean.mesh.wc14

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import mpas_tools.mesh.creation.mesh_definition_tools as mdt
import numpy as np
from geometric_features import read_feature_collection
from mpas_tools.cime.constants import constants
from mpas_tools.mesh.creation.signed_distance import (
    mask_from_geojson,
    signed_distance_from_geojson,
)
from mpas_tools.viz.colormaps import register_sci_viz_colormaps

from compass.mesh import QuasiUniformSphericalMeshStep


[docs] class WC14BaseMesh(QuasiUniformSphericalMeshStep): """ A step for creating WC14 mesh """ def setup(self): """ Add some input files """ inputs = ['coastline_CUSP.geojson', 'land_mask_Kamchatka.geojson', 'land_mask_Mexico.geojson', 'namelist.split_explicit_ab2', 'region_Arctic_Ocean.geojson', 'region_Bering_Sea.geojson', 'region_Bering_Sea_reduced.geojson', 'region_Central_America.geojson', 'region_Gulf_of_Mexico.geojson', 'region_Gulf_Stream_extension.geojson'] for filename in inputs: self.add_input_file(filename=filename, package=self.__module__) super().setup()
[docs] def build_cell_width_lat_lon(self): """ Create cell width array for this mesh on a regular latitude-longitude grid Returns ------- cellWidth : numpy.array m x n array of cell width in km lon : numpy.array longitude in degrees (length n and between -180 and 180) lat : numpy.array longitude in degrees (length m and between -90 and 90) """ dlon = 0.1 dlat = dlon earth_radius = constants['SHR_CONST_REARTH'] print(f'\nCreating cellWidth on a lat-lon grid of: ' f'{dlon:.2f} x {dlat:.2f} degrees') print('This can be set higher for faster test generation\n') nlon = int(360. / dlon) + 1 nlat = int(180. / dlat) + 1 lon = np.linspace(-180., 180., nlon) lat = np.linspace(-90., 90., nlat) km = 1.0e3 print('plotting ...') plt.switch_backend('Agg') fig = plt.figure() plt.clf() fig.set_size_inches(10.0, 14.0) register_sci_viz_colormaps() # Create cell width vs latitude for Atlantic and Pacific basins EC60to30 = mdt.EC_CellWidthVsLat(lat) EC60to30Narrow = mdt.EC_CellWidthVsLat(lat, latPosEq=8.0, latWidthEq=3.0) # Expand from 1D to 2D _, cellWidth = np.meshgrid(lon, EC60to30Narrow) _plot_cartopy(2, 'narrow EC60to30', cellWidth, '3Wbgy5') plotFrame = 3 # global settings for regionally refines mesh highRes = 14.0 # [km] fileName = 'region_Central_America' transitionWidth = 800.0 * km transitionOffset = 0.0 fc = read_feature_collection('{}.geojson'.format(fileName)) signedDistance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=0.25) mask = 0.5 * (1 + np.tanh((transitionOffset - signedDistance) / (transitionWidth / 2.))) cellWidth = 30.0 * mask + cellWidth * (1 - mask) fileName = 'coastline_CUSP' distanceToTransition = 600.0 * km # transitionWidth is distance from 0.07 to 0.03 of transition within # tanh transitionWidth = 600.0 * km transitionOffset = distanceToTransition + transitionWidth / 2.0 fc = read_feature_collection('{}.geojson'.format(fileName)) signedDistance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=0.25) mask = 0.5 * (1 + np.tanh((transitionOffset - signedDistance) / (transitionWidth / 2.))) cellWidth = highRes * mask + cellWidth * (1 - mask) _plot_cartopy(plotFrame, fileName + ' mask', mask, 'Blues') _plot_cartopy(plotFrame + 1, 'cellWidth ', cellWidth, '3Wbgy5') plotFrame += 2 fileName = 'region_Gulf_of_Mexico' transitionOffset = 600.0 * km transitionWidth = 600.0 * km fc = read_feature_collection('{}.geojson'.format(fileName)) signedDistance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=0.25) maskSmooth = 0.5 * (1 + np.tanh((transitionOffset - signedDistance) / (transitionWidth / 2.))) maskSharp = 0.5 * (1 + np.sign(-signedDistance)) fc = read_feature_collection('land_mask_Mexico.geojson') signedDistance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=0.25) landMask = 0.5 * (1 + np.sign(-signedDistance)) mask = maskSharp * landMask + maskSmooth * (1 - landMask) cellWidth = highRes * mask + cellWidth * (1 - mask) _plot_cartopy(plotFrame, fileName + ' mask', mask, 'Blues') _plot_cartopy(plotFrame + 1, 'cellWidth ', cellWidth, '3Wbgy5') plotFrame += 2 fileName = 'region_Bering_Sea' transitionOffset = 0.0 * km transitionWidth = 600.0 * km fc = read_feature_collection('{}.geojson'.format(fileName)) signedDistance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=0.25) maskSmoothEast = 0.5 * ( 1 + np.tanh((transitionOffset - signedDistance) / (transitionWidth / 2.))) fc = read_feature_collection('region_Bering_Sea_reduced.geojson') signedDistance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=0.25) maskSmoothWest = 0.5 * ( 1 + np.tanh((transitionOffset - signedDistance) / (transitionWidth / 2.))) fc = read_feature_collection('land_mask_Kamchatka.geojson') maskWest = mask_from_geojson(fc, lon, lat) mask = maskSmoothWest * maskWest + maskSmoothEast * (1 - maskWest) cellWidth = highRes * mask + cellWidth * (1 - mask) _plot_cartopy(plotFrame, fileName + ' mask', mask, 'Blues') _plot_cartopy(plotFrame + 1, 'cellWidth ', cellWidth, '3Wbgy5') plotFrame += 2 fileName = 'region_Arctic_Ocean' transitionOffset = 0.0 * km transitionWidth = 600.0 * km fc = read_feature_collection('{}.geojson'.format(fileName)) signedDistance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=0.25) mask = 0.5 * (1 + np.tanh((transitionOffset - signedDistance) / (transitionWidth / 2.))) cellWidth = highRes * mask + cellWidth * (1 - mask) _plot_cartopy(plotFrame, fileName + ' mask', mask, 'Blues') _plot_cartopy(plotFrame + 1, 'cellWidth ', cellWidth, '3Wbgy5') plotFrame += 2 fileName = 'region_Gulf_Stream_extension' transitionOffset = 0.0 * km transitionWidth = 600.0 * km fc = read_feature_collection('{}.geojson'.format(fileName)) signedDistance = signed_distance_from_geojson(fc, lon, lat, earth_radius, max_length=0.25) mask = 0.5 * (1 + np.tanh((transitionOffset - signedDistance) / (transitionWidth / 2.))) cellWidth = highRes * mask + cellWidth * (1 - mask) _plot_cartopy(plotFrame, fileName + ' mask', mask, 'Blues') _plot_cartopy(plotFrame + 1, 'cellWidth ', cellWidth, '3Wbgy5') plotFrame += 2 ax = plt.subplot(6, 2, 1) ax.plot(lat, EC60to30, label='original EC60to30') ax.plot(lat, EC60to30Narrow, label='narrow EC60to30') ax.grid(True) plt.title('Grid cell size [km] versus latitude') plt.legend(loc="upper left") plt.savefig('mesh_construction.png', dpi=300) return cellWidth, lon, lat
def _plot_cartopy(nPlot, varName, var, map_name): ax = plt.subplot(6, 2, nPlot, projection=ccrs.PlateCarree()) ax.set_global() im = ax.imshow(var, origin='lower', transform=ccrs.PlateCarree(), extent=[-180, 180, -90, 90], cmap=map_name, zorder=0) ax.add_feature(cfeature.LAND, edgecolor='black', zorder=1) gl = ax.gridlines( crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='-', zorder=2) ax.coastlines() gl.top_labels = False gl.bottom_labels = False gl.right_labels = False gl.left_labels = False plt.colorbar(im, shrink=.9) plt.title(varName)