<jupyter_start><jupyter_text>Accessing USGS Gap Land Cover data with the Planetary Computer STAC API  Environment setupThis notebook works with or without a Planetary Computer subscription key, but you will be given more permissive access to the data with a key. The Planetary Computer Hub is pre-configured to use your key.<jupyter_code>import planetary_computer

# Set the environment variable PC_SDK_SUBSCRIPTION_KEY, or set it here.
# The Hub sets PC_SDK_SUBSCRIPTION_KEY automatically.
# planetary_computer.settings.set_subscription_key("add your key here")<jupyter_output><empty_output><jupyter_text>Open and explore the USGS GAP collection[GAP/LANDFIRE](https://planetarycomputer.microsoft.com/dataset/gap) is a national land cover dataset from the [USGS](https://www.usgs.gov/core-science-systems/science-analytics-and-synthesis/gap).The Planetary Computer hosts the land cover dataset tiled into Cloud Optimized GeoTIFFs.<jupyter_code>from pystac_client import Client

catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
gap = catalog.get_child(id="gap")<jupyter_output><empty_output><jupyter_text>The GAP dataset has four distinct regions:- The continental United States (CONUS)- Alaska- Hawaii- Puerto RicoYou can see those four bounding boxes in the collection spatial extent:<jupyter_code>gap.extent.spatial.to_dict()<jupyter_output><empty_output><jupyter_text>Plotting land cover dataThe Mount Whitney area of California is an interesting area for investiating land cover classifications.The region includes developed areas, agriculture, forests, high alipine, and much more.Let's take a look at the region in the land cover dataset.<jupyter_code>bbox = [-118.35, 36.5, -118, 36.75]
search = catalog.search(collections="gap", bbox=bbox)
items = list(planetary_computer.sign(item) for item in search.get_items())
print(f"{len(items)} items found")<jupyter_output>1 items found<jupyter_text>In order to visualize the land cover data, we'll need to do a little work.The land cover data contains its own color table, so we need to read this color table and create a colormap that can be used with matplotlib.Also, because the source data is in a projected coordinate system, we'll need to transform our bounding box to that coordinate system for the merge.<jupyter_code>from matplotlib.colors import ListedColormap
from pystac.extensions.projection import ProjectionExtension
import rasterio
from rasterio.crs import CRS
from rasterio.warp import transform_bounds

assert items
item = items[0]

with rasterio.open(item.assets["data"].href) as dataset:
    colormap = dataset.colormap(1)

keys = list(colormap.keys())
keys.sort()
assert keys == list(
    range(0, len(keys))
)  # We make sure that every possible input value has an associated color.
color_list = list([float(v) / float(255) for v in colormap[key]] for key in keys)
cmap = ListedColormap(color_list)

projection = ProjectionExtension.ext(item)
dst_crs = CRS.from_string(projection.wkt2)
src_crs = CRS.from_string("EPSG:4326")  # WGS84 geodetic
projected_bbox = transform_bounds(src_crs, dst_crs, *bbox)<jupyter_output><empty_output><jupyter_text>Now we can plot the data, extracting the labels for the colors from the STAC item.<jupyter_code>from matplotlib import pyplot
from matplotlib.patches import Patch
import numpy
from rasterio.plot import show
from rasterio.merge import merge
from pystac.extensions.label import LabelExtension

assets = [item.assets["data"].href for item in items]
data, transform = merge(datasets=assets, bounds=projected_bbox)
data = numpy.squeeze(data)
color_data = numpy.moveaxis(
    cmap(data), [0, 1, 2], [1, 2, 0]
)  # rasterio.show.show wants bands in (band, x, y) order

label = LabelExtension.ext(item)
patches = [
    Patch(color=cmap(v), label=label.label_classes[0].classes[v])
    for v in numpy.unique(data)
]

fig, ax = pyplot.subplots(figsize=(12, 12))
ax.set_axis_off()


fig.legend(handles=patches, bbox_to_anchor=(1.4, 0.876))
_ = show(color_data, ax=ax, transform=transform, interpolation="none")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09079754601226994, 'fraction_numerical': 0.011533742331288344, 'mean_word_length': 6.125874125874126, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '478778', 'license': 'MIT', 'n_tokens_mistral': 1332, 'n_tokens_neox': 1261, 'n_words': 406, 'path': 'datasets/gap/gap-example.ipynb', 'repo_name': 'ianthomas23/PlanetaryComputerExamples', 'size': '4075'}