Last week, at the SciPy conference, I had chance to sprint on the pyugrid. The pyugrid module aims to be a Python API to utilize data written using the netCDF unstructured grid conventions. Right now it support only triangular meshes, but the conventions cover any type of unstructured grids.
Let's explore this module in this post. First: how to read an unstructured?
(And once hexagonal meshes are implemented I will be able to re-create Borges' Library of Babel and solved the crime using network theory and the connectivity array ;-)
import pyugrid
print(pyugrid.__version__)
url = ('http://comt.sura.org/thredds/dodsC/data/comt_1_archive/'
'inundation_tropical/VIMS_SELFE/Hurricane_Ike_3D_final_run_with_waves')
ugrid = pyugrid.UGrid.from_ncfile(url)
How about we check this dataset for consistency?
ugrid.check_consistent()
Oh, oh... I guess we should have worked harder during the sprint =)
OK, since this is a triangular mesh lets read read the topology into something
matplotlib users know (x, y, triangles).
lon = ugrid.nodes[:, 0]
lat = ugrid.nodes[:, 1]
triangles = ugrid.faces[:]
The dataset has nodes (or vertexes) representing a point in a 2D space, the faces (or polygons) correspond to a plane enclosed by a set of edges. Note that the dataset does not contain the edges! We need to compute those:
%time ugrid.build_edges()
~4 seconds is kind of slow. Since this grid has a 2D triangular mesh topology one might ask: how does the calculated edges compare with matplotlib triangulation? Let's check it out!
import matplotlib.tri as tri
%time triang = tri.Triangulation(lon, lat, triangles=triangles)
matplotlib create the triang object faster than pyugrid.
Let's take a closer look to see if the computed edges are similar.
ugrid.edges.shape == triang.edges.shape
ugrid.edges == triang.edges
OK. The edge matrix is different :-(, that is probably just the ordering of the elements though.
Let's take a quick look at the mesh.
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
def make_map(projection=ccrs.PlateCarree()):
fig, ax = plt.subplots(figsize=(8, 6),
subplot_kw=dict(projection=projection))
ax.coastlines(resolution='50m')
gl = ax.gridlines(draw_labels=True)
gl.xlabels_top = gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
return fig, ax
fig, ax = make_map()
kw = dict(marker='.', linestyle='-', alpha=0.25, color='darkgray')
ax.triplot(triang, **kw) # or lon, lat, triangules
ax.coastlines()
ax.set_extent([-84, -78, 25, 32])
Let's try to find the point nearest to one of the SECOORA buoys in the dataset.
duck = [-75.51997, 35.8112]
%%time
duck_face = ugrid.locate_face_simple(duck)
print(duck_face)
Ouch! That took way too long. Let's see if matplotlib is faster.
%%time
f = triang.get_trifinder()
print(f(duck[0], duck[1]))
Yep! matplotlib is not only faster but subsequent calls will instantaneous since the TriFinder object is cached.
How about finding the nodes?
%time duck_node = ugrid.locate_nodes(duck)
That was much faster! Locating nodes uses a cached KDTree for optimization.
Both pyugrid and matplotlib can return the face-face connectivity array (the neighbors of each triangle). The connectivity array can be extremely useful when searching for the shortest path between to grid points. Again, lets check if they are similar.
ugrid.build_face_face_connectivity()
(ugrid.face_face_connectivity == triang.neighbors).all()
neighbors = triang.neighbors[duck_face]
neighbors
This triangle has only two neighbors. Let's check if this is correct by plotting everything we found before ending this post.
import numpy as np
def plt_triangle(triang, face, ax=None, **kw):
if not ax:
fig, ax = plt.subplots()
ax.triplot(triang.x[triang.triangles[face]],
triang.y[triang.triangles[face]],
triangles=triang.triangles[face], **kw)
fig, ax = make_map()
ax.triplot(triang, **kw)
ax.plot(duck[0], duck[1], color='cornflowerblue',
marker='s', zorder=2)
ax.plot(lon[duck_node], lat[duck_node], color='crimson',
marker='o', zorder=2)
# Station.
ax.tripcolor(triang.x[triang.triangles[duck_face]], triang.y[triang.triangles[duck_face]],
triang.triangles[duck_face], facecolors=np.array([1]), alpha=0.25, zorder=1)
# Neighbors.
plt_triangle(triang, neighbors[0], ax=ax, color='darkolivegreen', zorder=1)
plt_triangle(triang, neighbors[-1], ax=ax, color='darkolivegreen', zorder=1)
ax.coastlines()
ax.set_extent([-76, -75.2, 35.5, 36])
Everything looks about right. We found the node (red circle) and the triangle (blue) nearest our buoy (blue square). We also found the neighboring triangles (green).
The best feature of pyugrid is that the library interprets the UGrid conventions for you when reading/saving unstructured grids (see this example). However, when it comes to manipulating the grid, pyugrid is lagging behind matplotlib.
BTW: matplotlib's TriAnalyzer packs some neat methods to analyze the triangular mesh.
HTML(html)
This post was written as an IPython notebook. It is available for download or as a static html.
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.