code: def twoSurfplots(self):
      
      # Could more elegantly just call surfplot twice
      # And also could include xyzinterp as an option inside surfplot.
      # Noted here in case anyone wants to take that on in the future...
  
      plt.subplot(211)
      plt.title('Load thickness, mantle equivalent [m]',fontsize=16)
      if self.latlon:
        plt.imshow(self.qs/(self.rho_m*self.g), extent=(0, self.dx*self.qs.shape[0], self.dy*self.qs.shape[1], 0))
        plt.xlabel('longitude [deg E]', fontsize=12, fontweight='bold')
        plt.ylabel('latitude [deg N]', fontsize=12, fontweight='bold')
      else:
        plt.imshow(self.qs/(self.rho_m*self.g), extent=(0, self.dx/1000.*self.qs.shape[0], self.dy/1000.*self.qs.shape[1], 0))
        plt.xlabel('x [km]', fontsize=12, fontweight='bold')
        plt.ylabel('y [km]', fontsize=12, fontweight='bold')
      plt.colorbar()
  
      plt.subplot(212)
      plt.title('Deflection [m]')
      if self.latlon:
        plt.imshow(self.w, extent=(0, self.dx*self.w.shape[0], self.dy*self.w.shape[1], 0))
        plt.xlabel('longitude [deg E]', fontsize=12, fontweight='bold')
        plt.ylabel('latitude [deg N]', fontsize=12, fontweight='bold')
      else:
        plt.imshow(self.w, extent=(0, self.dx/1000.*self.w.shape[0], self.dy/1000.*self.w.shape[1], 0))
        plt.xlabel('x [km]', fontsize=12, fontweight='bold')
        plt.ylabel('y [km]', fontsize=12, fontweight='bold')
      plt.colorbar()
docs: Plot multiple subplot figure for 2D array