Non Equispaced / Non Uniform DFT Bandwidth I need to construct Fourier transform of non-equispaced data. 

That is, I have signal $s(t)$, $t\in[0,T]$ sampled at non-equispaced points $t_k$, $k=0...N-1$ with sample values $s_k = s(t_k)$. For Fourier transform I use approximation of integral:
$$
   \hat S(\omega) = \int\limits_{-\infty}^\infty s(t)e^{-i\omega t}dt \approx \hat S_d(\omega) = \sum_{k=0}^{N-1}s_ke^{-i\omega t_k}\Delta t_k
\tag{1}
$$
where $\Delta t_k = (t_{k+1}-t_{k-1})/2$. As sampling points in frequency domain I choose $\omega_n = \frac{2\pi}{T}$.

My question is: as I can evaluate (1) for any $\omega$, what is the maximum $\omega$ such $\hat S_d(\omega)$ "adequately" represents $\hat S(\omega)$? What is the maximum of $n$ for $\omega_n$ I can use? For DFT we got Nyquist frequency. Do we have something similar for NDFT? Any references would be appreciated.

Note: I'm aware of such things as NDFT and NFFT. Though formula for NDFT, as presented in most papers, is 
$$
\hat S_d(\omega_n) = \sum_{k=0}^{N-1}s_ke^{-i\omega_n t_k}
$$
I strongly believe that I need to use formula (1) as I'm trying to build a periodogram. 

And I'm not interested in fast ways of computing NDFT yet, so I'm not considering NFFT.
Your formula isn't accurate.  Since you aren't trying for speed and internally  consider the interpolation between data points as a step/pulse then the formula should be.

$$
S(\omega)=\sum_{n=0}^{N-1}s_{k}sinc(\frac{\omega\triangle\left(t_{k}\right)}{2})e^{-wt_{k}}
$$ 
Having said that, what is "adequate"?  This is, more or less, subjective.  The original Nyquest criterion is based upon the result that any (freq &lt; sam/2) (I think there is some small quibbling on that) can be perfectly reconstructed; if it and the sampling lasts long enough.  It says nothing about out of band signals; but the actual calculations will show aliasing.
It behooves you to generate the signals and interference that you envision and run it through your proposed algorithm. An alternative is to treat your system as a filter and try to "deconvolve" it. This is hazardous but can be done.
For a in depth look (with the hazards) read "Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements"  SIAM.
Now MRI/CAT scanners use various means of reconstructed signals.  All of the mathematicians seem to think that their technical success is an example of extreme luck; although I am sure that the designers might take exception to that.  You might try borrowing some of their mathematical techniques though.
Now in general instrumentation what we concentrate on is "noise temperature" and such.  This can be translated as: how far will your system signal/noise (or signal/distortion) be from the optimal number.  If you can get within 1db of the optimal S/N value then you should probably give up.  Of course this implies that you can calculate the optimal S/N.  In most instruments this can be done and this criterion gives you a stopping point; a point to give up and look for something else to do :)