We report a $5sigma$ detection of the pairwise kinematic Sunyaev-Zeldovich (kSZ) effect, combining galaxy clusters in DESI imaging surveys and the Planck temperature maps. The detection is facilitated by both improvements in the data and in the analysis method. For the data, we adopt the recently released galaxy group catalog (Y20: cite{yang2020extended}) with $sim 10^6$ robustly-identified groups, and construct various galaxy cluster samples for the kSZ measurement. The Y20 catalogue also provides estimation of halo mass, which further improves the kSZ measurement by $sim 10%$. For the analysis method, we derive an optimal estimator of pairwise kSZ through the maximum likelihood analysis. It also handles potential systematic errors self-consistently. The baseline cluster sample, containing the $1.2times 10^5$ richest galaxy clusters of typical mass ~$ 10^{14} M_{odot}/h$ at typical redshift $0.2$-$0.5$, rules out the null hypothesis at $5sigma$. When fitting with a pairwise kSZ template from simulations, the signal is detected at $4.7sigma$ and the average optical depth is constrained as $bar{tau}_e=(1.66pm 0.35)times 10^{-4}$. We perform various internal checks, with different cluster selection criteria, different sky coverage and redshift range, different CMB maps, different filter sizes, different treatments of potential systematics and the covariance matrix. The kSZ effect is consistently detected with $2.5leq $S/N$leq 5.6$ and acceptable $chi^2_{rm min}$, across a variety of cluster samples. The S/N is limited by both the Planck resolution and the photo-z accuracy, and therefore can be significant improved with DESI spectroscopic redshift information and with other CMB experiments.