Source code for mpas_tools.ocean.inject_meshDensity

#!/usr/bin/env python
# Simple script to inject mesh density onto a mesh
# example usage:
#   ./inject_meshDensity.py cellWidthVsLatLon.nc base_mesh.nc
# where:
#   cellWidthVsLatLon.nc is a netcdf file with cellWidth
#   base_mesh.nc is the mpas netcdf file where meshDensity is added
# Mark Petersen, 7/24/2018

from __future__ import absolute_import, division, print_function, \
    unicode_literals

import numpy as np
import netCDF4 as nc4
import sys

from mpas_tools.mesh.interpolation import interp_bilin


[docs]def inject_meshDensity_from_file(cw_filename, mesh_filename, on_sphere=True): """ Add a ``meshDensity`` field into an MPAS mesh. The mesh density is defined as: meshDensity = (minCellWidth / cellWidth)**4 Parameters ---------- cw_filename : str The file name to read ``cellWidth`` and coordinates from mesh_filename : str The mesh file to add ``meshDensity`` to on_sphere : bool, optional Whether the mesh is spherical (as opposed to planar) """ print('Read cell width field from nc file regular grid...') ds = nc4.Dataset(cw_filename,'r') cellWidth = ds.variables['cellWidth'][:] if on_sphere: lon = ds.variables['lon'][:] lat = ds.variables['lat'][:] ds.close() inject_spherical_meshDensity(cellWidth, lon, lat, mesh_filename) else: x = ds.variables['x'][:] y = ds.variables['y'][:] ds.close() inject_spherical_meshDensity(cellWidth, x, y, mesh_filename)
[docs]def inject_spherical_meshDensity(cellWidth, lon, lat, mesh_filename): """ Add a ``meshDensity`` field into a spherical MPAS mesh. The mesh density is defined as: meshDensity = (minCellWidth / cellWidth)**4 Parameters ---------- cellWidth : ndarray m x n array of cell width in km lon : ndarray longitude in degrees (length n and between -180 and 180) lat : ndarray longitude in degrees (length m and between -90 and 90) mesh_filename : str The mesh file to add ``meshDensity`` to """ minCellWidth = cellWidth.min() meshDensityVsXY = (minCellWidth / cellWidth)**4 print(' minimum cell width in grid definition: {0:.0f} km'.format( minCellWidth)) print(' maximum cell width in grid definition: {0:.0f} km'.format( cellWidth.max())) print('Open unstructured MPAS mesh file...') ds = nc4.Dataset(mesh_filename, 'r+') meshDensity = ds.variables['meshDensity'] lonCell = ds.variables['lonCell'][:] latCell = ds.variables['latCell'][:] lonCell = np.mod(np.rad2deg(lonCell) + 180., 360.) - 180. latCell = np.rad2deg(latCell) print('Interpolating and writing meshDensity...') mpasMeshDensity = interp_bilin(lon, lat, meshDensityVsXY, lonCell, latCell) meshDensity[:] = mpasMeshDensity ds.close()
[docs]def inject_planar_meshDensity(cellWidth, x, y, mesh_filename): """ Add a ``meshDensity`` field into a planar MPAS mesh. The mesh density is defined as: meshDensity = (minCellWidth / cellWidth)**4 Parameters ---------- cellWidth : ndarray m x n array of cell width in km x, y : ndarray Planar coordinates in meters mesh_filename : str The mesh file to add ``meshDensity`` to """ minCellWidth = cellWidth.min() meshDensityVsXY = (minCellWidth / cellWidth)**4 print(' minimum cell width in grid definition: {0:.0f} km'.format(minCellWidth)) print(' maximum cell width in grid definition: {0:.0f} km'.format(cellWidth.max())) print('Open unstructured MPAS mesh file...') ds = nc4.Dataset(mesh_filename, 'r+') meshDensity = ds.variables['meshDensity'] xCell = ds.variables['xCell'][:] yCell = ds.variables['xCell'][:] print('Interpolating and writing meshDensity...') mpasMeshDensity = interp_bilin(x, y, meshDensityVsXY, xCell, yCell) meshDensity[:] = mpasMeshDensity ds.close()
if __name__ == "__main__": inject_meshDensity_from_file(cw_filename=sys.argv[1], mesh_filename=sys.argv[2])