Weibull plot to assess goodness of fit I would like show whether my data fit a Weibull distribution.
What kinds of graph should I plot?
Is the following graph appropriate one and what is this called?
How can I plot it in R?


The plot you display is sometimes called a Weibull plot.
See, for example, https://www.itl.nist.gov/div898/handbook/eda/section3/weibplot.htm
(which shows your exact plot)
and also https://en.wikipedia.org/wiki/Weibull_distribution#Weibull_plot
While it's not the only sort of plot you might want to do*, it's a reasonably common way of looking at goodness of fit for a Weibull distribution, and performs reasonably well. If a Weibull model is suitable, it should look close to straight. You should be aware that with Weibull data the lower left end may be considerably more "wiggly" than the upper right end; its behaviour is not symmetric. The log-data values have the long left tail of a flipped Gumbel distribution.
To my mind it makes more sense to interchange the axes from the way around your plot has it, so the random quantity is on the y-axis, just as is usually done with a normal Q-Q plot, and as is common practice with scatterplots more generally. I will do this in my answer, but it's easy to swap them back if you're determined to have the logged data on the x-axis.
It's important to generate Weibull data across a variety of shape parameter values and sample sizes to see what Weibull data can look like, as well as to generate some non-Weibull data (e.g. gamma, lognormal, chi). It's useful to see what the plot tends to look like at a few different sample sizes as well.
These plots are easy to do in R, following the information at the above two links.
Here's some sample code. This is just a quick outline of how to do it, not what I'd suggest for a proper function:
    x <- rweibull(67,4.1,5.25) # just some random Weibull data 
                               # to show an example

    p <- ppoints(x, .3)
    q <- log(-log(1-p))
    o <- order(x)
    plot(q,log(x[o]))
    coeffs <- lm(log(x[o])~q)$coefficients
    ests <- c(1/coeffs[2],exp(coeffs[1]),quantile(x,1-exp(-1)))
    names(ests) <- c("shape.line","scale.line","scale.quantile")
    ests

This gives a plot something like this one:

You could interchange the axes if you felt the need.
When working by hand, it would be easier to use a line joining two quantiles rather than a linear regression fit. Indeed, because of the tendency of the lower-left end to be noisy, it may be worth the slight robustness this brings. With actual Weibull data, it looks like taking the slope of the line joining the points for the  $0.167$ and $0.975$ quantiles may work well, but making the intercept of the drawn line correspond to the quantile-estimate of the scale parameter, at the $(1-1/e)$ quantile. It may be worth giving up a little efficiency to gain more robustness by moving the slope quantiles in somewhat further. Another alternative would be to use the quartiles, as is done for the normal Q-Q plot in R.
If you need the axes scaled in a similar way to the plot you show, you could replace the 'plot' command above with something like:
    qax=c(0.1,0.5,1,5,10,20,50,90,99,99.9)
    plot(q,x[o], ylab="x", xaxt="n", log="y", 
           xlim=range(log(-log(1-qax/100))), ylim=range(x))
    axis(1, at=log(-log(1-qax/100)), labels=qax)


but it might be worth pulling the range on the theoretical-quantiles axis (what your plot calls "Weibull probability") in to just fit the data, which you can work out (e.g. from the number of data-values).
There are numerous posts relating to Weibull distributions, Weibull plots and other ways of looking at Weibull goodness of fit on site. e.g. try searches, such as a search for Weibull plot.

*(e.g. you might do a Q-Q plot or a P-P plot using fitted parameters, or an ecdf with fitted Weibull cdf)
