Sampling the Cube

In the previous chapter, you learned how to ingest and visually verify a dataset. Now, we move from passive inspection to active selection. While conceptually a cube is a coherent structure of observations, its real power is unlocked when you can query it deliberately to solve scientific problems.

The shift in this chapter is fundamental: you will move away from “index guessing”—manually counting matrix positions—and toward a readable scientific workflow based on named dimensions, coordinate labels, and conditions. This approach makes your code easier to read, easier to debug, and much closer to the scientific questions you are actually asking.

We will progress from basic coordinate selection to advanced conditional masking and vectorized sampling using the air_temperature tutorial dataset.

1. Why querying the cube matters¶

A data cube is powerful because it stores thousands (or millions) of observations in one aligned structure. However, that structure is only useful if we can ask targeted questions of it. In spatial data science, analysis is almost always selective by nature.

A scientist, planner, or analyst rarely asks to see the “full array.” Instead, they “interrogate” the cube with questions involving specific constraints:

• Temporal Patterns: What was the air temperature on a specific morning in June 2013?

• Local Trends: How does the temperature at a specific weather station change through time?

• Regional Summaries: What is the mean value inside a specific administrative boundary or study area?

• Condition-based Filtering: Where in the region did the temperature exceed 30°C today?

Xarray provides the tools to move from a large, abstract cube to a small, interpretable subset that directly matches your research question. If the previous chapter taught you how to read the cube, this chapter teaches you how to interrogate it.

2. Position and label based selection¶

Xarray provides four distinct ways to index and select data. The most important distinction is between navigating by physical matrix indices and navigating by coordinate values. By using labeled dimensions, Xarray eliminates the mental overhead of remembering whether “axis 0” refers to time, latitude, or a specific sensor band.

.isel for position based selection¶

The method .isel (integer selection) is the direct descendant of NumPy style indexing. It is 0-based and allows you to select specific elements by their numerical position along a dimension.

Positional indexing is most useful when your query is relative to the array structure rather than the physical world, such as selecting the last image in a stack or a specific subset of pixels for a training window.

.sel for label based selection¶

Scientific analysis is rarely defined by matrix indices; instead, it is defined by real world parameters. The .sel method uses dimension labels to automatically look up the correct indices for you.

Because coordinates are backed by an index, these lookups are highly efficient. You can also use slice() with .sel to filter data between specific geographic or temporal bounds.

Nearest neighbor lookups¶

In real-world Geospatial Data Science, your target coordinates (e.g., a weather station location) might not perfectly align with the center of a grid cell. While standard indexing would fail, .sel supports inexact lookups via the method and tolerance parameters.

• method="nearest": Selects the value of the closest pixel without interpolation.

• tolerance: Limits the maximum distance for a valid match, preventing the selection of data that is too far from your target location.

Differences from NumPy behavior¶

While Xarray feels like NumPy, it behaves differently when you index with multiple arrays simultaneously.

In NumPy, indexing with two arrays (e.g., arr[[0,1], [0,1]]) performs pointwise indexing, returning elements at (0,0) and (1,1). In Xarray, the default behavior is orthogonal (or outer) indexing. This means Xarray selects the entire Cartesian product of the indices, resulting in a 2D sub array rather than a 1D list of points.

When to use which¶

To maintain clean and reproducible code, follow this practical rule:

• Use .isel when you care about array positions (e.g., the “first” or “last” element).

• Use .sel when you care about real coordinates or dates (e.g., “Seattle” or “2013-05-12”).

These methods work identically for both DataArray and Dataset objects, allowing you to index all variables in a dataset simultaneously.

3. Slicing time¶

Time is arguably the most important dimension in a data cube. Xarray is exceptionally elegant here because time coordinates are treated as first class datetime objects rather than simple strings or numerical indices. This temporal awareness allows you to navigate years, months, and specific moments with high precision.

Selecting specific moments¶

You can query the cube for an exact moment in time using a full timestamp string. In a 3D cube, selecting one moment across the time dimension returns a 2D spatial slice (a map).

Natural shorthand for months and years¶

One of Xarray’s most convenient features is its ability to interpret partial datetime strings. This allows you to grab entire blocks of data without needing to define start and end dates manually.

Temporal intervals with slice¶

To select a continuous range of time, use the Python slice function. This is the standard way to isolate specific seasons or study periods.

Components and fancy temporal indexing¶

Behind the scenes, Xarray provides a .dt accessor that allows you to reach into the datetime objects and extract specific components like month, day, year, or even season.

This enables “fancy indexing” where you select data based on a condition rather than a continuous range. For example, if you want to perform a climatological study on a specific month across multiple years, you can filter the cube by the month component.

This logic also works for grouping by the season component (DJF, MAM, JJA, SON), providing a bridge to the advanced temporal aggregation and seasonal analysis covered later in this course.

4. Slicing space¶

