Approximation of smooth diffeomorphisms by polynomial diffeomorphisms?

The answer is 'no', because polynomial mappings with polynomial inverses preserve volumes up to a constant multiple.

To see why this property holds, suppose that $p:\mathbb{R}^d\to\mathbb{R}^d$ is a polynomial mapping with polynomial inverse $q:\mathbb{R}^d\to\mathbb{R}^d$. Then $p$ and $q$ extend to $\mathbb{C}^d$ as polynomial maps with polynomial inverses. This means that the Jacobian determinant of $p$ on $\mathbb{C}^d$ is a complex polynomial with no zeros and hence must be a (nonzero) constant.

Now, consider a diffeomorphism $f:\mathbb{R}^d\to\mathbb{R}^d$ that is radial, i.e., $f(x) = m(|x|^2)x$ for some smooth function $m>0$. One can easily choose $m$ in such a way that $m(4)=1/2$ and $m(9)=4/3$, so that $f$ maps the ball of radius $2$ about the origin diffeomorphically onto the ball of radius $1$ about the origin while it maps the ball of radius $3$ about the origin diffeomorphically onto the ball of radius $4$ about the origin.

Let $\epsilon>0$ be very small and suppose that $\|f-p\|_{\infty;U} <\epsilon$ for $U$ chosen to be some very large ball centered on the origin. Then $p$ maps the sphere of radius $2$ about the origin to within an $\epsilon$-neighborhood of the sphere of radius $1$, while it maps the sphere of radius $3$ about the origin to within an $\epsilon$-neighborhood of the sphere of radius $4$. It's easy to see from this that $p$ cannot have constant Jacobian determinant.

Added remark: The group $\mathrm{SDiff}(\mathbb{R}^d)$ consisting of volume-preserving diffeomorphisms of $\mathbb{R}^d$ is a 'Lie group' in Sophus Lie's original sense (i.e., a group of diffeomorphisms defined by the satisfaction of a system of differential equations; in this case, that the Jacobian determinant be equal to $1$).

The subgroup $\mathcal{SP}(\mathbb{R}^d)\subset \mathrm{SDiff}(\mathbb{R}^d)$ consisting of volume-preserving polynomial diffeomorphisms with polynomial inverses however, is not a 'Lie subgroup' in Lie's original sense when $d>1$, as it cannot be defined by the satisfaction of a system of differential equations: It contains all of the mappings of the form $p(x) = x + a\,(b{\cdot}x)^m$ where $a,b\in\mathbb{R}^d$ satisfy $a\cdot b = 0$ and $m>1$ is an integer (indeed, $p^{-1}(y) = y - a\,(b{\cdot}y)^m$), plus, it contains $\mathrm{SL}(d,\mathbb{R})$ and the subgroup consisting of the translations. Using this, it is easy to show that, for any $f\in\mathrm{SDiff}(\mathbb{R}^d)$ and for any integer $k$, there exists a $p\in \mathcal{SP}(\mathbb{R}^d)$ such that $f$ and $p$ have the same Taylor series at the origin up to and including order $k$. Thus, $\mathcal{SP}(\mathbb{R}^d)$ cannot be defined by a system of differential equations (in Lie's sense).

Using this Taylor approximation property, one can prove that $\mathcal{SP}(\mathbb{R}^d)$, like $\mathrm{SDiff}(\mathbb{R}^d)$, acts transitively on $n$-tuples of distinct points in $\mathbb{R}^d$ for any integer $n$. Whether one can prove that $\mathcal{SP}(\mathbb{R}^d)$ can 'uniformly approximate' $\mathrm{SDiff}(\mathbb{R}^d)$ on compact sets is an interesting question.

An illustration for one of the examples in the answer by Robert Bryant. It is supposed to convey the feeling of something extremely rigid, unyielding and inflexible.

Image of the square $[-1,1]\times[-1,1]$ under the map $(x,y)\mapsto(x-y^2-2x^2y-x^4,y+x^2)$ (composite of $(x,y)\mapsto(x-y^2,y)$ with $(x,y)\mapsto(x,y+x^2)$).