diff --git a/README.md b/README.md index ead391c..1b59684 100644 --- a/README.md +++ b/README.md @@ -25,9 +25,10 @@ manipulation, processing, and visualization. ## Example -Here's a quick example of loading this +Here's a quick example of loading and plotting this [Landsat 9 scene from the city of Manaus, Brazil](https://doi.org/10.6084/m9.figshare.24167235.v1), -which is where the Amazonas and Negro river merge: +which is where the Solimões (brown water) and Negro (black water) rivers merge +to form the Amazon river: ```python import xlandsat as xls @@ -51,7 +52,7 @@ plt.axis("scaled") plt.show() ``` -RGB image showing the city on the left and the black waters of the Negro river merging with the brown waters of the Amazonas river +RGB image showing the city on the left and the black waters of the Negro river merging with the brown waters of the Solimões river ## Project goals diff --git a/doc/composites.rst b/doc/composites.rst index df1cec7..3ea0779 100644 --- a/doc/composites.rst +++ b/doc/composites.rst @@ -1,24 +1,28 @@ .. _composites: -Composites -========== +Making composites +================= -Plotting individual bands is good but we usually want to make some composite -images to visualize information from multiple bands at once. -For that, we have to create **composites**. -We provide the :func:`xlandsat.composite` function to make this process easier. +Plotting individual bands is nice but we usually want to make some composite +images, both RGB and false-color, to visualize information from multiple bands +at once. +We provide the :func:`xlandsat.composite` function to make this process easier +and produce composites that are compatible with both :mod:`xarray` and +:mod:`matplotlib`. -As an example, let's load an example scene from the -`Brumadinho tailings dam disaster `__: +As an example, let's load our example scene from Manaus, Brazil, which is where +the Solimões (brown water) and Negro (black water) rivers merge to form the +Amazon river: .. jupyter-execute:: import xlandsat as xls import matplotlib.pyplot as plt - path = xls.datasets.fetch_brumadinho_after() + path = xls.datasets.fetch_manaus() scene = xls.load_scene(path) + scene RGB composites @@ -34,10 +38,12 @@ bands aren't specified. rgb = xls.composite(scene) rgb -The composite has a similar layout as the bands of a scene but with an extra -``"channel"`` dimension corresponding to red, green, blue, and -alpha/transparency. The values are scaled to the [0, 255] range and the -composite is an array of unsigned 8-bit integers. +The composite is also an :class:`xarray.DataArray` and is similar to the bands. +It has the same easting and northing dimensions but also an extra ``"channel"`` +dimension, which corresponds to red, green, blue, and alpha/transparency. This +extra dimension is what combines the three bands into a single color image. The +values are scaled to the [0, 255] range and the composite is an array of +unsigned 8-bit integers. .. admonition:: Transparency :class: note @@ -48,11 +54,21 @@ composite is an array of unsigned 8-bit integers. blue). -Plotting a composite --------------------- +Plotting composites +------------------- -Composites can be plotted using :meth:`xarray.DataArray.plot.imshow` (using -:meth:`~xarray.DataArray.plot` won't work and will display histograms instead). +Composites can be plotted using the :meth:`xarray.DataArray.plot.imshow` method: + +.. jupyter-execute:: + + rgb.plot.imshow() + + +With no arguments, xarray takes care of creating the new figure and adding a +lot of the different plot elements, like axis labels and units. +If we want more control over the plot, we must create a matplotlib figure and +axes separately and tell :meth:`~xarray.DataArray.plot.imshow` to plot on those +instead: .. jupyter-execute:: @@ -61,17 +77,23 @@ Composites can be plotted using :meth:`xarray.DataArray.plot.imshow` (using rgb.plot.imshow(ax=ax) # The "long_name" of the composite is the band combination - ax.set_title(f"Composite: {rgb.attrs['long_name']}") + ax.set_title(rgb.attrs["long_name"].title()) - # Make sure pixels are square and don't have any distortions from plotting + # Make sure pixels are square when plotting to avoid distortions ax.set_aspect("equal") plt.show() -Well, this looks bad because that bright cloud is making it so the ground -pixels have only a small share of the full range of available values. This can -be mitigated by rescaling the intensity of the image to a smaller range of -reflectance values. +Well, this looks bad because some very bright pixels in the city are making the +majority of the other pixels have only a small share of the full range of +available values. This can be mitigated by rescaling the intensity of the image +to a smaller range of reflectance values. + +.. note:: + + Using :meth:`xarray.DataArray.plot` instead of + :meth:`xarray.DataArray.plot.imshow` won't work and will display histograms + instead. Rescaling intensity (AKA contrast stretching) @@ -84,16 +106,21 @@ the following values by trial and error until it looks nice: .. jupyter-execute:: - rgb = xls.composite(scene, rescale_to=[0.03, 0.2]) + rgb = xls.composite(scene, rescale_to=[0.01, 0.2]) + # Pretty much the same plotting code fig, ax = plt.subplots(1, 1, figsize=(10, 6)) rgb.plot.imshow(ax=ax) - ax.set_title(f"Rescaled composite: {rgb.attrs['long_name']}") + ax.set_title(f"Rescaled {rgb.attrs['long_name'].title()}") ax.set_aspect("equal") plt.show() -Notice that we can more clearly see the colors of the ground but we lose a lot -of detail in the clouds. +Notice that we can more clearly see the different colors of the forest and the +rivers. +However, it can still be a bit hard to distinguish between some of the water +bodies and the forest in the right side of the scene. +Other band combinations can generate composites that better highlight these +features. .. note:: @@ -106,23 +133,48 @@ of detail in the clouds. Color infrared composites ------------------------- -Another common type of composite is the color infrared (CIR) composites. These -change the bands used to NIR, red, and green and serve primarily to distinguish -healthy vegetation from other objects in the scene. Let's make one by specifying -this band combination to :func:`xlandsat.composite` to see if we can more clearly spot the dam flood. +Another common type of composite is the color infrared (CIR) composites, +which uses the NIR, red, and green bands as the red, green, blue channels. +The added information of the NIR band helps highlight vegetation, which can +help us distinguish between the water and forest on the right. +Let's make one by specifying +this band combination to :func:`xlandsat.composite` to see what happens: .. jupyter-execute:: - cir = xls.composite(scene, bands=("nir", "red", "green"), rescale_to=[0, 0.4]) + cir = xls.composite( + scene, rescale_to=[0, 0.4], bands=("nir", "red", "green"), + ) fig, ax = plt.subplots(1, 1, figsize=(10, 6)) cir.plot.imshow(ax=ax) - ax.set_title(f"Composite: {rgb.attrs['long_name']}") + ax.set_title(cir.attrs["long_name"].title()) ax.set_aspect("equal") plt.show() -The flood region can be clearly spotted in the image above as the brown/gray -blog in the center. +In this composite, the contrast between the forest and water bodies on the +right are much clearer! + + +Other composites +---------------- + +You can make pretty much any composite you'd like by passing the correct band +combination to :func:`xlandsat.composite`. +For example, let's make one with NIR as red, SWIR1 as green, and red as blue: + +.. jupyter-execute:: + + custom = xls.composite( + scene, rescale_to=[0, 0.4], bands=("nir", "swir1", "red"), + ) + + fig, ax = plt.subplots(1, 1, figsize=(10, 6)) + custom.plot.imshow(ax=ax) + ax.set_title(custom.attrs["long_name"].title()) + ax.set_aspect("equal") + plt.show() -**With this, you can now make composites using any other band combination you -may want!** +This particular composite is great at distinguishing between built structures +in the city and along the canals (light green), the water ways (dark blue and +purple), and the forest (orange). diff --git a/xlandsat/datasets.py b/xlandsat/datasets.py index 7494c26..fd638ee 100644 --- a/xlandsat/datasets.py +++ b/xlandsat/datasets.py @@ -259,7 +259,8 @@ def fetch_manaus(untar=False): Download a sample scene from Manaus, Brazil Manaus is located in the Brazilian Amazon. The scene shows a part of the - city and the meeting of the Amazon and Negro rivers. + city and the meeting of the Solimões and Negro rivers to form the Amazon + river. This is a cropped version of a Landsat 9 scene from 2023/07/23, during the annual Amazon river floods.