Spatial selection mirrors the logic of temporal slicing but operates on the horizontal axes of the cube. By using coordinate labels instead of axis numbers, you avoid the mental overhead of remembering if “axis 0” refers to latitude or longitude.

Horizontal cross sections¶

The simplest spatial slice is a single row or column. This reduces a 3D cube into a 2D “vertical” slice of time and space.

Bounding boxes¶

A rectangular subset, or bounding box, is the most common way to isolate a specific study area. You define these windows using the slice() function within .sel().

Point extraction¶

Selecting a single cell by its coordinates is the spatial equivalent of selecting one exact timestamp. If you do not specify a time constraint, this operation isolates a location across the entire temporal dimension, effectively turning a spatial cube into a 1D time series.

Reusable study areas¶

In professional Geospatial Data Science workflows, it is best practice to define your study area once and persist it as a new object. This makes your code more readable and prevents you from accidentally running heavy calculations on the full global dataset when you only need a local subset.

By persisting subsets, you move from general exploration to a focused analysis of a specific geographic region.

4. Conditional selection and masking¶

While .sel() and .isel() navigate the cube using its “outer” labels (where the data is), you often need to interrogate the “inner” values (what the data says). This is the role of boolean masking, a technique that enables filtering values based on specific scientific criteria rather than just their location.

The concept of a boolean mask¶

A Boolean mask is a binary array of the same shape as your data, containing only True or False values. It acts as a spatial filter: when applied to a DataArray, only the elements corresponding to True are retained.

Filtering with .where()¶

Xarray provides the .where() method to apply these masks. Unlike coordinate selection, .where() typically preserves the shape of the original data. Values that do not meet the condition are replaced with a mask value—by default, NaN (Not a Number).

If you need to replace masked values with a specific placeholder instead of NaN, you can provide a second argument:

Clipping with drop=True¶

By default, .where() keeps the entire grid, even if large portions are empty. To shrink the cube to the smallest possible bounding box that contains your valid data, use the drop=True parameter. This is particularly useful for “clipping” a dataset to a specific feature or anomaly.

Complex logical conditions¶

You can specify multiple masking conditions within a single .where() statement using bitwise operators: & (and) and | (or). This allows for the selection of specific value ranges or complex geographic windows.

Discrete selection with .isin()¶

When your data contains discrete categories rather than continuous measurements—such as land cover types or quality status flags—the .isin() method is the preferred tool. It allows you to check for a list of specific values simultaneously.

Strategic considerations for workflows¶

Masking is central to professional Geospatial Data Science. Typical applications include:

• Cloud Removal: Masking pixels flagged as cloudy in satellite metadata.

• Hydrological Modeling: Masking out water bodies or oceans from a terrain analysis.

• Anomaly Detection: Isolating only the most extreme temperature or precipitation values.

Comparison: Selection vs Masking¶

Concept Check: The Subsetting Showdown¶

You need to remove all temperatures below 0°C from your dataset, leaving those pixels empty for future statistical calculations. Which method is the right tool for the job?

1. .isel()

2. .sel()

3. .where()

5. Extracting a point time series¶

Extracting a time series is the bridge between spatial maps and temporal trends. While maps show the distribution of a variable at one moment, a time series reveals how a specific location evolves throughout the entire duration of the dataset.

From 3D cube to 1D signal¶

In a data cube with dimensions (time, lat, lon), selecting one specific coordinate reduces the spatial dimensions to single values. This transforms the 3D structure into a one-dimensional signal indexed by time.

Handling coordinate grid mismatches¶

In real-world Geospatial Data Science, the coordinates of your observation sites (e.g., weather stations or cities) rarely align perfectly with the centers of a raster grid. If you attempt to select a coordinate that does not exist exactly in the index, Xarray will raise an error unless you specify a lookup strategy.

The method="nearest" parameter tells Xarray to find the single closest available grid cell.

Continuous estimation with interpolation¶

While nearest neighbor selection is fast and preserves original data values, it can lead to “blocky” results if you are sampling many points close together. To look at inferred values between existing grid cells, you should use the interpolation method.

Instead of snapping to the nearest pixel, interpolation uses neighboring observations (typically via a linear or cubic method) to estimate the value at your exact requested location.

Visualizing temporal trends¶

One of the most efficient features of Xarray is its ability to automatically detect the dimensions of a result. Because the extracted object is 1D, calling .plot() immediately generates a line chart with correctly formatted time labels.

By selecting one spatial location, you have successfully “interrogated” the 3D cube to extract a continuous signal ready for trend analysis, seasonality checks, or climatological modeling.

6. Sampling multiple locations and regional summaries¶

While single point extraction is a vital diagnostic tool, scientific analysis often requires comparing several distinct sites or generating a representative summary for an entire geographic area. This section moves from simple subsetting to pointwise indexing and spatial reductions.

Vectorized pointwise indexing¶

