Reproject xarray Dataset with pyproj

Re-projecting is the process of changing the representation of locations from one coordinate system to another. All projections of locations on the Earth into a two-dimensional plane are distortions, the projection that is best for your application may be different from the projection associated with the data you import.

Sometimes, this is the case for xarray Datasets, which may come with an unwanted projection. However, we could easily reproject the dataset into another projections with the help of pyproj. For example, let's change the Coordinate Reference System (CRS) from a projection called epsg:32618 to WGS84, one of the most commonly used CRS in latitude-longitude.

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.

%matplotlib inline

import numpy as np
import xarray as xr
import cartopy.crs as ccrs
from pyproj import Transformer
import matplotlib.pyplot as plt

Read a projected data (epsg:32618)

url = 'https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif'
da = xr.open_rasterio(url)
da.crs

[out]: '+init=epsg:32618'

Transform projection to WGS84 (epsg: 4326)

%%time
xv, yv = np.meshgrid(da.x, da.y)

transformer = Transformer.from_crs("epsg:32618", 
                                   "epsg:4326",
                                   always_xy=True,
                                  ) 

lon, lat = transformer.transform(xv, yv)
da.coords['lon'] = (('y', 'x'), lon)
da.coords['lat'] = (('y', 'x'), lat)
da.attrs['crs']  = '+init=epsg:4326'

Visualize

greyscale0 = da.mean(dim='band')

# Plot on a map
ax = plt.subplot(projection=ccrs.PlateCarree())
greyscale0.plot(ax=ax, 
                x='lon', y='lat', 
                transform=ccrs.PlateCarree(),
                #cmap='Greys_r', 
                add_colorbar=True
               )

ax.coastlines('10m', color='r')

In many cases, transforming projections are carried out as a neccesary step before interpolating the reprojected dataset into an irregular projection coordinate in lon/lat.

References

Travis E, Oliphant. A guide to NumPy, USA: Trelgol Publishing, (2006).

Stéfan van der Walt, S. Chris Colbert and Ga?l Varoquaux. The NumPy Array: A Structure for Efficient Numerical Computation, Computing in Science & Engineering, 13, 22-30 (2011), DOI:10.1109/MCSE.2011.37

Fernando Pérez and Brian E. Granger. IPython: A System for Interactive Scientific Computing, Computing in Science & Engineering, 9, 21-29 (2007), DOI:10.1109/MCSE.2007.53

John D. Hunter. Matplotlib: A 2D Graphics Environment, Computing in Science & Engineering, 9, 90-95 (2007), DOI:10.1109/MCSE.2007.55

https://xarray.pydata.org/en/stable/

https://github.com/pyproj4/pyproj