text: # -*- coding: utf-8 -*-
  """
  .. _recipes.rasterio:
  
  =================================
  Parsing rasterio's geocoordinates
  =================================
  
  
  Converting a projection's cartesian coordinates into 2D longitudes and
  latitudes.
  
  These new coordinates might be handy for plotting and indexing, but it should
  be kept in mind that a grid which is regular in projection coordinates will
  likely be irregular in lon/lat. It is often recommended to work in the data's
  original map projection.
  """
  
  import os
  import urllib.request
  import numpy as np
  import xarray as xr
  import cartopy.crs as ccrs
  import matplotlib.pyplot as plt
  from rasterio.warp import transform
  
  
  # Download the file from rasterio's repository
  url = 'https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif'
  urllib.request.urlretrieve(url, 'RGB.byte.tif')
  
  # Read the data
  da = xr.open_rasterio('RGB.byte.tif')
  
  # Compute the lon/lat coordinates with rasterio.warp.transform
  ny, nx = len(da['y']), len(da['x'])
  x, y = np.meshgrid(da['x'], da['y'])
  
  # Rasterio works with 1D arrays
  lon, lat = transform(da.crs, {'init': 'EPSG:4326'},
                       x.flatten(), y.flatten())
  lon = np.asarray(lon).reshape((ny, nx))
  lat = np.asarray(lat).reshape((ny, nx))
  da.coords['lon'] = (('y', 'x'), lon)
  da.coords['lat'] = (('y', 'x'), lat)
  
  # Plot on a map
  ax = plt.subplot(projection=ccrs.PlateCarree())
  da.plot.imshow(ax=ax, x='lon', y='lat', rgb='band',
                 transform=ccrs.PlateCarree())
  ax.coastlines('10m', color='r')
  plt.show()
  
  # Delete the file
  os.remove('RGB.byte.tif')
  