matplotlib unstructered quadrilaterals instead of triangles

I've two netcdf files containing both unstructured grids. The first grid has 3 vertices per face and the second has 4 vertices per face.

For the grid containing 3 vertices per face I can use matplotlib.tri for visualization (like triplot_demo.py:

import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
xy = np.asarray([
    [-0.101, 0.872], [-0.080, 0.883], [-0.069, 0.888], [-0.054, 0.890],
    [-0.045, 0.897], [-0.057, 0.895], [-0.073, 0.900], [-0.087, 0.898],
    [-0.090, 0.904], [-0.069, 0.907], [-0.069, 0.921], [-0.080, 0.919],
    [-0.073, 0.928], [-0.052, 0.930], [-0.048, 0.942], [-0.062, 0.949],
    [-0.054, 0.958], [-0.069, 0.954], [-0.087, 0.952], [-0.087, 0.959],
    [-0.080, 0.966], [-0.085, 0.973], [-0.087, 0.965], [-0.097, 0.965],
    [-0.097, 0.975], [-0.092, 0.984], [-0.101, 0.980], [-0.108, 0.980],
    [-0.104, 0.987], [-0.102, 0.993], [-0.115, 1.001], [-0.099, 0.996],
    [-0.101, 1.007], [-0.090, 1.010], [-0.087, 1.021], [-0.069, 1.021],
    [-0.052, 1.022], [-0.052, 1.017], [-0.069, 1.010], [-0.064, 1.005],
    [-0.048, 1.005], [-0.031, 1.005], [-0.031, 0.996], [-0.040, 0.987],
    [-0.045, 0.980], [-0.052, 0.975], [-0.040, 0.973], [-0.026, 0.968],
    [-0.020, 0.954], [-0.006, 0.947], [ 0.003, 0.935], [ 0.006, 0.926],
    [ 0.005, 0.921], [ 0.022, 0.923], [ 0.033, 0.912], [ 0.029, 0.905],
    [ 0.017, 0.900], [ 0.012, 0.895], [ 0.027, 0.893], [ 0.019, 0.886],
    [ 0.001, 0.883], [-0.012, 0.884], [-0.029, 0.883], [-0.038, 0.879],
    [-0.057, 0.881], [-0.062, 0.876], [-0.078, 0.876], [-0.087, 0.872],
    [-0.030, 0.907], [-0.007, 0.905], [-0.057, 0.916], [-0.025, 0.933],
    [-0.077, 0.990], [-0.059, 0.993]])
x = np.degrees(xy[:, 0])
y = np.degrees(xy[:, 1])

triangles = np.asarray([
    [65, 44, 20],
    [65, 60, 44]])

triang = tri.Triangulation(x, y, triangles)

plt.figure()
plt.gca().set_aspect('equal')
plt.triplot(triang, 'go-', lw=1.0)
plt.title('triplot of user-specified triangulation')
plt.xlabel('Longitude (degrees)')
plt.ylabel('Latitude (degrees)')

plt.show()

-- indices of the related point annotated afterwards

BUT how to visualize the unstructured grid containing 4 vertices per face (quadrilaterals)? Following the previous exapmle, my faces looks like:

quatrang = np.asarray([
    [65, 60, 44, 20]])

Obviously trying tri.Triangulation doesn't work:

quatr = tri.Triangulation(x, y, quatrang)

ValueError: triangles must be a (?,3) array

I cannot find anything in the matplotlib libraries regarding 4 vertices per face. Any help is greatly appreciated..

EDIT: Changed the question based upon a minimal, complete and verifiable example

As commented already, since there is no Quatrangulation or simiar, there is no standard way to plot a a similar plot as triplot with four points per shape in matplotlib. Of course you could triangulate your mesh again to obtain 2 triangles per quadrilateral. Or, you can plot a PolyCollection of the shapes, given their coordinates in space. The following shows the latter, defining a quatplot function which takes the coordinates and the indices of the vertices as input and draws a PolyCollection of those to the axes.

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.collections

xy = np.asarray([
    [-0.101, 0.872], [-0.080, 0.883], [-0.069, 0.888], [-0.054, 0.890],
    [-0.090, 0.904], [-0.069, 0.907], [-0.069, 0.921], [-0.080, 0.919],
    [-0.080, 0.966], [-0.085, 0.973], [-0.087, 0.965], [-0.097, 0.965],
    [-0.104, 0.987], [-0.102, 0.993], [-0.115, 1.001], [-0.099, 0.996],
    [-0.052, 1.022], [-0.052, 1.017], [-0.069, 1.010], [-0.064, 1.005],
    [-0.045, 0.980], [-0.052, 0.975], [-0.040, 0.973], [-0.026, 0.968],
    [ 0.017, 0.900], [ 0.012, 0.895], [ 0.027, 0.893], [ 0.019, 0.886],
    [ 0.001, 0.883], [-0.012, 0.884], [-0.029, 0.883], [-0.038, 0.879],
    [-0.030, 0.907], [-0.007, 0.905], [-0.057, 0.916], [-0.025, 0.933],
    [-0.077, 0.990], [-0.059, 0.993]])
x = np.degrees(xy[:, 0])
y = np.degrees(xy[:, 1])

quatrang = np.asarray([
    [19,13,10,22], [35,7,3,28]])

def quatplot(x,y, quatrangles, ax=None, **kwargs):
    if not ax: ax=plt.gca()
    xy = np.c_[x,y]
    verts=xy[quatrangles]
    pc = matplotlib.collections.PolyCollection(verts, **kwargs)
    ax.add_collection(pc)
    ax.autoscale()


plt.figure()
plt.gca().set_aspect('equal')

quatplot(x,y, quatrang, ax=None, color="crimson", facecolor="None")
plt.plot(x,y, marker="o", ls="", color="crimson")

plt.title('quatplot of user-specified quatrangulation')
plt.xlabel('Longitude (degrees)')
plt.ylabel('Latitude (degrees)')

for i, (xi,yi) in enumerate(np.degrees(xy)):
    plt.text(xi,yi,i, size=8)

plt.show()