In professional workflows, you often need to sample the data cube at a specific collection of coordinates (for example, at the locations of several weather stations or urban centers).

Xarray handles this through Vectorized Indexing. If you provide selection coordinates as Xarray DataArrays that share a dimension name (such as "points"), Xarray will extract exactly those pairs rather than a 2D rectangular grid. This is significantly more efficient than looping through locations one by one.

Comparing site time series¶

Because the resulting object preserves the new "points" dimension, you can immediately visualize the differences between your selected locations using a line plot. This allows for a rapid assessment of whether different sites follow similar temporal trends.

Regional summaries and reductions¶

Sometimes individual pixels are too noisy or specific for your analysis. In these cases, you want to reduce a 3D spatial subset into a single representative 1D time series for a region.

This process is called a Spatial Reduction. By calculating the mean across both horizontal dimensions (lat and lon), you effectively collapse the spatial variability into a single regional signal.

Point, box, and polygon sampling¶

In this chapter, we have focused on point sampling (zero dimensional) and bounding box sampling (two dimensional rectangular). These tools cover the vast majority of exploratory sampling tasks.

For more complex analysis involving irregular geographic shapes (like political boundaries or watersheds), you will eventually use Zonal Statistics. This requires an additional step of rasterizing your polygons to align them with the cube, a topic we will explore in the next section on integrating vector and raster data.

7. Exercise: Scale Comparison¶

In this exercise, you will evaluate how well a single observation site represents its broader surrounding region. This task requires combining point selection, spatial slicing, and dimension reduction.

Task¶

1. Ingestion: Open the air_temperature tutorial dataset and extract the air variable.

2. Point Sampling: Select the coordinate at lat=40.0 and lon=260.0 to extract its full time series.

3. Regional Slicing: Define a spatial bounding box centered around that same location using lon=slice(255, 265) and lat=slice(45, 35).

4. Reduction: Calculate the mean air temperature of that box for every time step.

5. Visualization: Plot both the point time series and the regional mean on a single figure with appropriate labels and a legend.

6. Interpretation: Write two short sentences comparing the behavior of the point and the regional mean. Consider which series is more stable and why.

Starter code¶

Sample solution

import xarray as xr
import matplotlib.pyplot as plt

ds = xr.tutorial.open_dataset("air_temperature")
air = ds["air"]

# Point time series extraction using label-based selection
point_ts = air.sel(lat=40.0, lon=260.0)

# Bounding box selection and spatial mean reduction
# We reduce 'lat' and 'lon' simultaneously to get a 1D time series
region = air.sel(lon=slice(255, 265), lat=slice(45, 35))
region_mean = region.mean(dim=["lat", "lon"])

# Plotting the results together
fig, ax = plt.subplots(figsize=(10, 4))
point_ts.plot(ax=ax, label="Single Point (40N, 260E)", alpha=0.6)
region_mean.plot(ax=ax, label="Regional Mean (10° Box)", linewidth=2)

ax.set_title("Point Time Series vs. Regional Mean")
ax.set_ylabel("Temperature (K)")
ax.legend()
plt.show()

Interpretation The point time series is more variable because it reflects the specific atmospheric fluctuations of a single grid cell. The regional mean is smoother because it averages over many neighboring cells, which cancels out local extremes and reduces overall noise in the signal.

8. Summary¶

Sampling is where the data cube transforms from a collection of observations into a targeted analytical tool. In this chapter, you moved from passive inspection to active interrogation of the cube.

You learned that Xarray supports three primary paths for navigating data:

1. Position based selection (.isel): Navigating the cube using 0-based integer indices, serving as the direct successor to NumPy indexing.

2. Label based selection (.sel): Querying the cube using meaningful, human-readable labels like specific dates or physical coordinates.

3. Condition based masking (.where): Filtering the data based on internal values, such as cloud flags or temperature thresholds, rather than its location.

Beyond basic selection, you also learned how to:

• Slice time and space: Use partial date strings or coordinate ranges with slice() to isolate specific intervals and geographic bounding boxes.

• Perform nearest neighbor lookups: Handle mismatches between your study coordinates and the raster grid using the nearest method and distance tolerances.

• Extract point time series: Reduce a 3D cube into a 1D temporal signal for site-specific analysis.

• Leverage vectorized indexing: Sample multiple distinct locations simultaneously by providing coordinate arrays with shared dimension names.

• Generate regional summaries: Reduce complex spatial subsets into interpretable indicators by calculating means across horizontal dimensions.

The key lesson of this chapter is analytical clarity. By using named dimensions and coordinates, Xarray allows you to move through a cube in a way that matches how scientists actually think about their data: by dates, places, regions, and conditions, rather than by anonymous integer axes.

Now that you can isolate exactly the right data, the next chapter covers how to combine, align, and perform map algebra between multiple different cubes to build complex Geospatial Data Science workflows.