As far as I see, if we already found the extremum points, then we just need to evaluate the function at the point itself and a few points around it to guess/estimate the behavior. It is in my opinion the most simple way for functions with explicitly known expression.
Of course, there's a downside: if the function is badly behaved (oscillations, singularities, etc.) this method may not give anything reliable, no matter how many points we use.
For a more reliable way see @phaedo's answer.
But the OP's function is nice enough, so it makes sense to try this simple method.
Also, the saddlepoints can give a bit of trouble, depending on where we check.
Let us pick the point from the OP $(-2,0)$ and check the function values around it:
$$f(-2,0)=\sqrt{56 \cdot 4+16 \cdot 2-31}+1+8 \cdot 2=32$$
$$f(-2.1,0)=33.59747 \ldots$$
$$f(-1.9,0)=30.39718 \ldots$$
$$f(-2,0.1)=f(-2,-0.1)=31.99767 \ldots$$
This is neither a maximum nor a minimun, so I would say, saddlepoint. Note that the function is even in $y$, so checking around $y=0$ is simple.
Now for $(16/7,0)$:
$$f(16/7,0)=-2.285714 \ldots$$
$$f(17/7,0)=-2.290772 \ldots$$
$$f(15/7,0)=-2.291607 \ldots$$
$$f(16/7,1/7)=f(16/7,-1/7)=-2.296764 \ldots$$
This point looks like a maximum to me, as all the values around it are smaller (note the negative signs).
Of course, as I said, this method may lead to uncertainties for less well behaved functions, and numerical errors may also play a role, so if you want to be certain, you have to use the Hessian matrix of second derivatives.
Using more than $5$ points (for example, $9$ points) should make the method more reliable.
Here, I made a little R program which uses the nine point method, where you can input a two variable function, a point you want to check, and the distance around it where we check the values. It wouldn't know a saddle point though, it only recognizes a maximum or a minimum:
PointCheck <- function(x0,y0,d,f){
f0 <- f(x0,y0);
f1 <- f(x0-d,y0);
f2 <- f(x0-d,y0-d);
f3 <- f(x0,y0-d);
f4 <- f(x0+d,y0-d);
f5 <- f(x0+d,y0);
f6 <- f(x0+d,y0+d);
f7 <- f(x0,y0+d);
f8 <- f(x0-d,y0+d);
c1 <- f0 > f1;
c2 <- f0 > f2;
c3 <- f0 > f3;
c4 <- f0 > f4;
c5 <- f0 > f5;
c6 <- f0 > f6;
c7 <- f0 > f7;
c8 <- f0 > f8;
if(c1+c2+c3+c4+c5+c6+c7+c8 == 8){return(paste(c("(",x0,y0,") is a maximum"), collapse=" "))}
if(c1+c2+c3+c4+c5+c6+c7+c8 == 0){return(paste(c("(",x0,y0,") is a minimum"), collapse=" "))}
if(c1+c2+c3+c4+c5+c6+c7+c8 > 0 && c1+c2+c3+c4+c5+c6+c7+c8 < 8){return(paste(c("(",x0,y0,") is neither a maximum, nor a minimum"), collapse=" "))}}
The output for $(-2,0)$:
x0 <- -2;
y0 <- 0;
d <- 0.05;
f <- function(x,y){sqrt(56*x^2-7*y^2-16*x-31)-8*x+1}
PointCheck(x0,y0,d,f)
[1] "( -2 0 ) is neither a maximum, nor a minimum"
The output for $(16/7,0)$:
x0 <- 16/7;
y0 <- 0;
d <- 0.05;
f <- function(x,y){sqrt(56*x^2-7*y^2-16*x-31)-8*x+1}
PointCheck(x0,y0,d,f)
[1] "( 2.28571428571429 0 ) is a maximum"
Note, that $d$ is also a good estimate for the error of this method, so if we try to check a point closer than $d$ to our critical point, the program wouldn't recognize the difference.
In fact, this function is a little tricky, because numerically it's hard to see if the second point is a maximum or not, because the plot for the whole range looks like this:

Now, plotting only the right side for $y=0$ gives us a more clear picture of a local maximum:
