Contour irregular data within polygon

Question

I need to create filled contour plots of sea surface temperature (SST) data within a polygon, however I am not sure the best way to do this. I have three 1D arrays containing data for X, Y, and SST which I plot using the following to create the attached plot:

p=PatchCollection(mypatches,color='none', alpha=1.0,edgecolor="purple",linewidth=2.0)
levels=np.arange(SST.min(),SST.max(),0.2)
datamap=mymap.scatter(x,y,c=SST, s=55, vmin=-2,vmax=3,alpha=1.0)

I would like to be able to plot these data as filled contours (contourf instead of scatter) that are constrained (clipped) within the polygon boundaries (the purple line). Suggestions for how to achieve this are greatly appreciated.

Update: I had originally tried griddata, but could not get it to work properly. However, based on the answer provided by @eatHam I decided to try again. I could not get my scipy griddata to work as it kept hanging at the gridding when selecting method 'cubic', however when I switched to matplotlib.mlab.griddata and used 'linear' interpolation it worked. The suggestion for masking the boundaries provided a very coarse and not as exact solution as i would prefer.

I searched for options on how to clip contours in matplotlib and I found an answer by @pelson at this link. I tried the suggested solution implied in: "The contour set itself does not have a set_clip_path method but you can iterate over each of the contour collections and set their respective clip paths". My new and final solution looks like this (see plot below):

  p=PatchCollection(mypatches,color='none', alpha=1.0,edgecolor="purple",linewidth=2.0)
  levels=np.arange(SST.min(),SST.max(),0.2)
  grid_x, grid_y = np.mgrid[x.min()-0.5*(x.min()):x.max()+0.5*(x.max()):200j,
                          y.min()-0.5*(y.min()):y.max()+0.5*(y.max()):200j]
  grid_z = griddata(x,y,SST,grid_x,grid_y)

  cs=mymap.contourf(grid_x, grid_y, grid_z)

  for poly in mypatches:
      for artist in ax.get_children():
          artist.set_clip_path(poly)

      ax.add_patch(poly)
  mymap.drawcountries()
  mymap.drawcoastlines()
  mymap.fillcontinents(color='lightgrey',lake_color='none')
  mymap.drawmapboundary(fill_color='none')

This solution could also be improved particularily in terms of extrapolating the extreme edges in the North. Suggestions for how to really 'fill-in' the full polygon are appreciated. I also would like to understand why mlab worked and scipy not.

Acrospire

What interpolation methods have you tried? – Morph

the scipy and mlab version of griddata and different methods so that isn't surprising. I think the issue with interpolating outside of the data you have is more of a fundamental problem and will involve questionable methods to get around ;) – Morph

Answer 1

I would interpolate the data using scipy.griddata. You can set the region outside of your area (mypatches) to np.nan. And then just use pyplot.contour to plot it.

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

def sst_data(x, y):
    return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2

                         #replace with ...
x = np.random.rand(1000) #... your x
y = np.random.rand(1000) #... your y
sst = sst_data(x, y)     #... your sst

# interpolate to a grid
grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j] 
grid_z = griddata((x,y), sst, (grid_x, grid_y), method='cubic')

# mask out the area outside of your region
nr, nc = grid_z.shape
grid_z[-nr//3:, -nc//3:] = np.nan

plt.contourf(grid_x, grid_y, grid_z)

plt.show()

EDIT: Changed variable name in the plt.contourf() call (was ..(grid_z, grid_y, grid_z))

Blackmon

Your suggestion got me on the right track so I accepted your answer. – Acrospire