After writing this, I noticed that Donu Arapura made essentially the same point in a comment. But perhaps my added detail will be of use.

There are two missing ingredients in what you wrote.

First, when you identify 1-forms with vector fields, you have implicitly chosen a bundle isomorphism $T\mathbb{R}^n \cong T^*\mathbb{R}^n$. Such a choice is the same thing as a choice of a Riemannian metric.

Second, when you identify n-forms on $\mathbb{R}^n$ with functions you are choosing an isomorphism between the exterior bundle $\bigwedge^n(\mathbb{R}^n)$ and a line bundle over $\mathbb{R}^n$. This is the same thing as a choice of orientation.

It is generally fine to leave the choice of standard Euclidean metric and standard orientation on $\mathbb{R}^n$ as implicit, but in this case these two choices are enough to determine a specific way to identify the rest of the De Rham complex. There are undoubtedly other ways to do it (particularly in low dimensions), but the way I will explain generalizes to any dimension (indeed, any oriented Riemannian manifold), it gives the right answers on $\mathbb{R}^3$, and it is compatible with Stokes' theorem.

Here's how it works. The Riemannian metric and the orientation determine a volume form $vol_n$ on $\mathbb{R}^n$; if $x_1, \ldots, x_n$ is a global oriented orthonormal coordinate system then we can take $vol_n = dx_1 \wedge \ldots \wedge dx_n$. There is an isomorphism $\star: \bigwedge^k(\mathbb{R}^n) \to \bigwedge^{n-k}(\mathbb{R}^n)$ (called the Hodge star operator) which is determined by the equation $\alpha \wedge \star(\alpha) = vol_n$. We will use this isomorphism to reduce high degree forms to low degree forms.

On $\mathbb{R}^2$ it's all very simple. Using the standard (oriented) coordinates $(x,y)$, the volume form is $dx \wedge dy$. The Euclidean metric identifies a vector field $(f, g)$ with the 1-form $f dx + g dy$. We have $d(f dx + g dy) = (f_x dx + f_y dy) \wedge dx + (g_x dx + g_y dy) \wedge dy = (g_x - f_y) dx \wedge dy$. The Hodge star operator identifies this 2-form with the 0-form $g_x - f_y$. Note that this is indeed the correct expression from the point of view of Stokes' theorem on the plane (AKA Green's theorem).

After writing this, I noticed that Donu Arapura made essentially the same point in a comment. But perhaps my added detail will be of use.

There are two missing ingredients in what you wrote.

First, when you identify 1-forms with vector fields, you have implicitly chosen a bundle isomorphism $T\mathbb{R}^n \cong T^*\mathbb{R}^n$. Such a choice is the same thing as a choice of a Riemannian metric.

Second, when you identify n-forms on $\mathbb{R}^n$ with functions you are choosing an isomorphism between the exterior bundle $\bigwedge^n(\mathbb{R}^n)$ and a line bundle over $\mathbb{R}^n$. This is the same thing as a choice of orientation.

It is generally fine to leave the choice of standard Euclidean metric and standard orientation on $\mathbb{R}^n$ as implicit, but in this case these two choices are enough to determine a specific way to identify the rest of the De Rham complex. There are undoubtedly other ways to do it (particularly in low dimensions), but the way I will explain generalizes to any dimension (indeed, any oriented Riemannian manifold), it gives the right answers on $\mathbb{R}^3$, and it is compatible with Stokes' theorem.

Here's how it works. The Riemannian metric and the orientation determine a volume form $vol_n$ on $\mathbb{R}^n$; if $x_1, \ldots, x_n$ is a global oriented orthonormal coordinate system then we can take $vol_n = dx_1 \wedge \ldots \wedge dx_n$. There is an isomorphism $\star: \bigwedge^k(\mathbb{R}^n) \to \bigwedge^{n-k}(\mathbb{R}^n)$ (called the Hodge star operator) which is determined by the equation $\alpha \wedge \star(\alpha) = vol_n$. We will use this isomorphism to reduce high degree forms to low degree forms.

On $\mathbb{R}^2$ it's all very simple. Using the standard (oriented) coordinates $(x,y)$, the volume form is $dx \wedge dy$. The Euclidean metric identifies a vector field $(f, g)$ with the 1-form $f dx + g dy$. We have $d(f dx + g dy) = (f_x dx + f_y dy) \wedge dx + (g_x dx + g_y dy) \wedge dy = (g_x - f_y) dx \wedge dy$. The Hodge star operator identifies this 2-form with the 0-form $g_x - f_y$. Note that this is indeed the correct expression from the point of view of Stokes' theorem on the plane (AKA Green's theorem).

ADDED One more comment. To give a geometric interpretation of your other choice of map from vector fields to functions, note that a choice of Riemannian metric on $\mathbb{R}^n$ yields a metric on the exterior algebra bundle of $\mathbb{R}^n$ and hence an inner product on the space of forms. If $V$ is a vector field and $\alpha$ is the corresponding 1-form, then the divergence of $V$ is nothing more than the function $d^\ast \alpha$ where $d^\ast$ is the adjoint of $d$ relative to this inner product. Up to a sign which I do not remember this adjoint is given by $d^\ast = \star d \star$, and this explains why the divergence on $\mathbb{R}^3$ can be regarded either as $d$ acting on 2-forms or as $d^\ast$ acting on 1-forms.