identify level shifts in a time series I have a time series as follows:

I want to identify the locations of level shifts in this time series. Are there R packages available to do the job?
 A: The kind of analysis for finding such shifts is called changepoint analysis (see other questions tagged like this). This can be achieved my using maximum likelihood estimation, where for $m$ changepoints the likelihood function is
$$ 
L(m, \tau_{1:m}, \theta_{1:(m+1)}) = \prod_{i=1}^{m+1} p(y_{(\tau_{i-1} + 1):\tau_i}\mid \theta_i)
$$
where $y_1,\dots,y_n$ is your data, $1 < \tau_1 <\dots<\tau_m<n$ are the boundary points marking the changes, and probability distributions $p$ are parametrized by $\theta_i$ for each $i$-th segment. 
In your case, you will be dealing with changes in the mean, but in the last segment there is a visible change in mean and variance.
As about software, there is the changepoint package for R. If you need to forecast given the series as well, there is the Prophet package for both R and Python, that automatically detects the changepoints and uses them as a part of the time-series model.

Rebecca Killick and Idris A. Eckley. (2013) changepoint: An R Package
  for Changepoint Analysis. (online paper)
Eckley, I.A., Fearnhead, P. and Killick, R. (2011) Analysis of
  changepoint models. [in:] Bayesian Time Series Models, eds. D.
  Barber, A.T. Cemgil and S. Chiappa, Cambridge University Press.

A: In R there are a wide range of methods/packages available for this task. A selection of them are summarized in CRAN's Task View on time series: https://cran.r-project.org/web/views/TimeSeries.html. Here is the section on change detection:

Change point detection is provided in strucchange and strucchangeRcpp
(using linear regression models) and in trend (using nonparametric
tests). The changepoint package provides many popular changepoint
methods, and ecp does nonparametric changepoint detection for
univariate and multivariate series. changepoint.np implements the
nonparametric PELT algorithm, changepoint.mv detects changepoints in
multivariate time series, while changepoint.geo implements the
high-dimensional changepoint detection method GeomCP. Factor-augmented
VAR (FAVAR) models are estimated by a Bayesian method with FAVAR.
InspectChangepoint uses sparse projection to estimate changepoints in
high-dimensional time series. Rbeast provides Bayesian change-point
detection and time series decomposition. breakfast includes methods
for fast multiple change-point detection and estimation.

