index.qmd

---
pagetitle: "Raster data handling with Python"
author: "Jan Verbesselt, Jorge Mendes de Jesus, Aldo Bergsma, Dainius Masiliunas, David Swinkels, Judith Verstegen, Corné Vreugdenhil, Arno Timmer"
date: today
engine: knitr
format:
  html:
    theme: simplex
    highlight-style: zenburn
    toc: true
    toc-location: right
    css: styles.css
    lightbox: true
execute:
  eval: true
  echo: true
---

```{r}
#| include: false
# This is necessary to show outputs for students, 
# the include is false hides this for the students. 
reticulate::py_require("numpy")
reticulate::py_require("matplotlib")
reticulate::py_require("owslib")
reticulate::py_require("geopandas")
reticulate::py_require("rasterio")
reticulate::py_require("rasterstats")
reticulate::py_require("affine")
reticulate::py_require("osmnx")
```


[[WUR Geoscripting](https://geoscripting-wur.github.io/)]{.page-header-title} <img src="https://geoscripting-wur.github.io/figs/WUR_RGB_standard_2025.svg" alt="WUR logo" class="page-header-logo"/>

# Raster data handling with Python

## Introduction

Today we will work with Python packages for spatial raster analysis. Python has some dedicated packages to handle rasters:

* [OWSLib](https://geopython.github.io/OWSLib/) allows us to download geospatial raster data from Web Coverage Services
* [GDAL](https://gdal.org/api/python_bindings.html) is a powerful library for reading, writing and warping raster datasets
* [Rasterio](https://Rasterio.readthedocs.io/en/latest/) reads and writes geospatial raster data
* [rasterstats](https://pythonhosted.org/rasterstats/) summarizes geospatial raster datasets based on vector geometries
* [NumPy](http://www.numpy.org/) is a fundamental package for scientific computing, such as array (thus raster) calculations

## Learning objectives

- Be able to read spatial raster data from web services and files
- Be able to write spatial raster formats to disk
- Know how to apply basic operations on raster data, such as arithmetics
- Be able to plot spatial raster data with Matplotlib

## Setting up the Python Environment
Like in the previous tutorials, we will create a pixi environment. Do this by creating a directory for this tutorial, copying the following text in a `pixi.toml` file and run `pixi install`. 

```toml
[workspace]
authors = ["Your Name <your.email@example.com>"]
channels = ["conda-forge"]
name = "raster"
platforms = ["linux-64"]
version = "0.1.0"

[tasks]

[dependencies]
python = "*"
numpy = "*"
matplotlib = "*"
owslib = "*"
gdal = "*"
geopandas = "*"
rasterio = "*"
rasterstats = "*"
affine = "*"
ipykernel = "*"
```

Create the necessary directories for this tutorial (for this we use the nifty [Pathlib package](https://docs.python.org/3/library/pathlib.html), the standard path handling library in python!  ):

```{python}
from pathlib import Path

base = Path.home() # This references the location of your script
(base / "data").mkdir(exist_ok=True)
```


From now on, if you see the format combining strings and variables together with `/`'s like: `base / 'some_dir_name'`, it is Pathlib based path handling. 

Use this directory for the rest of this tutorial. 

# Reading raster data and accessing metadata
There are many ways to work with raster data, you have been introduced to some in R and bash. In Python we can do similar things, often using the same underlying software (gdal for example).

Before we will look into how to deal with local files and data processing, we will introduce web services and will show you how to access and store remote (online hosted) data locally. This is useful, because we can rely on the owner of the data (governments for example) updating the data, making our analysis pipelines make use of up to date data. Think about a yearly update of some dataset, where the source data is renewed 'automagically' each year.

Conveniently, the data that we will access and store locally, will be used throughout this tutorial to introduce libraries that can be used for reading, writing and processing raster data in Python.

Lastly, we will look a more recent development in the area of analyzing (sometimes very large and multidimensional) remotely stored raster datasets. 

## From a Web Coverage Service

A Web Coverage Service (WCS) loads raster data in a similar way as Web Feature Services (WFS) load vector data. [Web Coverage Services](https://www.ogc.org/standards/wcs) are a standard by the Open Geospatial Consortium and allow the downloading of geospatial raster data with multiple types of format encoding: GeoTIFF, netCDF, JPEG2000 etc. A [Web Map Service](https://www.ogc.org/standards/wms) [WMS] is another dataservice for raster data. An important difference between them is that a WCS serves the raw raster data (such as elevation in meters for Digital Elevation Models), and a WMS serves visualizations of data, for example a nice elevation basemap with 3 bands for Red Green and Blue (RGB) instead of the raw data.

Today we will work with elevation rasters. More specifically, we will have a look at the WCS of the AHN dataset. AHN stands for "Actueel Hoogtebestand Nederland" and is a Digital Elevation Model [DEM] that covers the Netherlands. Access the web coverage service to have a look at the contents. For this we will use OWSLib:

```{python}
from owslib.wcs import WebCoverageService

# Access the WCS by proving the url and optional arguments
wcs = WebCoverageService('https://service.pdok.nl/rws/ahn/wcs/v1_0?SERVICE=WCS&request=GetCapabilities', version='1.0.0')

# Print to check the contents of the WCS
print(list(wcs.contents))
```

Running the last line of code shows that the Web Coverage Service of the AHN3 contains two rasters with the identifiers `dtm_05m` and `dsm_05m`, corresponding to a  0.5m resolution DSM, and 0.5m resolution DTM. Raster data of AHN has the projected coordinate system `RD_New (EPSG: 28992)`. The difference between a DEM, DSM and DTM is explained in this question on [GIS StackExchange](https://gis.stackexchange.com/questions/5701/what-is-the-difference-between-dem-dsm-and-dtm/5704) (Thanks [@underdark](https://gis.stackexchange.com/users/187/underdark)!).

We can also check what types of operations are available for this WCS:

```{python}
# Get all operations and print the name of each of them
print([op.name for op in wcs.operations])
```

You will see that the Web Coverage Service allows accessing the data (`GetCoverage`), the metadata (`DescribeCoverage`), and the capabilities (`GetCapabilities`). These are all standard protocols defined by the OGC. 

Several functions are available to access specific metadata of each individual raster, for example:

```{python}
# Take the 0.5m DSM as an example
cvg = wcs.contents['dsm_05m']

# Print supported reference systems, the bounding box defined in WGS 84 coordinates, and supported file formats
print(f'Available CRS options: {cvg.supportedCRS}')
print(f'Bounding box is WGS84: {cvg.boundingBoxWGS84}')
print(f'Supported formats: {cvg.supportedFormats}')
```

Let us have a look at the data itself. Calling the GetCoverage function will send a request to the server where the data is stored and download the data. As we do not want to overload the web service, we call this function once, download the data and store it locally. 

Download the Digital Surface Model [DSM], which we can access through the 'dsm_05m' identifier, and Digital Terrain Model [DTM], which is the 'dtm_05m' version, to a local file.

```{python}
# Define a bounding box in the available crs (see before) by picking a point and drawing a 1x1 km box around it
x, y = 174100, 444100
bbox = (x - 500, y - 500, x + 500, y + 500)

# Request the DSM data from the WCS
response = wcs.getCoverage(identifier='dsm_05m', bbox=bbox, format='GEOTIFF',
                           crs='urn:ogc:def:crs:EPSG::28992', resx=0.5, resy=0.5)

# Write the data to a local file in the 'data' directory
with open(base / "data" / "AHN3_05m_DSM.tif", 'wb') as file:
    _ = file.write(response.read()) # file.write returns the number of bytes written

# Do the same for the DTM
response = wcs.getCoverage(identifier='dtm_05m', bbox=bbox, format='GEOTIFF',
                           crs='urn:ogc:def:crs:EPSG::28992', resx=0.5, resy=0.5)

with open(base / "data" / "AHN3_05m_DTM.tif", 'wb') as file:
    _ = file.write(response.read())
```

That's a very short introduction to webservices and accessing with Python. The datasets that you now stored locally we will use in the rest of the tutorial. Before continuing, please check if these steps were successful (*Check if files have been written into corresponding directory*).


## From a file with GDAL
[GDAL](https://gdal.org/) is the fundament under most spatial operations, reading and writing data and a lot more. GDAL is a C++ library that can handle raster and vector geospatial data formats. There are APIs to access GDAL functionality in Python, but also in  Java, R and C. When opening a raster file in GDAL, the object has a hierarchical structure starting at the Dataset level. A Dataset has a Geotransform (metadata) and can contain one or more Bands. Each Band has a Data array and potentially Overviews.

![GDAL class structure, adapted from Garrard, 2016, Geoprocessing with Python.](./images/gdal.png){fig-alt="gdal structure" width="90%"}

As you can see handling spatial data with GDAL can become a bit confusing. But since it is so fundamental it is good to get familiar with the very basics of GDAL. Luckily there are packages that build upon GDAL that make life a bit easier, but they all use GDAL under the hood. More on those packages later, first let's have a look at opening geotiff in GDAL.

Let us open the file we just saved. You will see you first get the dataset (even though there is only one), before the data array can be accessed.

```{python, eval=FALSE}
from osgeo import gdal

# Open dataset, gdal automatically selects the correct driver
ds = gdal.Open(str(base / "data" / "AHN3_05m_DSM.tif"))

# Get the band (band number 1)
band = ds.GetRasterBand(1)

# Get the data array
data = band.ReadAsArray()
print(data)

# Delete objects to close the file
ds = None
```

::: {.callout-note appearance="simple"}
**Question 1**: What happening with `ds = None` at the end of your script? Is it important? What may happen if you do not do that?
<details>
  <summary>*Click for answer*</summary>
  Setting `ds` to None, makes python forget the dataset, closing the dataset, it`s a C thing. Keeping files open may leave you vulnerable to losing data, (Geo)Pandas manages resources under the hood so you don't explicitly need to close files, but for the case of GDAL, and as you will later see, Rasterio, it's important to close your files. Another way to handle this is using contexts, for this we would open them with a context manager `with open ...`. More on that later.
</details>
:::

## From a file with Rasterio
GDAL can be tricky to work with but it allows for a lot of flexibility. There are tools that build upon GDAL, 
making integration with other software more intuitive. One of them is Rasterio. As [Rasterio explains on their website](https://Rasterio.readthedocs.io/en/latest/intro.html),
the goal of Rasterio is to the standard geospatial abstraction library that can use modern python features, and frees users from C related pitfalls.
It allows to read and write the same raster formats as GDAL can, provides raster processing functions, 
and can be integrated with other packages using NumPy and GeoJSON, and allows for easy visualization based upon MatPlotLib.

The rest of the tutorial below is a complete route of handling a raster dataset. We will use the DEMs from the WCS for our study area, handle it with Rasterio,
calculate new information (CHM), overlay it with vector data representing buildings and visualize it. 
Let us read in the raster data we just stored from the WCS with Rasterio and plot it with `rasterio.plot`:

```{python}
import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt

# Open the two rasters
dsm = rasterio.open(base / "data" / "AHN3_05m_DSM.tif", driver="GTiff")
dtm = rasterio.open(base / "data" / "AHN3_05m_DTM.tif", driver="GTiff")

# Metadata functions from Rasterio
print(dsm.meta)
print(dtm.meta)

# Plot with rasterio.plot, which provides Matplotlib functionality
plt.figure(figsize=(5, 5), dpi=300)  # adjust size and resolution
show(dsm, title='Digital Surface Model', cmap='gist_ncar')
```

::: {.callout-note appearance="simple"}
**Question 2**: Adjust the code above to take a look at the DTM. Note the gaps that appear. What are these gaps?
<details>
  <summary>*Click for answer*</summary>
  These are buildings.
</details>
:::

The metadata shows the driver, datatype, nodata value, width of raster in number of cells, height of raster in number of cells, number of raster bands in the dataset, coordinate reference system, and transformation values. 

Rasterio allows us to access the data values in a raster as NumPy arrays. We can use the `.read()` method to access these:

```{python}
# Rasterio object
print(type(dsm))

# Read, show object type and data
dsm_data = dsm.read(1)
print(type(dsm_data))
print(dsm_data)
```

[NumPy](https://numpy.org/doc/stable/user/whatisnumpy.html) is a very common Python Library providing support for large, multi-dimensional arrays and fast mathematical operations. Because it is so common, other libraries will be able to read NumPy arrays. This makes it possible to read raster data into other libraries, such as the machine learning in scikit-learn.


# Processing raster data
Now that we have some options to read data into python, let's have a look at what we can do with it. 

## Creating a Canopy Height Model

A Canopy Height Model (CHM) gives an indication of the height of trees and/or buildings. It can be created by subtracting a Digital Terrain Model from a Digital Surface Model. In the resulting raster, each cell value represents the height above the underlying surface topography.

```{python}
import numpy as np

# Access the data from the two rasters
dsm_data = dsm.read()
dtm_data = dtm.read()

# Set our nodata to np.nan (this is important for later)
dsm_data[dsm_data == dsm.nodata] = np.nan
dtm_data[dtm_data == dtm.nodata] = np.nan

# Close open rasterio datasets so their files can be removed
dsm.close()
dtm.close()
```

Similar to a gdal dataset, we have to close rasterio datasets. If you don't do this, you are firstly not able to delete the files but it can also lead to other strange behaviour.

Earlier, we noticed that the DTM included gaps. Let's first fill these gaps using the `fillnodata()` function from `rasterio.fill`. For more information, see the [documentation](https://rasterio.readthedocs.io/en/latest/api/rasterio.fill.html).

```{python}
from rasterio.fill import fillnodata

# Create a mask to specify which pixels to fill (0=fill, 1=do not fill)
dtm_mask = dtm_data.copy()
dtm_mask[~np.isnan(dtm_data)] = 1
dtm_mask[np.isnan(dtm_data)] = 0

# Fill missing values
dtm_data = fillnodata(dtm_data, mask=dtm_mask)
```

Now, we can create a CHM ([remember](https://gis.stackexchange.com/questions/5701/what-is-the-difference-between-dem-dsm-and-dtm/5704?)): 

```{python}
# Subtract the NumPy arrays 
chm = dsm_data - dtm_data

# Check the resulting array
print(chm)

# Copy metadata of one of the rasters (does not matter which one)
kwargs = dsm.meta 

# Save the chm as a raster
with rasterio.open(base / "data" / "AHN3_05m_CHM.tif", 'w', **kwargs) as file:
    file.write(chm.astype(rasterio.float32))

```

The actual raster arithmetic is happening in the first line: `chm = dsm_data - dtm_data`. dsm_data and dtm_data are plain numpy arrays, without spatial information. If we want to write this data back to a `.tif`, we need to do include this spatial information. All this information is stored in the metadata of a rasterio object. Since we subtracted two rasters with the same spatial dimensions (the same number of rows, columns and the same CRS for example), we can re-use this metadata to write a raster with the new data. Run the following code block to see what exactly is in this metadata: 

```{python}
# We use pprint instead of print to pretty print a dictionary 
import pprint
pprint.pp(dsm.meta)
```

::: {.callout-note appearance="simple"}
**Question 3**: Where is the CHM the highest in the study area? Is it what you expected?
<details>
  <summary>*Click for answer*</summary>
  Think about where you have the most forests on campus.
</details>
:::

We have now applied the basic raster arithmetic and created a Canopy Height Model!

## Computing heights of buildings

Using our CHM, let's determine the average heights of the buildings in our study area. The first step is knowing where the buildings are. For this we use the building data from the BAG Web Feature Service that we also used in the vector tutorial. Note that we make use of the `bbox` from an earlier codeblock for this.

```{python}
import geopandas as gpd
import json
from owslib.wfs import WebFeatureService

# Get the WFS of the BAG
wfsUrl = 'https://service.pdok.nl/lv/bag/wfs/v2_0'
wfs = WebFeatureService(url=wfsUrl, version='2.0.0')
layer = list(wfs.contents)[0]

# Get the features for the study area
# notice that we now get them as json, in contrast to before
response = wfs.getfeature(typename=layer, bbox=bbox, outputFormat='json')
data = json.loads(response.read())

# Create GeoDataFrame, without saving first
buildings_gdf = gpd.GeoDataFrame.from_features(data['features'])

# Set crs to RD New
buildings_gdf.crs = 28992
```

Getting the average per building polygon can be done using a common operation called zonal statistics. There are several implementations of this in python, we will use the [Rasterstats](https://pythonhosted.org/rasterstats/) library. [`Rasterstats.zonal_stats`](https://pythonhosted.org/rasterstats/) has 3 inputs, vector data, rasterdata and the type of statistic that needs to be calculated. See the vector and raster datasources pages in the documentation to see what inputs are supported for those. We will use a `GeoDataFrame` and a `.tif` file for this task. Using the `geojson_out` parameter we make it output a [GeoJSON](http://geojson.org/).

```{python}
import rasterstats as rs

# Apply the zonal statistics function with gdf and tif as input
chm_buildings = rs.zonal_stats(buildings_gdf, base / "data" / "AHN3_05m_CHM.tif", prefix='CHM_', geojson_out=True)

# Convert GeoJSON to GeoDataFrame
buildings_gdf = gpd.GeoDataFrame.from_features(chm_buildings)

# Check the added attributes with a prefix 'CHM_'
print(buildings_gdf['CHM_mean'])
```

A quick visualization shows us the heights derived from the raster data on the map:

```{python}
# Create one plot with figure size 10 by 10
fig, ax = plt.subplots(1, figsize=(10, 10))

# Customize figure with title, legend, and facecolour
ax.set_title('Heights above ground (m) of buildings on the WUR campus')
buildings_gdf.plot(ax=ax, column='CHM_mean', k=6,
                   cmap=plt.cm.viridis, linewidth=1, edgecolor='black', legend=True)
ax.set_facecolor("lightgray")

# Visualize figure
plt.show()
```


## Other functionality
Note that this tutorial only scratches the surface of the possibilities of Rasterio. It can do most if not all things you did in `R` in the Vector - Raster tutorial. Rasterio for example also allows you to do [masking](https://rasterio.readthedocs.io/en/latest/topics/masking-by-shapefile.html), [reprojecting](https://rasterio.readthedocs.io/en/latest/topics/reproject.html), and [resampling](https://rasterio.readthedocs.io/en/latest/topics/resampling.html).

# Writing raster data to a file
As you've seen before, to store the NumPy array as a raster file, Rasterio needs the accompanying metadata. It is possible to use the metadata of an existing raster (which we did before), but it is also possible to create it from scratch. 

To create metadata from scratch, the CRS can be defined with a function from Rasterio and the transformation can be defined using Affine. Affine is a Python module that facilitates [affine transformations](https://www.quora.com/In-an-intuitive-explanation-what-is-an-affine-transformation-of-image), i.e. scaling, rotating, mirroring or skewing of images/rasters/arrays.

Rasterio can write most [raster formats from GDAL](https://www.gdal.org/formats_list.html). [The developers recommend using GeoTiff driver](https://github.com/mapbox/rasterio/issues/731) for writing as it is the best-tested and best-supported format. 

```{python}
import affine

# Specify the components of the crs (we know them from the DSM)
kwargs = {'driver': 'GTiff',
          'dtype': 'float32',
          'nodata': np.nan,
          'width': 2000,
          'height': 2000,
          'count': 1,
          'crs': rasterio.crs.CRS.from_epsg(28992),
          'transform': affine.Affine(0.5, 0.0, 173600.0, 0.0, -0.5, 444600.0)}

# Write the raster file
with rasterio.open(base / "data" / "AHN3_05m_CHM_affine.tif", 'w', **kwargs) as file:
    file.write(chm.astype(rasterio.float32))
```

# More on raster data visualization
Raster data can be visualized by passing NumPy arrays to Matplotlib directly or by making use of a method in Rasterio that accesses Matplotlib for you. Using Matplotlib 
directly allows more flexibility, such as tweaking the legend, axis and labels, and is more suitable for professional purposes. The visualization using Rasterio requires 
less code and can give a quick idea of your raster data. We show both approaches below. Let's first make a visualization of the DSM using Matplotlib:

```{python}
# Create one plot with figure size 10 by 10
fig, ax = plt.subplots(figsize=(10, 10), dpi=200)

# rasterio's imshow() is the main raster plotting method in Matplotlib
# Ensure an equal scale in the x and y direction
dsmplot = ax.imshow(dsm_data[0], cmap='Oranges', extent=bbox, aspect='equal')

# Title (do not do this for a scientific report, use a caption instead)
ax.set_title("Digital Surface Model - WUR Campus", fontsize=14)

# Add a legend (colourbar) with label
cbar = fig.colorbar(dsmplot, fraction=0.035, pad=0.01)
cbar.ax.get_yaxis().labelpad = 15
cbar.ax.set_ylabel('Height (m)', rotation=90)

# Hide the axes
ax.set_axis_off()
plt.show()
```

If you do not like the orange colourmap of Matplotlib, it is also possible to pick [another colourmap](https://Matplotlib.org/examples/color/colormaps_reference.html).

The second approach with Rasterio only requires one line of code to make a plot. By creating subplots, the figures can be combined (this can be done with Matplotlib directly as well).

```{python}
# Figure with 1 row with 3 columns (1,3), so 3 subplots. 
# They are unpacked directly, plt.subplot returns fig and a tuple of the 3 axes
fig, (axdsm, axdtm, axchm) = plt.subplots(1, 3, figsize=(15, 7), dpi=200)

# Populate the three subplots with raster data
show(dsm_data, ax=axdsm, title='DSM')
show(dtm_data, ax=axdtm, title='Filled DTM')
show(chm, ax=axchm, title='CHM')
plt.show()
```

Rasterio can also create simple histograms by calling functions of Matplotlib. Let's see if the 3 different elevation models are different comparing their histograms. For a bit of efficiency, we will stack the three models and store them as a three band raster. We will open it again and look at its histogram. 

```{python}
from rasterio.plot import show_hist
import numpy as np

# Stack the three arrays along a new axis to create a 3-band raster
# Ensure all arrays have the same shape and are 2D
three_band_array = np.stack([dsm_data, dtm_data, chm], axis=0).squeeze()

# Copy metadata and update for 3 bands
three_band_meta = dsm.meta.copy()
three_band_meta.update(count=3)

# Write the 3-band raster to file
with rasterio.open(base / "data" / "AHN3_05m_3band.tif", 'w', **three_band_meta) as dst:
  dst.write(three_band_array.astype(rasterio.float32))

# Open it again and plot the histogram
with rasterio.open(base / "data" / "AHN3_05m_3band.tif") as raster:
  show_hist(raster, bins=50, lw=0.0, stacked=False, alpha=0.3, 
            histtype='stepfilled', title="Histogram", label = ['dsm', 'dtm', 'chm'])

```


::: {.callout-note appearance="simple"}
**Question 4**: What is represented on the x and y axis? The default axis labels are DN (x) and Frequency (y); if you were to change them, what labels would you pick to better reflect the content of the plots?
<details>
  <summary>*Click for answer*</summary>
  The y axis represents the count of pixels. Meanwhile the x axis represents the pixel's DN (digital value), in this tutorial since we are looking at elevation this value is actually meters. For example, in the first plot (DSM) you can see that most pixel values are in the 10 to 15 meter range
</details>
:::


## STAC 
In the examples above, we have looked at single images or relatively small areas. In the past years the amount of data that we are capturing has been growing. Web services are not suitable for sharing very large quantities of multimodal or multidimensional data. Searching and accessing these volumes of data can be done using the Spatial Temporal Asset Catalogue (STAC). 

This is a short (and incomplete) description of what STAC is and how to use it. For a more elaborate tutorial visit the dedicated [STAC tutorial](https://geoscripting-wur.github.io/STAC/). 

STAC is a standard with a strong community working on enabling easier access to data about our planet. It provides a standardized structure for accessing spatiotemporal datasets. The starting point is usually a STAC catalog, for example [https://stac.ecodatacube.eu/](https://stac.ecodatacube.eu/?.language=en). This website is a visual representation of a `catalog.json` file, that can be accessed by clicking on the source button in the top right. It is possible to search the catalog programmatically, and is discussed in the [STAC Tutorial](https://geoscripting-wur.github.io/STAC/). An alternative is to browse through the viewer and find relevant datasets there.

A *STAC catalog* consists of a links to *collections*, *items* or other catalogs. A STAC collection is a collection of STAC items with similar traits. For example [`https://s3.ecodatacube.eu/arco/stac/oc_iso.10694.1995.mg.cm3/collection.json`](https://s3.ecodatacube.eu/arco/stac/oc_iso.10694.1995.mg.cm3/collection.json) (that can be found through the ecodatacube catalog). This collection, is a collection of STAC items regarding the soil organic carbon density (SOCD). In this case it links to  datasets at different depths, at different timeframes. 

Generally these separate datasets can be found by a url stored in the items. These urls point, in the case of raster data, to a Cloud Optimized GeoTiff (COG). This is a Tiff file as we know it, it also ends at `.tiff`, but it contains some indexing and other metadata that  makes the file searchable remotely instead of having to download the entire file before being able to for example mask a part of it. 

```{python}
import rasterio
from rasterio.mask import mask
from shapely import box
import matplotlib.pyplot as plt

# Redefine the bounding box in the available crs (see before), now in EPSG:3035 to match the raster
x, y = 4024900, 3215843.5
bbox = box(x - 3579, y - 2134.5, x + 3579, y + 2134.5) # In EPSG:3035 - same as raster

# get url from item asset link
asset_url = 'https://s3.ecodatacube.eu/arco/ndvi_glad.landsat.ard2.seasconv.m.yearly_p25_30m_s_20220101_20221231_eu_epsg.3035_v20231127.tif'

# open raster directly from S3 (requires rasterio with HTTP enabled, which is default)
with rasterio.open(asset_url) as src:
    # crop the raster with the polygon
    out_image, out_transform = mask(src, [bbox], crop=True)

# plot the cropped raster
plt.imshow(out_image[0], cmap="viridis")
plt.title("Cropped raster over Wageningen")
plt.colorbar(label="Value")
plt.show()
```

# More info
  * [Tutorial working with rasters in Python with Rasterio](https://geohackweek.github.io/raster/04-workingwithrasters/)
  * [Tutorial working with raster in Python with GDAL (for Python 2)](https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html)
  * [Landsat satellite images](https://earthexplorer.usgs.gov/)
  * [Resampled landsat satellite images](http://espa.cr.usgs.gov/index/)
  * [Sentinel satellite images](https://scihub.copernicus.eu/dhus/#/home)
  


```{python}
#| include: false

# Clean up directories after rendering
import shutil

# Close open rasterio datasets so their files can be removed
dsm.close()
dtm.close()

shutil.rmtree(base / "data")

```