GH5005 fix documentation on open_rasterio (pydata#5021)

* GH5005 fix documentation on open_rasterio

* Example gets the transformation from attrs instead of da.transform

* Doctest example in open_rasterio

* Doctest for open_rasterio, reverting commit

* Applying blackdoc

* Remove unnecessary imports

Co-authored-by: keewis <keewis@users.noreply.github.com>

* Getting transform from attribute

Co-authored-by: keewis <keewis@users.noreply.github.com>

Co-authored-by: keewis <keewis@users.noreply.github.com>

xarray/backends/rasterio_.py

@@ -174,11 +174,46 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, loc
 
     You can generate 2D coordinates from the file's attributes with::
 
-        from affine import Affine
-        da = xr.open_rasterio('path_to_file.tif')
-        transform = Affine.from_gdal(*da.attrs['transform'])
-        nx, ny = da.sizes['x'], da.sizes['y']
-        x, y = np.meshgrid(np.arange(nx)+0.5, np.arange(ny)+0.5) * transform
+        >>> from affine import Affine
+        >>> da = xr.open_rasterio(
+        ...     "https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif"
+        ... )
+        >>> da
+        <xarray.DataArray (band: 3, y: 718, x: 791)>
+        [1703814 values with dtype=uint8]
+        Coordinates:
+          * band     (band) int64 1 2 3
+          * y        (y) float64 2.827e+06 2.826e+06 2.826e+06 ... 2.612e+06 2.612e+06
+          * x        (x) float64 1.021e+05 1.024e+05 1.027e+05 ... 3.389e+05 3.392e+05
+        Attributes:
+            transform:      (300.0379266750948, 0.0, 101985.0, 0.0, -300.041782729805...
+            crs:            +init=epsg:32618
+            res:            (300.0379266750948, 300.041782729805)
+            is_tiled:       0
+            nodatavals:     (0.0, 0.0, 0.0)
+            scales:         (1.0, 1.0, 1.0)
+            offsets:        (0.0, 0.0, 0.0)
+            AREA_OR_POINT:  Area
+        >>> transform = Affine(*da.attrs["transform"])
+        >>> transform
+        Affine(300.0379266750948, 0.0, 101985.0,
+               0.0, -300.041782729805, 2826915.0)
+        >>> nx, ny = da.sizes["x"], da.sizes["y"]
+        >>> x, y = transform * np.meshgrid(np.arange(nx) + 0.5, np.arange(ny) + 0.5)
+        >>> x
+        array([[102135.01896334, 102435.05689001, 102735.09481669, ...,
+                338564.90518331, 338864.94310999, 339164.98103666],
+               [102135.01896334, 102435.05689001, 102735.09481669, ...,
+                338564.90518331, 338864.94310999, 339164.98103666],
+               [102135.01896334, 102435.05689001, 102735.09481669, ...,
+                338564.90518331, 338864.94310999, 339164.98103666],
+               ...,
+               [102135.01896334, 102435.05689001, 102735.09481669, ...,
+                338564.90518331, 338864.94310999, 339164.98103666],
+               [102135.01896334, 102435.05689001, 102735.09481669, ...,
+                338564.90518331, 338864.94310999, 339164.98103666],
+               [102135.01896334, 102435.05689001, 102735.09481669, ...,
+                338564.90518331, 338864.94310999, 339164.98103666]])
 
     Parameters
     ----------