corteva / rioxarray

geospatial xarray extension powered by rasterio
https://corteva.github.io/rioxarray
Other
504 stars 80 forks source link

Migrate xarray examples #676

Open dcherian opened 1 year ago

dcherian commented 1 year ago

From the [visualization gallery](https://docs.xarray.dev/en/stable/examples/visualization_gallery.html#imshow()-and-rasterio-map-projections)

The data is here


imshow() and rasterio map projections

Using rasterio's projection information for more accurate plots.

This example extends recipes.rasterio and plots the image in the original map projection instead of relying on pcolormesh and a map transformation.

# opening raster files requires the rioxarray package
da = xr.open_rasterio("RGB.byte")

# The data is in UTM projection. We have to set it manually until
# https://github.com/SciTools/cartopy/issues/813 is implemented
crs = ccrs.UTM("18")

# Plot on a map
ax = plt.subplot(projection=crs)
da.plot.imshow(ax=ax, rgb="band", transform=crs)
ax.coastlines("10m", color="r")

Parsing rasterio 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 (see recipes.rasterio_rgb).

from pyproj import Transformer
import numpy as np

# opening raster files requires the rioxarray package
da = xr.tutorial.open_rasterio("RGB.byte")

x, y = np.meshgrid(da["x"], da["y"])
transformer = Transformer.from_crs(da.crs, "EPSG:4326", always_xy=True)
lon, lat = transformer.transform(x, y)
da.coords["lon"] = (("y", "x"), lon)
da.coords["lat"] = (("y", "x"), lat)

# Compute a greyscale out of the rgb image
greyscale = da.mean(dim="band")

# Plot on a map
ax = plt.subplot(projection=ccrs.PlateCarree())
greyscale.plot(
    ax=ax,
    x="lon",
    y="lat",
    transform=ccrs.PlateCarree(),
    cmap="Greys_r",
    shading="auto",
    add_colorbar=False,
)
ax.coastlines("10m", color="r")