Is the correlation coefficient useful for measuring linearity with only 10 observations? I need to measure the linearity between two variables but I only have 10 samples.
Is it relevant to use correlation in this case?
If not, what else could I use?
 A: Correlation does not measure linearity--full stop.
I will provide a brief explanation and then a pair of examples showing two similar-looking datasets of ten observations each with identical (and high) correlation coefficients--but one of which exhibits a linear relationship and the other one shows significant evidence of nonlinearity.

Correlation can be interpreted only by assuming an underlying linear relationship.  If you're inclined to object to this pronouncement, please contemplate the Anscombe quartet before proceeding.  This collection of bivariate datasets (of only 11 points each--quite relevant here) shows how strikingly different relations, some of them obviously nonlinear, can correspond to the very same value of the correlation coefficient.
A conventional way to assess the linearity of a relation among bivariate data $\mathcal{X}=(x_i,y_i),\ i=1,2,\ldots,n$ is to compare two models.  The base model posits a linear relation.  Because these are data, we have to prepare for the possibility of departures from this relation.  Such departures are the "errors" in the responses $\varepsilon_i,$ whence the model is of the form
$$y_i = \beta_0 + \beta_1 x_i + \varepsilon_i\tag{1}$$
for unknown parameters $\beta_0$ and $\beta_1.$  Fitting this model is called "regression."  It results in estimates for the parameters, correspondingly named $\hat\beta_0$ and $\hat\beta_1.$  They are functions of the dataset $\mathcal X.$
At this point you have to contemplate what "nonlinearity" means to you.  In some applications this might take on a specific form.  For instance, the alternative model might posit a sinusoidal pattern in the errors.  The simplest, most general alternative introduces one more parameter in the form
$$y_i = \gamma_0 + \gamma_1 x_i + \gamma_2 f(x_i) + \delta_i\tag{2}$$
where the $\delta_i$ represent the errors and $f$ is some function of the explanatory variable chosen beforehand to represent likely or meaningful departures from linearity.  When this is a "better" fit to the data in some meaningful sense, you should decide to use it: that is, you will act as if the relationship between the $x_i$ and $y_i$ is not linear.
A conventional way to carry out this program uses Ordinary Least Squares (OLS) regression and chooses $f(x)=x^2:$ a quadratic term.  Assuming OLS is suitable for your objectives and your data, this technology offers a simple hypothesis test for comparing the two models, variously known as a t test (for the quadratic coefficient), an added-variable $F$ test, or an Analysis of Variance.  They all produce equivalent information.

Example 1
Let me illustrate with randomly generated data.  Here is a plot of ten $(x_i,y_i)$ points.  On it I have superimposed the OLS fit for model $(1);$ that is, the graph of the function $x\to \hat\beta_0 + \hat\beta_1 x^2.$

The correlation coefficient is $0.873:$ pretty high.  But does it, or does it not, indicate a linear relation?  Let's fit model $(2).$  Here is a summary of the standard OLS output.
Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  1.08750    0.78832   1.380   0.2102  
x           -0.30873    0.32923  -0.938   0.3796  
I(x^2)       0.07142    0.02917   2.449   0.0442 *

The t-test for the quadratic term yields a p-value of $0.0442$ (shown at bottom right).  If you take $5\%$ to be your level of significance, then because $0.0442 \lt 5/100$ you would take this to be evidence of nonlinearity.
The ANOVA table compares the fits of models $(1)$ and $(2)$ directly.
Model 1: y ~ x
Model 2: y ~ x + I(x^2)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)  
1      8 5.8381                             
2      7 3.1446  1    2.6934 5.9956 0.0442 *

It does this by comparing the "residual sums of squares" RSS in the two models.  Note that its p-value is the same as before: it's in this sense that this test is equivalent to the preceding one.
Here is a graph of the quadratic fit in red:

The preceding tests have established that the red curve is a sufficiently good improvement over the black line, as a description of the overall $(x_i,y_i)$ relationship, that it may be the better model.
Example 2
Here is a plot of another dataset.  The $x_i$ values are the same but the $y_i$ values have changed.  Nevertheless, the correlation coefficient is exactly the same as before.

Here is the quadratic fit:

The ANOVA produces a p-value of only $0.3382:$ that wouldn't generally be considered sufficient evidence to pick the quadratic model over the linear one.  Subject to all the choices, judgments, and assumptions made to get to this point, it would be valid to conclude that these data reflect a linear relation.

For those who would like to explore these data further, here is the R code that produced them:
x <- 1:10
set.seed(17)
y <- 0.2 * x + rnorm(length(x))
z <- y + 2.535153*(x/5-1)^2

The coefficient of 2.535153 was chosen to make the correlation coefficients (of x vs y and x vs z) identical.
