Yes. You do need to think about what you are doing.
I'll create a quasi-triangular region as an example.
px = [0, 1, .5, .49];
py = [0, 0, 1, 1-0.1];
plot(px([1 2 3 4 1]),py([1 2 3 4 1]),'r-')
So in the sequence I listed them, then connecting back to the beginning, the area is:
polyarea(px([1 2 3 4 1]),py([1 2 3 4 1]))
This is exactly as I would expect, just a bit less than 1/2, which would be the area of the triangle that forms the convex hull. However, in a different sequence, the included area is seen as:
polyarea(px([1 4 2 3 1]),py([1 4 2 3 1]))
plot(px([1 4 2 3 1]),py([1 4 2 3 1]),'g-')
So the same points, but in a different sequence, we got a completely different result, and that is expected.
Can you know which sequence was intended? Of course not. That is impossible, nor should a tool like polyarea be expected to know which you intended.
In the case of a convex polygon, IF the region is KNOWN to be convex, then you can just use a sort on polar angle.
t = rand(1,50)*2*pi;
px = cos(t);
py = sin(t);
As it turns out, I would have expected an area that was reasonably close to pi.
mux = mean(px);
muy = mean(py);
[th,r] = cart2pol(px,py);
[~,ind] = sort(th);
So the area is now seen as quite close to pi. Which is correct? As I said, IF WE KNEW the region was a convex polygon, then the latter is easy. But given only the points from some general possibly non-convex domain in unsorted order, the problem is essentially impossible to solve without knowing the desired sequence of the points.
Is the area that of the blue curve in the last plot, or the area "inside" the red curve? See that when you have a Catholic polygon, that is a polygon that crosses itself multiple time,s then some parts of the polygon have positive area, some have a negative area. This effectively has some of those areas cancel each other out. (Think right hand rule, for those who know what a cross product is and what it tells you.)