diff --git a/compass/meta.yaml b/compass/meta.yaml index 4fe58e1f9..6586de39e 100644 --- a/compass/meta.yaml +++ b/compass/meta.yaml @@ -1,5 +1,5 @@ {% set name = "compass" %} -{% set version = "0.1.9" %} +{% set version = "0.1.10" %} {% set build = 0 %} {% if mpi == "nompi" %} @@ -31,7 +31,7 @@ requirements: run: - python - geometric_features 0.1.11 - - mpas_tools 0.0.11 + - mpas_tools 0.0.12 - jigsaw 0.9.12 - jigsawpy 0.2.1 - metis diff --git a/conda_package/docs/api.rst b/conda_package/docs/api.rst index 1619cf9d6..084061284 100644 --- a/conda_package/docs/api.rst +++ b/conda_package/docs/api.rst @@ -28,20 +28,20 @@ MPAS mesh tools .. autosummary:: :toctree: generated/ - build_mesh.build_mesh - coastal_tools.coastal_refined_mesh - inject_bathymetry.inject_bathymetry - inject_meshDensity.inject_meshDensity - inject_preserve_floodplain.inject_preserve_floodplain + build_mesh.build_spherical_mesh + build_mesh.build_planar_mesh jigsaw_driver.jigsaw_driver + jigsaw_to_netcdf.jigsaw_to_netcdf mesh_definition_tools.mergeCellWidthVsLat mesh_definition_tools.EC_CellWidthVsLat mesh_definition_tools.RRS_CellWidthVsLat mesh_definition_tools.AtlanticPacificGrid mpas_to_triangle.mpas_to_triangle open_msh.readmsh + signed_distance.signed_distance_from_geojson + signed_distance.mask_from_geojson + signed_distance.distance_from_geojson triangle_to_netcdf.triangle_to_netcdf - jigsaw_to_netcdf.jigsaw_to_netcdf .. currentmodule:: mpas_tools.mesh.conversion @@ -99,6 +99,49 @@ Ocean Tools make_moc_basins_and_transects build_moc_basins +.. currentmodule:: mpas_tools.ocean.coastal_tools + +.. autosummary:: + :toctree: generated/ + + coastal_refined_mesh + create_background_mesh + extract_coastlines + distance_to_coast + compute_cell_width + save_matfile + CPP_projection + smooth_coastline + get_data_inside_box + get_indices_inside_quad + get_convex_hull_coordinates + plot_coarse_coast + plot_region_box + +.. currentmodule:: mpas_tools.ocean.inject_bathymetry + +.. autosummary:: + :toctree: generated/ + + inject_bathymetry + +.. currentmodule:: mpas_tools.ocean.inject_meshDensity + +.. autosummary:: + :toctree: generated/ + + inject_meshDensity_from_file + inject_spherical_meshDensity + inject_planar_meshDensity + +.. currentmodule:: mpas_tools.ocean.inject_preserve_floodplain + +.. autosummary:: + :toctree: generated/ + + inject_preserve_floodplain + + Visualization ============= diff --git a/conda_package/docs/define_base_mesh.py b/conda_package/docs/define_base_mesh.py deleted file mode 100755 index a9aa150fb..000000000 --- a/conda_package/docs/define_base_mesh.py +++ /dev/null @@ -1,21 +0,0 @@ -#!/usr/bin/env python -""" -% Create cell width array for this mesh on a regular latitude-longitude grid. -% Outputs: -% cellWidth - m x n array, entries are desired cell width in km -% lat - latitude, vector of length m, with entries between -90 and 90, degrees -% lon - longitude, vector of length n, with entries between -180 and 180, degrees -""" -import numpy as np - - -def cellWidthVsLatLon(): - - ddeg = 10 - constantCellWidth = 240 - - lat = np.arange(-90, 90.01, ddeg) - lon = np.arange(-180, 180.01, ddeg) - - cellWidth = constantCellWidth * np.ones((lat.size, lon.size)) - return cellWidth, lon, lat diff --git a/conda_package/docs/mesh_creation.rst b/conda_package/docs/mesh_creation.rst index 4fd444967..c1e5b4597 100644 --- a/conda_package/docs/mesh_creation.rst +++ b/conda_package/docs/mesh_creation.rst @@ -7,6 +7,8 @@ Mesh Creation Building a Mesh =============== +This is all out of date and needs refining!!! + The :py:func:`mpas_tools.mesh.creation.build_mesh.build_mesh` function is used create an MPAS mesh using the `JIGSAW `_ and `JIGSAW-Python `_ packages. diff --git a/conda_package/docs/versions.rst b/conda_package/docs/versions.rst index f3f9be704..dfa5ecaff 100644 --- a/conda_package/docs/versions.rst +++ b/conda_package/docs/versions.rst @@ -8,9 +8,12 @@ Documentation On GitHub ================ =============== `stable`_ `master`_ `v0.0.11`_ `0.0.11`_ +`v0.0.12`_ `0.0.12`_ ================ =============== .. _`stable`: ../stable/index.html .. _`master`: https://github.com/MPAS-Dev/MPAS-Tools/tree/master .. _`v0.0.11`: ../0.0.11/index.html -.. _`0.0.11`: https://github.com/MPAS-Dev/MPAS-Analysis/tree/0.0.11 \ No newline at end of file +.. _`0.0.11`: https://github.com/MPAS-Dev/MPAS-Analysis/tree/0.0.11 +.. _`v0.0.12`: ../0.0.12/index.html +.. _`0.0.12`: https://github.com/MPAS-Dev/MPAS-Analysis/tree/0.0.12 diff --git a/conda_package/mpas_tools/__init__.py b/conda_package/mpas_tools/__init__.py index 714773297..225b5f6d6 100644 --- a/conda_package/mpas_tools/__init__.py +++ b/conda_package/mpas_tools/__init__.py @@ -1,2 +1,2 @@ -__version_info__ = (0, 0, 11) +__version_info__ = (0, 0, 12) __version__ = '.'.join(str(vi) for vi in __version_info__) diff --git a/conda_package/mpas_tools/mesh/creation/__init__.py b/conda_package/mpas_tools/mesh/creation/__init__.py index 5653559e7..ca4707ba0 100644 --- a/conda_package/mpas_tools/mesh/creation/__init__.py +++ b/conda_package/mpas_tools/mesh/creation/__init__.py @@ -1,5 +1,6 @@ +from mpas_tools.mesh.creation.build_mesh import build_spherical_mesh, \ + build_planar_mesh from mpas_tools.mesh.creation.triangle_to_netcdf import triangle_to_netcdf from mpas_tools.mesh.creation.jigsaw_to_netcdf import jigsaw_to_netcdf -from mpas_tools.mesh.creation.inject_meshDensity import inject_meshDensity from mpas_tools.mesh.creation.mpas_to_triangle import mpas_to_triangle from mpas_tools.mesh.creation.open_msh import readmsh diff --git a/conda_package/mpas_tools/mesh/creation/build_mesh.py b/conda_package/mpas_tools/mesh/creation/build_mesh.py index 0e4ffa9e0..f851e116b 100644 --- a/conda_package/mpas_tools/mesh/creation/build_mesh.py +++ b/conda_package/mpas_tools/mesh/creation/build_mesh.py @@ -1,18 +1,7 @@ -""" -This script performs the first step of initializing the global ocean. This -includes: -Step 1. Build cellWidth array as function of latitude and longitude -Step 2. Build mesh using JIGSAW -Step 3. Convert triangles from jigsaw format to netcdf -Step 4. Convert from triangles to MPAS mesh -Step 5. Create vtk file for visualization -""" - from __future__ import absolute_import, division, print_function, \ unicode_literals import xarray -import argparse import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt @@ -21,193 +10,135 @@ from mpas_tools.mesh.conversion import convert from mpas_tools.io import write_netcdf -from mpas_tools.viz.paraview_extractor import extract_vtk from mpas_tools.mesh.creation.jigsaw_driver import jigsaw_driver from mpas_tools.mesh.creation.jigsaw_to_netcdf import jigsaw_to_netcdf -from mpas_tools.mesh.creation.inject_bathymetry import inject_bathymetry -from mpas_tools.mesh.creation.inject_meshDensity import inject_meshDensity -from mpas_tools.mesh.creation.inject_preserve_floodplain import \ - inject_preserve_floodplain from mpas_tools.viz.colormaps import register_sci_viz_colormaps -import sys -import os -# add the current working directory to the system path -sys.path.append(os.getcwd()) -import define_base_mesh -# okay, now we don't want to get anything else from CWD -del sys.path[-1] - -def build_mesh( - preserve_floodplain=False, - floodplain_elevation=20.0, - do_inject_bathymetry=False, - geometry='sphere', - plot_cellWidth=True): + +def build_spherical_mesh(cellWidth, lon, lat, earth_radius, + out_filename='base_mesh.nc', plot_cellWidth=True): """ Build an MPAS mesh using JIGSAW with the given cell sizes as a function of - latitude and longitude (on a sphere) or x and y (on a plane). + latitude and longitude. + + The result is a mesh file stored in ``out_filename`` as well as several + intermediate files: ``mesh.log``, ``mesh-HFUN.msh``, ``mesh.jig``, + ``mesh-MESH.msh``, ``mesh.msh``, and ``mesh_triangles.nc``. + + 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) + + earth_radius : float + Earth radius in meters - The user must define a local python module ``define_base_mesh`` that - provides a function that returns a 2D array ``cellWidth`` of cell sizes in - kilometers. + out_filename : str, optional + The file name of the resulting MPAS mesh - If ``geometry = 'sphere'``, this function is called ``cellWidthVsLatLon()`` - and also returns 1D ``lon`` and ``lat`` arrays. + plot_cellWidth : bool, optional + Whether to produce a plot of ``cellWidth``. If so, it will be written to + ``cellWidthGlobal.png``. + """ - If ``geometry = 'plane'`` (or any value other than `'sphere'``), the - function is called ``cellWidthVsXY()`` and returns 4 arrays in addition to - ``cellWidth``: 1D ``x`` and ``y`` arrays defining planar coordinates in - meters; as well as ``geom_points``, list of point coordinates for bounding - polygon for the planar mesh; and ``geom_edges``, list of edges between - points in ``geom_points`` that define the bounding polygon. + da = xarray.DataArray(cellWidth, + dims=['lat', 'lon'], + coords={'lat': lat, 'lon': lon}, + name='cellWidth') + cw_filename = 'cellWidthVsLatLon.nc' + da.to_netcdf(cw_filename) + if plot_cellWidth: + register_sci_viz_colormaps() + 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='3Wbgy5', + 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='-', zorder=2) + gl.top_labels = False + gl.right_labels = False + 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() + plt.savefig('cellWidthGlobal.png', bbox_inches='tight') + plt.close() + + print('Step 1. Generate mesh with JIGSAW') + jigsaw_driver(cellWidth, lon, lat, on_sphere=True, + earth_radius=earth_radius) + + print('Step 2. Convert triangles from jigsaw format to netcdf') + jigsaw_to_netcdf(msh_filename='mesh-MESH.msh', + output_name='mesh_triangles.nc', on_sphere=True) - The result is ``base_mesh.nc`` as well as several intermediate files: - ``mesh.log``, ``mesh-HFUN.msh``, ``mesh.jig``, ``mesh-MESH.msh``, - ``mesh.msh``, and ``mesh_triangles.nc``. + print('Step 3. Convert from triangles to MPAS mesh') + write_netcdf(convert(xarray.open_dataset('mesh_triangles.nc')), + out_filename) - The ``extract_vtk()`` function is used to produce a VTK file in the - ``base_mesh_vtk`` directory that can be viewed in ParaVeiw. + +def build_planar_mesh(cellWidth, x, y, geom_points, geom_edges, + out_filename='base_mesh.nc'): + """ + Build a planar MPAS mesh Parameters ---------- - preserve_floodplain : bool, optional - Whether a flood plain (bathymetry above z = 0) should be preserved in - the mesh. If so, a field ``cellSeedMask`` is added to the MPAS mesh - indicating positive elevations that should be preserved. + cellWidth : ndarray + m x n array of cell width in km - floodplain_elevation : float, optional - The elevation in meters to which the flood plain is preserved. + x, y : ndarray + arrays defining planar coordinates in meters - do_inject_bathymetry : bool, optional - Whether one of the default bathymetry datasets, ``earth_relief_15s.nc`` - or ``topo.msh``, should be added to the MPAS mesh in the field - ``bottomDepthObserved``. If so, a local link to one of these file names - must exist. + geom_points : ndarray + list of point coordinates for bounding polygon for the planar mesh - geometry : {'sphere', 'plane'}, optional - Whether the mesh is spherical or planar + geom_edges : ndarray + list of edges between points in ``geom_points`` that define the + bounding polygon - plot_cellWidth : bool, optional - If ``geometry = 'sphere'``, whether to produce a plot of ``cellWidth``. - If so, it will be written to ``cellWidthGlobal.png``. + out_filename : str, optional + The file name of the resulting MPAS mesh """ - - if geometry == 'sphere': - on_sphere = True - else: - on_sphere = False - - print('Step 1. Build cellWidth array as function of horizontal coordinates') - if on_sphere: - cellWidth, lon, lat = define_base_mesh.cellWidthVsLatLon() - da = xarray.DataArray(cellWidth, - dims=['lat', 'lon'], - coords={'lat': lat, 'lon': lon}, - name='cellWidth') - cw_filename = 'cellWidthVsLatLon.nc' - da.to_netcdf(cw_filename) - plot_cellWidth = True - if plot_cellWidth: - register_sci_viz_colormaps() - 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='3Wbgy5', - 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='-', zorder=2) - gl.top_labels = False - gl.right_labels = False - 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() - plt.savefig('cellWidthGlobal.png', bbox_inches='tight') - plt.close() - - else: - cellWidth, x, y, geom_points, geom_edges = define_base_mesh.cellWidthVsXY() - da = xarray.DataArray(cellWidth, - dims=['y', 'x'], - coords={'y': y, 'x': x}, - name='cellWidth') - cw_filename = 'cellWidthVsXY.nc' - da.to_netcdf(cw_filename) - - print('Step 2. Generate mesh with JIGSAW') - if on_sphere: - jigsaw_driver(cellWidth, lon, lat) - else: - jigsaw_driver( - cellWidth, - x, - y, - on_sphere=False, - geom_points=geom_points, - geom_edges=geom_edges) - - print('Step 3. Convert triangles from jigsaw format to netcdf') + da = xarray.DataArray(cellWidth, + dims=['y', 'x'], + coords={'y': y, 'x': x}, + name='cellWidth') + cw_filename = 'cellWidthVsXY.nc' + da.to_netcdf(cw_filename) + + print('Step 1. Generate mesh with JIGSAW') + jigsaw_driver( + cellWidth, + x, + y, + on_sphere=False, + geom_points=geom_points, + geom_edges=geom_edges) + + print('Step 2. Convert triangles from jigsaw format to netcdf') jigsaw_to_netcdf(msh_filename='mesh-MESH.msh', - output_name='mesh_triangles.nc', on_sphere=on_sphere) + output_name='mesh_triangles.nc', on_sphere=False) - print('Step 4. Convert from triangles to MPAS mesh') + print('Step 3. Convert from triangles to MPAS mesh') write_netcdf(convert(xarray.open_dataset('mesh_triangles.nc')), - 'base_mesh.nc') - - print('Step 5. Inject correct meshDensity variable into base mesh file') - inject_meshDensity(cw_filename=cw_filename, - mesh_filename='base_mesh.nc', on_sphere=on_sphere) - - if do_inject_bathymetry: - print('Step 6. Injecting bathymetry') - inject_bathymetry(mesh_file='base_mesh.nc') - - if preserve_floodplain: - print('Step 7. Injecting flag to preserve floodplain') - inject_preserve_floodplain(mesh_file='base_mesh.nc', - floodplain_elevation=floodplain_elevation) - - print('Step 8. Create vtk file for visualization') - extract_vtk(ignore_time=True, lonlat=True, dimension_list=['maxEdges='], - variable_list=['allOnCells'], filename_pattern='base_mesh.nc', - out_dir='base_mesh_vtk') - - print("***********************************************") - print("** The global mesh file is base_mesh.nc **") - print("***********************************************") - - -def main(): - parser = argparse.ArgumentParser() - parser.add_argument('--preserve_floodplain', action='store_true', - help='Whether a flood plain (bathymetry above z = 0) ' - 'should be preserved in the mesh') - parser.add_argument('--floodplain_elevation', action='store', - type=float, default=20.0, - help='The elevation in meters to which the flood plain ' - 'is preserved, default is 20 m') - parser.add_argument('--inject_bathymetry', action='store_true', - help='Whether one of the default bathymetry datasets, ' - 'earth_relief_15s.nc or topo.msh, should be added ' - 'to the MPAS mesh') - parser.add_argument('--geometry', default='sphere', - help='Whether the mesh is on a sphere or a plane, ' - 'default is a sphere') - parser.add_argument('--plot_cellWidth', action='store_true', - help='Whether to produce a plot of cellWidth') - cl_args = parser.parse_args() - build_mesh(cl_args.preserve_floodplain, cl_args.floodplain_elevation, - cl_args.inject_bathymetry, cl_args.geometry, - cl_args.plot_cellWidth) + out_filename) + diff --git a/conda_package/mpas_tools/mesh/creation/inject_meshDensity.py b/conda_package/mpas_tools/mesh/creation/inject_meshDensity.py deleted file mode 100644 index 5224e36d3..000000000 --- a/conda_package/mpas_tools/mesh/creation/inject_meshDensity.py +++ /dev/null @@ -1,82 +0,0 @@ -#!/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 -from scipy import interpolate -import netCDF4 as nc4 -import sys - - -def inject_meshDensity(cw_filename, mesh_filename, on_sphere=True): - 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'][:] - else: - x = ds.variables['x'][:] - y = ds.variables['y'][:] - ds.close() - - if on_sphere: - # Add extra column in longitude to interpolate over the Date Line - cellWidth = np.concatenate( - (cellWidth, cellWidth[:, 0:1]), axis=1) - LonPos = np.deg2rad(np.concatenate( - (lon.T, lon.T[0:1] + 360))) - LatPos = np.deg2rad(lat.T) - # set max lat position to be exactly at North Pole to avoid interpolation - # errors - LatPos[np.argmax(LatPos)] = np.pi / 2.0 - minCellWidth = cellWidth.min() - meshDensityVsXY = (minCellWidth / cellWidth)**4 - print(' minimum cell width in grid definition: {0:.0f} km'.format(minCellWidth/1000.)) - print(' maximum cell width in grid definition: {0:.0f} km'.format(cellWidth.max()/1000.)) - - if on_sphere: - X, Y = np.meshgrid(LonPos, LatPos) - else: - X, Y = np.meshgrid(x, y) - - print('Open unstructured MPAS mesh file...') - ds = nc4.Dataset(mesh_filename, 'r+') - meshDensity = ds.variables['meshDensity'] - - print('Preparing interpolation of meshDensity from native coordinates to mesh...') - meshDensityInterp = interpolate.LinearNDInterpolator( - np.vstack((X.ravel(), Y.ravel())).T, meshDensityVsXY.ravel()) - - print('Interpolating and writing meshDensity...') - if on_sphere: - meshDensity[:] = meshDensityInterp( - np.vstack( - (np.mod( - ds.variables['lonCell'][:] + - np.pi, - 2 * - np.pi) - - np.pi, - ds.variables['latCell'][:])).T) - else: - meshDensity[:] = meshDensityInterp( - np.vstack( - (ds.variables['xCell'][:], - ds.variables['yCell'][:]) - ).T) - - ds.close() - - -if __name__ == "__main__": - - inject_meshDensity(cw_filename=sys.argv[1], mesh_filename=sys.argv[2]) diff --git a/conda_package/mpas_tools/mesh/creation/jigsaw_driver.py b/conda_package/mpas_tools/mesh/creation/jigsaw_driver.py index 0960b27e1..a24393485 100644 --- a/conda_package/mpas_tools/mesh/creation/jigsaw_driver.py +++ b/conda_package/mpas_tools/mesh/creation/jigsaw_driver.py @@ -5,7 +5,8 @@ import jigsawpy -def jigsaw_driver(cellWidth, x, y, on_sphere=True, geom_points=None, geom_edges=None): +def jigsaw_driver(cellWidth, x, y, on_sphere=True, earth_radius=6371.0e3, + geom_points=None, geom_edges=None): ''' A function for building a jigsaw mesh @@ -15,14 +16,20 @@ def jigsaw_driver(cellWidth, x, y, on_sphere=True, geom_points=None, geom_edges= The size of each cell in the resulting mesh as a function of space x, y : ndarray - The x and y coordinates of each point in the cellWidth array (lon and lat for spherical mesh) + The x and y coordinates of each point in the cellWidth array (lon and + lat for spherical mesh) - on_sphere : logical + on_sphere : logical, optional Whether this mesh is spherical or planar - geom_points : list of point coordinates for bounding polygon for planar mesh + earth_radius : float, optional + Earth radius in meters - geom_edges : list of edges between points in geom_points that define the bounding polygon + geom_points : ndarray, optional + list of point coordinates for bounding polygon for planar mesh + + geom_edges : ndarray, optional + list of edges between points in geom_points that define the bounding polygon ''' # Authors @@ -53,7 +60,7 @@ def jigsaw_driver(cellWidth, x, y, on_sphere=True, geom_points=None, geom_edges= geom = jigsawpy.jigsaw_msh_t() if on_sphere: geom.mshID = 'ELLIPSOID-MESH' - geom.radii = 6371.*numpy.ones(3, float) + geom.radii = earth_radius*1e-3*numpy.ones(3, float) else: geom.mshID = 'EUCLIDEAN-MESH' geom.vert2 = geom_points diff --git a/conda_package/mpas_tools/mesh/creation/signed_distance.py b/conda_package/mpas_tools/mesh/creation/signed_distance.py new file mode 100644 index 000000000..f0debc69e --- /dev/null +++ b/conda_package/mpas_tools/mesh/creation/signed_distance.py @@ -0,0 +1,200 @@ +from __future__ import absolute_import, division, print_function, \ + unicode_literals + +import numpy as np +import pyflann +from scipy import spatial +import timeit +from rasterio.features import rasterize +from affine import Affine +import shapely.geometry +from geometric_features.plot import subdivide_geom +from mpas_tools.mesh.creation.util import lonlat2xyz + + +def signed_distance_from_geojson(fc, lon_grd, lat_grd, earth_radius, + max_length=None): + """ + Get the distance for each point on a lon/lat grid from the closest point + on the boundary of the geojson regions. + + Parameters + ---------- + fc : geometrics_features.FeatureCollection + The regions to be rasterized + lon_grd : numpy.ndarray + A 1D array of evenly spaced longitude values + lat_grd : numpy.ndarray + A 1D array of evenly spaced latitude values + earth_radius : float + Earth radius in meters + max_length : float, optional + The maximum distance (in degrees) between points on the boundary of the + geojson region. If the boundary is too coarse, it will be subdivided. + + Returns + ------- + signed_distance : numpy.ndarray + A 2D field of distances (negative inside the region, positive outside) + to the shape boundary + """ + distance = distance_from_geojson(fc, lon_grd, lat_grd, earth_radius, + nn_search='flann', max_length=max_length) + + mask = mask_from_geojson(fc, lon_grd, lat_grd) + + signed_distance = (-2.0 * mask + 1.0) * distance + return signed_distance + + +def mask_from_geojson(fc, lon_grd, lat_grd): + """ + Make a rasterized mask on a lon/lat grid from shapes (geojson multipolygon + data). + + Parameters + ---------- + fc : geometrics_features.FeatureCollection + The regions to be rasterized + lon_grd : numpy.ndarray + A 1D array of evenly spaced longitude values + lat_grd : numpy.ndarray + A 1D array of evenly spaced latitude values + + Returns + ------- + mask : numpy.ndarray + A 2D mask with the shapes rasterized (0.0 outside, 1.0 inside) + """ + + print("Mask from geojson") + print("-----------------") + + nlon = len(lon_grd) + nlat = len(lat_grd) + + dlon = (lon_grd[-1] - lon_grd[0])/(nlon-1) + dlat = (lat_grd[-1] - lat_grd[0])/(nlat-1) + + transform = Affine(dlon, 0.0, lon_grd[0], + 0.0, dlat, lat_grd[0]) + + shapes = [] + for feature in fc.features: + # a list of feature geometries and mask values (always 1.0) + shape = shapely.geometry.shape(feature['geometry']) + # expand a bit to make sure we hit the edges of the domain + shape = shape.buffer(dlat) + shapes.append((shapely.geometry.mapping(shape), 1.0)) + + mask = rasterize(shapes, out_shape=(nlat, nlon), transform=transform) + if np.abs(lon_grd[0] - (-180.)) < 1e-3 and \ + np.abs(lon_grd[-1] - (180.)) < 1e-3: + # the extra column at the periodic boundary needs to be copied + print(' fixing periodic boundary') + mask[:, -1] = mask[:, 0] + return mask + + +def distance_from_geojson(fc, lon_grd, lat_grd, earth_radius, nn_search, + max_length=None): + # {{{ + """ + Get the distance for each point on a lon/lat grid from the closest point + on the boundary of the geojson regions. + + Parameters + ---------- + fc : geometrics_features.FeatureCollection + The regions to be rasterized + lon_grd : numpy.ndarray + A 1D array of evenly spaced longitude values + lat_grd : numpy.ndarray + A 1D array of evenly spaced latitude values + earth_radius : float + Earth radius in meters + nn_search: {'kdtree', 'flann'} + The method used to find the nearest point on the shape boundary + max_length : float, optional + The maximum distance (in degrees) between points on the boundary of the + geojson region. If the boundary is too coarse, it will be subdivided. + + Returns + ------- + distance : numpy.ndarray + A 2D field of distances to the shape boundary + """ + print("Distance from geojson") + print("---------------------") + + print(" Finding region boundaries") + boundary_lon = [] + boundary_lat = [] + for feature in fc.features: + # get the boundary of each shape + shape = shapely.geometry.shape(feature['geometry']).boundary + if max_length is not None: + # subdivide the shape if it's too coarse + geom_type = shape.geom_type + shape = subdivide_geom(shape, geom_type, max_length) + x, y = shape.coords.xy + boundary_lon.extend(x) + boundary_lat.extend(y) + + boundary_lon = np.array(boundary_lon) + boundary_lat = np.array(boundary_lat) + + # Remove point at +/- 180 lon and +/- 90 lat because these are "fake". + # Need a little buffer (0.01 degrees) to avoid misses due to rounding. + mask = np.logical_not(np.logical_or( + np.logical_or(boundary_lon <= -179.99, boundary_lon >= 179.99), + np.logical_or(boundary_lat <= -89.99, boundary_lat >= 89.99))) + + boundary_lon = boundary_lon[mask] + boundary_lat = boundary_lat[mask] + + print(" Mean boundary latitude: {0:.2f}".format(np.mean(boundary_lat))) + + # Convert coastline points to x,y,z and create kd-tree + npoints = len(boundary_lon) + boundary_xyz = np.zeros((npoints, 3)) + boundary_xyz[:, 0], boundary_xyz[:, 1], boundary_xyz[:, 2] = \ + lonlat2xyz(boundary_lon, boundary_lat, earth_radius) + flann = None + tree = None + if nn_search == "kdtree": + tree = spatial.KDTree(boundary_xyz) + elif nn_search == "flann": + flann = pyflann.FLANN() + flann.build_index( + boundary_xyz, + algorithm='kdtree', + target_precision=1.0, + random_seed=0) + else: + raise ValueError('Bad nn_search: expected kdtree or flann, got ' + '{}'.format(nn_search)) + + # Convert backgound grid coordinates to x,y,z and put in a nx_grd x 3 array + # for kd-tree query + Lon_grd, Lat_grd = np.meshgrid(lon_grd, lat_grd) + X_grd, Y_grd, Z_grd = lonlat2xyz(Lon_grd, Lat_grd, earth_radius) + pts = np.vstack([X_grd.ravel(), Y_grd.ravel(), Z_grd.ravel()]).T + + # Find distances of background grid coordinates to the coast + print(" Finding distance") + start = timeit.default_timer() + distance = None + if nn_search == "kdtree": + distance, _ = tree.query(pts) + elif nn_search == "flann": + _, distance = flann.nn_index(pts, checks=2000, random_seed=0) + distance = np.sqrt(distance) + end = timeit.default_timer() + print(" Done") + print(" {0:.0f} seconds".format(end-start)) + + # Make distance array that corresponds with cell_width array + distance = np.reshape(distance, Lon_grd.shape) + + return distance diff --git a/conda_package/mpas_tools/mesh/creation/util.py b/conda_package/mpas_tools/mesh/creation/util.py index b9948127a..f8488ae4a 100644 --- a/conda_package/mpas_tools/mesh/creation/util.py +++ b/conda_package/mpas_tools/mesh/creation/util.py @@ -1,9 +1,20 @@ import collections +import numpy as np point = collections.namedtuple('Point', ['x', 'y', 'z']) -def circumcenter(on_sphere, x1, y1, z1, x2, y2, z2, x3, y3, z3): # {{{ +def circumcenter(on_sphere, x1, y1, z1, x2, y2, z2, x3, y3, z3): + """ + Compute the circumcenter of the triangle (possibly on a sphere) + with the three given vertices in Cartesian coordinates. + + Returns + ------- + center : point + The circumcenter of the triangle with x, y and z attributes + + """ p1 = point(x1, y1, z1) p2 = point(x2, y2, z2) p3 = point(x3, y3, z3) @@ -36,4 +47,30 @@ def circumcenter(on_sphere, x1, y1, z1, x2, y2, z2, x3, y3, z3): # {{{ # yv = yv / 3.0 return point(xv, yv, zv) -# }}} \ No newline at end of file + +def lonlat2xyz(lon, lat, R=6378206.4): + """ + Convert from longitude and latitude to Cartesian coordinates + + Parameters + ---------- + lon : float or numpy.ndarray + longitude + lat : float or numpy.ndarray + latitude + R : float, optional + Earth radius in meters + + Returns + ------- + x, y, z: float or numpy.array + Cartesian coordinates + """ + + lon = np.deg2rad(lon) + lat = np.deg2rad(lat) + x = R*np.multiply(np.cos(lon), np.cos(lat)) + y = R*np.multiply(np.sin(lon), np.cos(lat)) + z = R*np.sin(lat) + + return x, y, z diff --git a/conda_package/mpas_tools/mesh/interpolation.py b/conda_package/mpas_tools/mesh/interpolation.py new file mode 100644 index 000000000..32aae3268 --- /dev/null +++ b/conda_package/mpas_tools/mesh/interpolation.py @@ -0,0 +1,70 @@ +import numpy as np + + +def interp_bilin(x, y, field, xCell, yCell): + """ + Perform bilinear interpolation of ``field`` on a tensor grid to cell centers + on an MPAS mesh. ``xCell`` and ``yCell`` must be bounded by ``x`` and ``y``, + respectively. + + If x and y coordinates are longitude and latitude, respectively, it is + recommended that they be passed in degrees to avoid round-off problems at + the north and south poles and at the date line. + + Parameters + ---------- + x : ndarray + x coordinate of the input field (length n) + + y : ndarray + y coordinate fo the input field (length m) + + field : ndarray + a field of size m x n + + xCell : ndarray + x coordinate of MPAS cell centers + + yCell : ndarray + y coordinate of MPAS cell centers + + Returns + ------- + mpasField : ndarray + ``field`` interpoyed to MPAS cell centers + """ + + assert np.all(xCell >= x[0]) + assert np.all(xCell <= x[-1]) + assert np.all(yCell >= y[0]) + assert np.all(yCell <= y[-1]) + + # find float indices into the x and y arrays of cells on the MPAS mesh + xFrac = np.interp(xCell, x, np.arange(len(x))) + yFrac = np.interp(yCell, y, np.arange(len(y))) + + # xIndices/yIndices are the integer indices of the lower bound for bilinear + # interpoyion; xFrac/yFrac are the fraction of the way ot the next index + xIndices = np.array(xFrac, dtype=int) + xFrac -= xIndices + yIndices = np.array(yFrac, dtype=int) + yFrac -= yIndices + + # If points are exactly at the upper index, this is going to give us a bit + # of trouble so we'll move them down one index and adjust the fraction + # accordingly + mask = xIndices == len(x) + xIndices[mask] -= 1 + xFrac[mask] += 1. + + mask = yIndices == len(y) + yIndices[mask] -= 1 + yFrac[mask] += 1. + + mpasField = \ + (1. - xFrac) * (1. - yFrac) * field[yIndices, xIndices] + \ + xFrac * (1. - yFrac) * field[yIndices, xIndices + 1] + \ + (1. - xFrac) * yFrac * field[yIndices + 1, xIndices] + \ + xFrac * yFrac * field[yIndices + 1, xIndices + 1] + + return mpasField diff --git a/conda_package/mpas_tools/ocean/__init__.py b/conda_package/mpas_tools/ocean/__init__.py index e69de29bb..1fc1dcb22 100644 --- a/conda_package/mpas_tools/ocean/__init__.py +++ b/conda_package/mpas_tools/ocean/__init__.py @@ -0,0 +1,7 @@ +from mpas_tools.ocean.build_mesh import build_spherical_mesh, \ + build_planar_mesh +from mpas_tools.ocean.inject_bathymetry import inject_bathymetry +from mpas_tools.ocean.inject_meshDensity import inject_meshDensity_from_file, \ + inject_spherical_meshDensity, inject_planar_meshDensity +from mpas_tools.ocean.inject_preserve_floodplain import \ + inject_preserve_floodplain diff --git a/conda_package/mpas_tools/ocean/build_mesh.py b/conda_package/mpas_tools/ocean/build_mesh.py new file mode 100644 index 000000000..08b432a06 --- /dev/null +++ b/conda_package/mpas_tools/ocean/build_mesh.py @@ -0,0 +1,146 @@ +from __future__ import absolute_import, division, print_function, \ + unicode_literals + +from mpas_tools.mesh.creation import build_spherical_mesh as create_spherical_mesh +from mpas_tools.mesh.creation import build_planar_mesh as create_planar_mesh +from mpas_tools.cime.constants import constants +from mpas_tools.ocean.inject_bathymetry import inject_bathymetry +from mpas_tools.ocean.inject_meshDensity import inject_spherical_meshDensity, \ + inject_planar_meshDensity +from mpas_tools.ocean.inject_preserve_floodplain import inject_preserve_floodplain +from mpas_tools.viz.paraview_extractor import extract_vtk + + +def build_spherical_mesh(cellWidth, lon, lat, out_filename='base_mesh.nc', + plot_cellWidth=True, vtk_dir='base_mesh_vtk', + preserve_floodplain=False, floodplain_elevation=20.0, + do_inject_bathymetry=False): + """ + Build an MPAS mesh using JIGSAW with the given cell sizes as a function of + latitude and longitude + + The result is a mesh file stored in ``out_filename`` as well as several + intermediate files: ``mesh.log``, ``mesh-HFUN.msh``, ``mesh.jig``, + ``mesh-MESH.msh``, ``mesh.msh``, and ``mesh_triangles.nc``. + + 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) + + out_filename : str, optional + The file name of the resulting MPAS mesh + + plot_cellWidth : bool, optional + Whether to produce a plot of ``cellWidth``. If so, it will be written to + ``cellWidthGlobal.png``. + + vtk_dir : str, optional + The name of the directory where mesh data will be extracted for viewing + in ParaVeiw. + + preserve_floodplain : bool, optional + Whether a flood plain (bathymetry above z = 0) should be preserved in + the mesh + + floodplain_elevation : float, optional + The elevation in meters to which the flood plain is preserved + + do_inject_bathymetry : bool, optional + Whether one of the default bathymetry datasets, ``earth_relief_15s.nc`` + or ``topo.msh``, should be added to the MPAS mesh + """ + + # from https://github.com/E3SM-Project/E3SM/blob/master/cime/src/share/util/shr_const_mod.F90#L20 + # earth_radius = constants['SHR_CONST_REARTH'] + earth_radius = 6371.0e3 + + create_spherical_mesh(cellWidth, lon, lat, earth_radius, out_filename, + plot_cellWidth) + + print('Step 4. Inject meshDensity into the mesh file') + inject_spherical_meshDensity(cellWidth, lon, lat, + mesh_filename=out_filename) + + _shared_steps(out_filename, vtk_dir, preserve_floodplain, + floodplain_elevation, do_inject_bathymetry) + + +def build_planar_mesh(cellWidth, x, y, geom_points, geom_edges, + out_filename='base_mesh.nc', vtk_dir='base_mesh_vtk', + preserve_floodplain=False, floodplain_elevation=20.0, + do_inject_bathymetry=False): + """ + Build a planar MPAS mesh + + Parameters + ---------- + cellWidth : ndarray + m x n array of cell width in km + + x, y : ndarray + arrays defining planar coordinates in meters + + geom_points : ndarray + list of point coordinates for bounding polygon for the planar mesh + + geom_edges : ndarray + list of edges between points in ``geom_points`` that define the + bounding polygon + + out_filename : str, optional + The file name of the resulting MPAS mesh + + vtk_dir : str, optional + The name of the directory where mesh data will be extracted for viewing + in ParaVeiw. + + preserve_floodplain : bool, optional + Whether a flood plain (bathymetry above z = 0) should be preserved in + the mesh + + floodplain_elevation : float, optional + The elevation in meters to which the flood plain is preserved + + do_inject_bathymetry : bool, optional + Whether one of the default bathymetry datasets, ``earth_relief_15s.nc`` + or ``topo.msh``, should be added to the MPAS mesh + """ + create_planar_mesh(cellWidth, x, y, geom_points, geom_edges, out_filename) + + print('Step 4. Inject meshDensity into the mesh file') + inject_planar_meshDensity(cellWidth, x, y, mesh_filename=out_filename) + + _shared_steps(out_filename, vtk_dir, preserve_floodplain, + floodplain_elevation, do_inject_bathymetry) + + +def _shared_steps(out_filename, vtk_dir, preserve_floodplain, + floodplain_elevation, do_inject_bathymetry): + step = 5 + + if do_inject_bathymetry: + print('Step {}. Injecting bathymetry'.format(step)) + inject_bathymetry(mesh_file=out_filename) + step += 1 + + if preserve_floodplain: + print('Step {}. Injecting flag to preserve floodplain'.format(step)) + inject_preserve_floodplain(mesh_file=out_filename, + floodplain_elevation=floodplain_elevation) + step += 1 + + print('Step {}. Create vtk file for visualization'.format(step)) + extract_vtk(ignore_time=True, lonlat=True, dimension_list=['maxEdges='], + variable_list=['allOnCells'], filename_pattern=out_filename, + out_dir=vtk_dir) + + print("***********************************************") + print("** The global mesh file is {} **".format(out_filename)) + print("***********************************************") diff --git a/conda_package/mpas_tools/mesh/creation/coastal_tools.py b/conda_package/mpas_tools/ocean/coastal_tools.py similarity index 82% rename from conda_package/mpas_tools/mesh/creation/coastal_tools.py rename to conda_package/mpas_tools/ocean/coastal_tools.py index 7fbaceebe..ec8e4a7dd 100644 --- a/conda_package/mpas_tools/mesh/creation/coastal_tools.py +++ b/conda_package/mpas_tools/ocean/coastal_tools.py @@ -16,18 +16,11 @@ import matplotlib.pyplot as plt from scipy import spatial, io import timeit -import os -import cartopy import cartopy.crs as ccrs import cartopy.feature as cfeature -from rasterio.features import rasterize -from affine import Affine -import shapely.geometry -from geometric_features.plot import subdivide_geom import mpas_tools.mesh.creation.mesh_definition_tools as mdt +from mpas_tools.mesh.creation.util import lonlat2xyz plt.switch_backend('agg') -cartopy.config['pre_existing_data_dir'] = \ - os.getenv('CARTOPY_DIR', cartopy.config.get('pre_existing_data_dir')) # Constants km = 1000.0 @@ -336,10 +329,10 @@ def create_background_mesh(grd_box, ddeg, mesh_type, dx_min, dx_max, # {{{ print("------------------------") # Create cell width background grid - lat_grd = np.arange(grd_box[2], grd_box[3], ddeg) - lon_grd = np.arange(grd_box[0], grd_box[1], ddeg) - nx_grd = lon_grd.size - ny_grd = lat_grd.size + ny_grd = int((grd_box[3]-grd_box[2])/ddeg) + 1 + nx_grd = int((grd_box[1]-grd_box[0])/ddeg) + 1 + lat_grd = grd_box[2] + ddeg*np.arange(ny_grd) + lon_grd = grd_box[0] + ddeg*np.arange(nx_grd) print(" Background grid dimensions:", ny_grd, nx_grd) # Assign background grid cell width values @@ -349,6 +342,8 @@ def create_background_mesh(grd_box, ddeg, mesh_type, dx_min, dx_max, # {{{ cell_width_lat = mdt.EC_CellWidthVsLat(lat_grd) elif mesh_type == 'RRS': cell_width_lat = mdt.RRS_CellWidthVsLat(lat_grd, dx_max / km, dx_min / km) + else: + raise ValueError('Unknown mesh_type {}'.format(mesh_type)) cell_width = np.tile(cell_width_lat, (nx_grd, 1)).T * km # Plot background cell width @@ -566,188 +561,6 @@ def distance_to_coast(coastlines, lon_grd, lat_grd, origin, nn_search, smooth_wi ############################################################## -def signed_distance_from_geojson(fc, lon_grd, lat_grd, max_length=None): # {{{ - """ - Get the distance for each point on a lon/lat grid from the closest point - on the boundary of the geojson regions. - - Parameters - ---------- - fc : geometrics_features.FeatureCollection - The regions to be rasterized - lon_grd : numpy.ndarray - A 1D array of evenly spaced longitude values - lat_grd : numpy.ndarray - A 1D array of evenly spaced latitude values - max_length : float, optional - The maximum distance (in degrees) between points on the boundary of the - geojson region. If the boundary is too coarse, it will be subdivided. - - Returns - ------- - signed_distance : numpy.ndarray - A 2D field of distances (negative inside the region, positive outside) - to the shape boundary - """ - distance = distance_from_geojson(fc, lon_grd, lat_grd, - nn_search='flann', max_length=max_length) - - mask = mask_from_geojson(fc, lon_grd, lat_grd) - - signed_distance = (-2.0 * mask + 1.0) * distance - return signed_distance - - -def mask_from_geojson(fc, lon_grd, lat_grd): # {{{ - """ - Make a rasterized mask on a lon/lat grid from shapes (geojson multipolygon - data). - - Parameters - ---------- - fc : geometrics_features.FeatureCollection - The regions to be rasterized - lon_grd : numpy.ndarray - A 1D array of evenly spaced longitude values - lat_grd : numpy.ndarray - A 1D array of evenly spaced latitude values - - Returns - ------- - mask : numpy.ndarray - A 2D mask with the shapes rasterized (0.0 outside, 1.0 inside) - """ - - print("Mask from geojson") - print("-----------------") - - nlon = len(lon_grd) - nlat = len(lat_grd) - - dlon = (lon_grd[-1] - lon_grd[0])/(nlon-1) - dlat = (lat_grd[-1] - lat_grd[0])/(nlat-1) - - transform = Affine(dlon, 0.0, lon_grd[0], - 0.0, dlat, lat_grd[0]) - - shapes = [] - for feature in fc.features: - # a list of feature geometries and mask values (always 1.0) - shape = shapely.geometry.shape(feature['geometry']) - # expand a bit to make sure we hit the edges of the domain - shape = shape.buffer(dlat) - shapes.append((shapely.geometry.mapping(shape), 1.0)) - - mask = rasterize(shapes, out_shape=(nlat, nlon), transform=transform) - if np.abs(lon_grd[0] - (-180.)) < 1e-3 and \ - np.abs(lon_grd[-1] - (180.)) < 1e-3: - # the extra column at the periodic boundary needs to be copied - print(' fixing periodic boundary') - mask[:, -1] = mask[:, 0] - return mask # }}} - - -def distance_from_geojson(fc, lon_grd, lat_grd, nn_search, max_length=None): - # {{{ - """ - Get the distance for each point on a lon/lat grid from the closest point - on the boundary of the geojson regions. - - Parameters - ---------- - fc : geometrics_features.FeatureCollection - The regions to be rasterized - lon_grd : numpy.ndarray - A 1D array of evenly spaced longitude values - lat_grd : numpy.ndarray - A 1D array of evenly spaced latitude values - nn_search: {'kdtree', 'flann'} - The method used to find the nearest point on the shape boundary - max_length : float, optional - The maximum distance (in degrees) between points on the boundary of the - geojson region. If the boundary is too coarse, it will be subdivided. - - Returns - ------- - distance : numpy.ndarray - A 2D field of distances to the shape boundary - """ - print("Distance from geojson") - print("---------------------") - - print(" Finding region boundaries") - boundary_lon = [] - boundary_lat = [] - for feature in fc.features: - # get the boundary of each shape - shape = shapely.geometry.shape(feature['geometry']).boundary - if max_length is not None: - # subdivide the shape if it's too coarse - geom_type = shape.geom_type - shape = subdivide_geom(shape, geom_type, max_length) - x, y = shape.coords.xy - boundary_lon.extend(x) - boundary_lat.extend(y) - - boundary_lon = np.array(boundary_lon) - boundary_lat = np.array(boundary_lat) - - # Remove point at +/- 180 lon and +/- 90 lat because these are "fake". - # Need a little buffer (0.01 degrees) to avoid misses due to rounding. - mask = np.logical_not(np.logical_or( - np.logical_or(boundary_lon <= -179.99, boundary_lon >= 179.99), - np.logical_or(boundary_lat <= -89.99, boundary_lat >= 89.99))) - - boundary_lon = boundary_lon[mask] - boundary_lat = boundary_lat[mask] - - print(" Mean boundary latitude: {0:.2f}".format(np.mean(boundary_lat))) - - # Convert coastline points to x,y,z and create kd-tree - npoints = len(boundary_lon) - boundary_xyz = np.zeros((npoints, 3)) - boundary_xyz[:, 0], boundary_xyz[:, 1], boundary_xyz[:, 2] = \ - lonlat2xyz(boundary_lon, boundary_lat) - flann = None - tree = None - if nn_search == "kdtree": - tree = spatial.KDTree(boundary_xyz) - elif nn_search == "flann": - flann = pyflann.FLANN() - flann.build_index( - boundary_xyz, - algorithm='kdtree', - target_precision=1.0, - random_seed=0) - else: - raise ValueError('Bad nn_search: expected kdtree or flann, got ' - '{}'.format(nn_search)) - - # Convert backgound grid coordinates to x,y,z and put in a nx_grd x 3 array - # for kd-tree query - Lon_grd, Lat_grd = np.meshgrid(lon_grd, lat_grd) - X_grd, Y_grd, Z_grd = lonlat2xyz(Lon_grd, Lat_grd) - pts = np.vstack([X_grd.ravel(), Y_grd.ravel(), Z_grd.ravel()]).T - - # Find distances of background grid coordinates to the coast - print(" Finding distance") - start = timeit.default_timer() - distance = None - if nn_search == "kdtree": - distance, _ = tree.query(pts) - elif nn_search == "flann": - _, distance = flann.nn_index(pts, checks=2000, random_seed=0) - distance = np.sqrt(distance) - end = timeit.default_timer() - print(" Done") - print(" {0:.0f} seconds".format(end-start)) - - # Make distance array that corresponds with cell_width array - distance = np.reshape(distance, Lon_grd.shape) - - return distance # }}} - - def compute_cell_width(D, cell_width, lon, lat, dx_min, trans_start, trans_width, restrict_box, # {{{ plot_option=False, plot_box=[], lon_grd=[], lat_grd=[], coastlines=[], call=None): @@ -862,17 +675,6 @@ def CPP_projection(lon, lat, origin): ############################################################## -def lonlat2xyz(lon,lat): - - R = 6378206.4 - x = R*np.multiply(np.cos(lon*deg2rad),np.cos(lat*deg2rad)) - y = R*np.multiply(np.sin(lon*deg2rad),np.cos(lat*deg2rad)) - z = R*np.sin(lat*deg2rad) - - return x,y,z - -############################################################## - def smooth_coastline(x, y, window): xs = np.copy(x) diff --git a/conda_package/mpas_tools/mesh/creation/inject_bathymetry.py b/conda_package/mpas_tools/ocean/inject_bathymetry.py similarity index 92% rename from conda_package/mpas_tools/mesh/creation/inject_bathymetry.py rename to conda_package/mpas_tools/ocean/inject_bathymetry.py index 5aa73d284..f3b577c4a 100644 --- a/conda_package/mpas_tools/mesh/creation/inject_bathymetry.py +++ b/conda_package/mpas_tools/ocean/inject_bathymetry.py @@ -5,6 +5,7 @@ unicode_literals from mpas_tools.mesh.creation.open_msh import readmsh +from mpas_tools.mesh.interpolation import interp_bilin import numpy as np from scipy import interpolate import netCDF4 as nc4 @@ -84,13 +85,8 @@ def interpolate_SRTM(lon_pts, lat_pts): idx = np.intersect1d(lon_idx, lat_idx) xpts = lon_pts[idx] ypts = lat_pts[idx] - xy_pts = np.vstack((xpts, ypts)).T - # Interpolate bathymetry onto points - bathy = interpolate.RegularGridInterpolator( - (xdata, ydata), zdata.T, bounds_error=False, fill_value=np.nan) - bathy_int = bathy(xy_pts) - bathymetry[idx] = bathy_int + bathymetry[idx] = interp_bilin(xdata, ydata, zdata, xpts, ypts) end = timeit.default_timer() print(end - start, " seconds") diff --git a/conda_package/mpas_tools/ocean/inject_meshDensity.py b/conda_package/mpas_tools/ocean/inject_meshDensity.py new file mode 100644 index 000000000..844fdee07 --- /dev/null +++ b/conda_package/mpas_tools/ocean/inject_meshDensity.py @@ -0,0 +1,141 @@ +#!/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 + + +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) + + Returns + ------- + + """ + 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) + + + +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() + + +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]) diff --git a/conda_package/mpas_tools/mesh/creation/inject_preserve_floodplain.py b/conda_package/mpas_tools/ocean/inject_preserve_floodplain.py similarity index 100% rename from conda_package/mpas_tools/mesh/creation/inject_preserve_floodplain.py rename to conda_package/mpas_tools/ocean/inject_preserve_floodplain.py diff --git a/conda_package/recipe/meta.yaml b/conda_package/recipe/meta.yaml index 6f7592067..ece711440 100644 --- a/conda_package/recipe/meta.yaml +++ b/conda_package/recipe/meta.yaml @@ -1,5 +1,5 @@ {% set name = "mpas_tools" %} -{% set version = "0.0.11" %} +{% set version = "0.0.12" %} package: name: {{ name|lower }} @@ -15,9 +15,8 @@ build: - translate_planar_grid = mpas_tools.translate:main - merge_grids = mpas_tools.merge_grids:main - split_grids = mpas_tools.split_grids:main - - build_mesh = mpas_tools.mesh.creation.build_mesh:main - - inject_bathymetry = mpas_tools.mesh.creation.inject_bathymetry:main - - inject_preserve_floodplain = mpas_tools.mesh.creation.inject_preserve_floodplain:main + - inject_bathymetry = mpas_tools.ocean.inject_bathymetry:main + - inject_preserve_floodplain = mpas_tools.ocean.inject_preserve_floodplain:main - mpas_to_triangle = mpas_tools.mesh.creation.mpas_to_triangle:main - triangle_to_netcdf = mpas_tools.mesh.creation.triangle_to_netcdf:main - jigsaw_to_netcdf = mpas_tools.mesh.creation.jigsaw_to_netcdf:main @@ -97,9 +96,6 @@ test: - paraview_vtk_field_extractor.py -f mesh_tools/mesh_conversion_tools/test/mesh.QU.1920km.151026.nc -v latCell,lonCell --ignore_time -o vtk_test - split_grids --help - merge_grids --help - # define_base_mesh.py must exist locally for build_mesh to work - - cp conda_package/mpas_tools/tests/define_base_mesh.py . - - build_mesh --help - inject_bathymetry mesh_tools/mesh_conversion_tools/test/mesh.QU.1920km.151026.nc - inject_preserve_floodplain --help - mpas_to_triangle --help diff --git a/conda_package/setup.py b/conda_package/setup.py index fb5ab3aa4..068767977 100755 --- a/conda_package/setup.py +++ b/conda_package/setup.py @@ -60,9 +60,8 @@ 'translate_planar_grid = mpas_tools.translate:main', 'merge_grids = mpas_tools.merge_grids:main', 'split_grids = mpas_tools.split_grids:main', - 'build_mesh = mpas_tools.mesh.creation.build_mesh:main', - 'inject_bathymetry = mpas_tools.mesh.creation.inject_bathymetry:main', - 'inject_preserve_floodplain = mpas_tools.mesh.creation.inject_preserve_floodplain:main', + 'inject_bathymetry = mpas_tools.ocean.inject_bathymetry:main', + 'inject_preserve_floodplain = mpas_tools.ocean.inject_preserve_floodplain:main', 'mpas_to_triangle = mpas_tools.mesh.creation.mpas_to_triangle:main', 'triangle_to_netcdf = mpas_tools.mesh.creation.triangle_to_netcdf:main', 'jigsaw_to_netcdf = mpas_tools.mesh.creation.jigsaw_to_netcdf:main']})