Skip to content

Fix rasterio example in docs #1881

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Feb 2, 2018
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 9 additions & 1 deletion doc/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@
import datetime
import importlib

allowed_failures = []

print("python exec:", sys.executable)
print("sys.path:", sys.path)
for name in ('numpy scipy pandas matplotlib dask IPython seaborn '
Expand All @@ -32,6 +34,11 @@
print("%s: %s, %s" % (name, module.__version__, fname))
except ImportError:
print("no %s" % name)
if name == 'rasterio':
# not having rasterio should not break the build process
allowed_failures = ['gallery/plot_rasterio_rgb.py',
'gallery/plot_rasterio.py'
]

import xarray
print("xarray: %s, %s" % (xarray.__version__, xarray.__file__))
Expand Down Expand Up @@ -62,7 +69,8 @@

sphinx_gallery_conf = {'examples_dirs': 'gallery',
'gallery_dirs': 'auto_gallery',
'backreferences_dir': False
'backreferences_dir': False,
'expected_failing_examples': allowed_failures
}

autosummary_generate = True
Expand Down
9 changes: 6 additions & 3 deletions doc/gallery/plot_rasterio.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
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.
original map projection (see :ref:`recipes.rasterio_rgb`).
"""

import os
Expand Down Expand Up @@ -44,10 +44,13 @@
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())
da.plot.imshow(ax=ax, x='lon', y='lat', rgb='band',
transform=ccrs.PlateCarree())
greyscale.plot(ax=ax, x='lon', y='lat', transform=ccrs.PlateCarree(),
cmap='Greys_r', add_colorbar=False)
ax.coastlines('10m', color='r')
plt.show()

Expand Down
44 changes: 44 additions & 0 deletions doc/gallery/plot_rasterio_rgb.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
# -*- coding: utf-8 -*-
"""
.. _recipes.rasterio_rgb:

============================
imshow() and map projections
============================

Using rasterio's projection information for more accurate plots.

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

import os
import urllib.request
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

# 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')

# Normalize the image
da = da / 255
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this is strictly necessary with RGB support in imshow, though I suppose it doesn't hurt.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's unfortunately necessary, see #1880


# The data is in UTM projection
# Until https://github.com/SciTools/cartopy/issues/813 is implemented
# we have to do this manually
crs = ccrs.UTM('18N')

# Plot on a map
ax = plt.subplot(projection=crs)
da.plot.imshow(ax=ax, transform=crs)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For pedagogical purposes, it might be nice to explicitly write rga='band'.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, will do

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

# Delete the file
os.remove('RGB.byte.tif')