This paper gives a classification of first order polynomial differential operators of form $\mathscr{X} = X_1(x_1,x_2)\delta_1 + X_2(x_1,x_2)\delta_2$, $(\delta_i = \partial/\partial x_i)$. The classification is given through the order of an operator that is defined in this paper. Let $X=\mathscr{X}y$ to be the differential polynomial associated with $\mathscr{X}$, the order of $\mathscr{X}$, $\mathrm{ord}(\mathscr{X})$, is defined as the order of a differential ideal $\Lambda$ of differential polynomials that is a nontrivial expansion of the ideal $\{X\}$ and with the lowest order. In this paper, we prove that there are only four possible values for the order of a differential operator, $0$, $1$, $2$, $3$, or $\infty$. Furthermore, when the order is finite, the expansion $\Lambda$ is generated by $X$ and a differential polynomial $A$, which can be obtained through a rational solution of a partial differential equation that is given explicitly in this paper. When the order is infinite, the expansion $\Lambda$ is just the unit ideal. In additional, if, and only if, the order of $\mathscr{X}$ is $0$, $1$, or $2$, the polynomial differential equation associating with $\mathscr{X}$ has Liouvillian first integrals. Examples and connections with Godbilon-Vey sequences are also discussed.