> For the complete documentation index, see [llms.txt](https://satdocs.hydrosat.com/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://satdocs.hydrosat.com/use-cases/industrial-tata-steel.md).

# Environmental compliance at Tata Steel (IJmuiden, Netherlands)

![Tata Steel IJmuiden](https://3702723385-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FBgOfOR7RYuK5Fm69JAWF%2Fuploads%2Fgit-blob-c7fcd947d53ba63ebaea3d4f75f1282702bce859%2Ftata_steel_hero.png?alt=media)

**In this case study, we’ll show how you can identify changes in industrial activity with a few simple analyses applied to Hydrosat data.**

### Introduction

Tata Steel IJmuiden is one of Europe's largest integrated steelworks and a major employer and exporter for the Netherlands, producing several million tonnes of steel a year at its site on the North Sea coast. In April 2026, the plant's DSP unit was shut down after air quality measurements showed chromium-6 emissions from one of its chimneys exceeding permissible limits, a serious environmental and public health concern. [Source: GMK Center](https://gmk.center/en/news/tata-steel-nederland-has-temporarily-shut-down-its-direct-sheet-plant-in-eemeden/).

Thermal satellite imagery offers a way to independently verify whether the DSP has genuinely gone offline, without relying solely on self-reported compliance data.

### Loading the L1B LIRI Collection

We'll start with a single Level-1B LIRI scene, `VZ02_L1b_20260302_132735`, taken before the DSP shutdown. Before diving into any analysis, we'll take a quick look at the images and their metadata.

```python
from pathlib import Path
import matplotlib.pyplot as plt
import rioxarray
import xarray as xr

# Define the scene ID and data directory for the imagery
DATA_DIR = Path("data")
SCENE_ID = "VZ02_L1b_20260302_132735"
```

```python
# Read LWIR bands
lwir1 = rioxarray.open_rasterio(DATA_DIR / SCENE_ID / "LWIR1.tiff", masked=True).squeeze("band", drop=True)
lwir2 = rioxarray.open_rasterio(DATA_DIR / SCENE_ID / "LWIR2.tiff", masked=True).squeeze("band", drop=True)

# Build an RGB preview from the individual RED/GREEN/BLUE bands
red = rioxarray.open_rasterio(DATA_DIR / SCENE_ID / "RED.tiff", masked=True).squeeze("band", drop=True)
green = rioxarray.open_rasterio(DATA_DIR / SCENE_ID / "GREEN.tiff", masked=True).squeeze("band", drop=True)
blue = rioxarray.open_rasterio(DATA_DIR / SCENE_ID / "BLUE.tiff", masked=True).squeeze("band", drop=True)
preview = xr.concat([red, green, blue], dim="band").assign_coords(band=["red", "green", "blue"])
```

#### Visualizing the LWIR Bands and Preview

`rioxarray` already masked out nodata pixels for us on load. L1B LIRI pixel values are scaled brightness temperature, so multiplying by the `0.01` scaling factor converts them to Kelvin, which we then convert to Fahrenheit for readability. We'll also show the full-resolution RGB preview alongside the two thermal bands for context.

```python
# Scaling factor for brightness temperature (from DN to Kelvin)
BT_SCALE = 0.01


# Convert brightness temperature from Kelvin to Fahrenheit
def kelvin_to_fahrenheit(da):
    return (da - 273.15) * 9 / 5 + 32


# Stretch each band independently to enhance contrast in RGB composite
def percentile_stretch(da, low=1, high=99):
    vmin, vmax = da.quantile(low / 100), da.quantile(high / 100)
    return ((da - vmin) / (vmax - vmin)).clip(0, 1)


# Convert LWIR bands to brightness temperature and stretch RGB preview
lwir1_bt = kelvin_to_fahrenheit(lwir1 * BT_SCALE)
lwir2_bt = kelvin_to_fahrenheit(lwir2 * BT_SCALE)
preview_stretched = preview.groupby("band").map(percentile_stretch)

# Plot LWIR brightness temperature and RGB preview
fig, axes = plt.subplots(1, 3, figsize=(19, 6))

for ax, da, label in zip(
    axes[:2], [lwir1_bt, lwir2_bt], ["LWIR1 (10.9 µm)", "LWIR2 (12.0 µm)"]
):
    da.plot.imshow(ax=ax, cmap="magma", add_colorbar=True, cbar_kwargs={"label": "°F"})
    ax.set_title(f"{label} Brightness Temperature")
    ax.set_aspect("equal")
    ax.axis("off")

preview_stretched.plot.imshow(ax=axes[2], rgb="band")
axes[2].set_title("RGB Preview")
axes[2].set_aspect("equal")
axes[2].axis("off")

fig.suptitle(SCENE_ID)
plt.tight_layout()
plt.show()
```

![](https://3702723385-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FBgOfOR7RYuK5Fm69JAWF%2Fuploads%2Fgit-blob-92484f81ffbdad669bd6893e6c38f774aa9c9f9c%2Foutput_8_0.png?alt=media)

Here we can see the broader context of the area, including the North Sea coastline and the harbor canal running through IJmuiden. Let's focus in on the steelworks for our next steps.

### Cropping to a Region of Interest

The full scene extends well beyond the steelworks. We'll crop LWIR1 to a region of interest (ROI) around the plant, loaded from a local file, so the same box can be reused across scenes with different projections or pixel grids.

```python
import json


# Function to extract ROI bounds from a GeoJSON file
def roi_bounds_from_geojson(path):
    coords = json.loads(Path(path).read_text())["features"][0]["geometry"]["coordinates"][0]
    lons, lats = zip(*coords)
    return (min(lons), min(lats), max(lons), max(lats))


# Extract ROI bounds from the GeoJSON file and clip the LWIR band to the ROI
ROI_BOUNDS = roi_bounds_from_geojson("assets/tata_steel.geojson")
lwir1_roi = lwir1_bt.rio.clip_box(*ROI_BOUNDS, crs="EPSG:4326")
```

```python
from matplotlib.patches import Rectangle
from rasterio.warp import transform_bounds

# Reproject the ROI bounds into the raster's CRS so the box lines up on the full scene
roi_minx, roi_miny, roi_maxx, roi_maxy = transform_bounds(
    "EPSG:4326", lwir1_bt.rio.crs, *ROI_BOUNDS
)

# Create a figure with two subplots: one for the full scene and one for the ROI
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

lwir1_bt.plot.imshow(ax=axes[0], cmap="magma", add_colorbar=True, cbar_kwargs={"label": "°F"})
axes[0].add_patch(
    Rectangle(
        (roi_minx, roi_miny),
        roi_maxx - roi_minx,
        roi_maxy - roi_miny,
        edgecolor="cyan",
        facecolor="none",
        linewidth=2,
    )
)
axes[0].set_title("LWIR1 (10.9 µm) - Full Scene")
axes[0].set_aspect("equal")
axes[0].axis("off")

lwir1_roi.plot.imshow(ax=axes[1], cmap="magma", add_colorbar=True, cbar_kwargs={"label": "°F"})
axes[1].set_title("LWIR1 (10.9 µm) - ROI")
axes[1].set_aspect("equal")
axes[1].axis("off")

fig.suptitle(f"{SCENE_ID}")
plt.tight_layout()
plt.show()
```

![](https://3702723385-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FBgOfOR7RYuK5Fm69JAWF%2Fuploads%2Fgit-blob-747a4b1e42111fb4ffb7e2fa8cc45e700d1693cf%2Foutput_12_0.png?alt=media)

Here, several distinct hotspots stand out clearly against the surrounding land: individual stacks and process units across the steelworks, each running well above ambient. [The facility plan](https://www.tatasteelnederland.com/sites/default/files/tata-steel-plattegrond.pdf) reveals that the DSP is located in the middle-left part of the cropped region. Since this image was taken before the DSP shutdown, its thermal signature is visible (the two hotspots near the shoreline). Let's take a look at how the overall thermal picture changes after the DSP was taken offline.

```python
# Close all opened raster files to free up resources
lwir1.close(); lwir2.close(); red.close(); green.close(); blue.close(); preview.close()
```

### Time Series: LWIR1 Before and After the DSP Shutdown

We have two LWIR1 acquisitions bracketing the shutdown: `20260302` (before, March 2026) and `20260619` (after, June 2026), both from the VZ02 sensor. We'll plot them side by side with a shared color scale for direct comparison.

```python
# Imagery scene IDs
TIMESERIES_SCENES = [
    "VZ02_L1b_20260302_132735",
    "VZ02_L1b_20260619_132916",
]


# Helper function to extract the date from a scene ID
def scene_date(scene_id):
    return scene_id.split("_")[2]


# Load and crop LWIR1 imagery for a given scene ID
def load_cropped_lwir1(scene_id):
    da = rioxarray.open_rasterio(
        Path("data") / scene_id / f"LWIR1.tiff", masked=True
    ).squeeze("band", drop=True)
    da = da.rio.clip_box(*ROI_BOUNDS, crs="EPSG:4326") * BT_SCALE
    return kelvin_to_fahrenheit(da)


# Load the LWIR1 time series for all scenes
timeseries = {scene_id: load_cropped_lwir1(scene_id) for scene_id in TIMESERIES_SCENES}
dates = sorted(timeseries, key=scene_date)
```

```python
from datetime import datetime

# Determine the min and max brightness temperature values for consistent color scaling
vmin = min(float(da.min()) for da in timeseries.values())
vmax = max(float(da.max()) for da in timeseries.values())

# Plot the LWIR1 time series
fig, axes = plt.subplots(1, len(dates), figsize=(7 * len(dates), 6))

for ax, scene_id in zip(axes, dates):
    da = timeseries[scene_id]
    da.plot.imshow(ax=ax, cmap="magma", vmin=vmin, vmax=vmax, cbar_kwargs={"label": "°F"})
    ax.set_title(datetime.strptime(scene_date(scene_id), "%Y%m%d").strftime("%B %d, %Y"))
    ax.axis("off")

fig.suptitle("LWIR1 Brightness Temperature")
plt.tight_layout()
plt.show()
```

![](https://3702723385-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FBgOfOR7RYuK5Fm69JAWF%2Fuploads%2Fgit-blob-253c9d8831a9747a898857d5710d081fda684ff9%2Foutput_17_0.png?alt=media)

When sharing a color scale, the March scene shows a handful of sharp, distinct hotspots against a cool background, while the June scene is warm across almost the entire site. This is a product of summer solar heating, making it hard to tell whether the DSP (or anything else) is active. Normalizing each frame relative to its own background can help separate genuine activity from ambient seasonal heating.

### Normalizing the Time Series

Different acquisitions can have different baseline brightness temperatures due to atmospheric and seasonal conditions unrelated to plant activity. Simply subtracting each frame's median still leaves frames with a lot of background texture (e.g. a hot summer day) looking uniformly "hot," since the deviations are on a different scale from frame to frame.

Instead, we use a **robust z-score**, which normalizes both the center and the spread of each frame using statistics that resist being skewed by a small number of extreme hot pixels:

$$z = \frac{x - \text{median}(x)}{1.4826 \times \text{MAD}(x)}$$

where `MAD(x) = median(|x - median(x)|)` is the median absolute deviation. Subtracting the median (instead of the mean) keeps the center from being pulled up by a handful of very hot pixels, and dividing by the (rescaled) MAD instead of the standard deviation keeps the scale from being inflated by those same outliers. The `1.4826` constant rescales MAD so it's comparable to a standard deviation for normally-distributed data, making the resulting z-score readable the same way as an ordinary one (e.g. `+3` ≈ "3 robust standard deviations above typical").

```python
# Compute robust z-score normalized timeseries for each scene
def robust_zscore(da):
    median = float(da.median())
    mad = float(abs(da - median).median()) * 1.4826  # scaled to be comparable to std
    return (da - median) / mad


# Apply robust z-score normalization to all scenes in the timeseries
normalized_timeseries = {
    scene_id: robust_zscore(da) for scene_id, da in timeseries.items()
}
```

```python
# Determine the global min and max for the normalized timeseries for consistent color scaling
norm_vmin = min(float(da.min()) for da in normalized_timeseries.values())
norm_vmax = max(float(da.max()) for da in normalized_timeseries.values())


# Plot the normalized timeseries
fig, axes = plt.subplots(1, len(dates), figsize=(7 * len(dates), 6))

for ax, scene_id in zip(axes, dates):
    da = normalized_timeseries[scene_id]
    da.plot.imshow(
        ax=ax,
        cmap="magma",
        vmin=norm_vmin,
        vmax=norm_vmax,
        cbar_kwargs={"label": "Robust z-score"},
    )
    ax.set_title(datetime.strptime(scene_date(scene_id), "%Y%m%d").strftime("%B %d, %Y"))
    ax.axis("off")

fig.suptitle("LWIR1 Brightness Temperature (robust z-score)")
plt.tight_layout()
plt.show()
```

![](https://3702723385-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FBgOfOR7RYuK5Fm69JAWF%2Fuploads%2Fgit-blob-457e552398f6ee0282db4cc923219ad7d322efc7%2Foutput_21_0.png?alt=media)

After normalization, the difference between the images becomes clear: by June, Tata Steel has shut off the DSP (as per the article), and the Warmband facility (the hot-rolled strip steel mill) also appears to be offline. With a few simple image processing steps, we have gone from raw imagery to actionable compliance insights.

### Thermal Profile Through the Hottest Point

Beyond a visual inspection, extracting pixel values can provide further insight into the thermal activity of a scene. A profile (a cross-sectional view that shows the change in pixel values along a linear path) through the hottest point tells us two distinct things at once:

1. Whether a given unit is active or idle (does the profile spike well above the surrounding background, or does it stay flat and close to ambient?)
2. How active it is (values along the profile indicate the intensity of the thermal signal, not just its presence)

We'll use the `20260302` acquisition, taken before the DSP shutdown, since it's where the plant's activity registers most clearly.

```python
from rasterio.warp import transform

# Get the data for the profile scene
profile_scene_id = "VZ02_L1b_20260302_132735"
profile_da = timeseries[profile_scene_id]

# Find the hotspot location in the profile scene
hotspot_idx = profile_da.argmax(dim=["y", "x"])
hotspot_row = int(hotspot_idx["y"])
hotspot_col = int(hotspot_idx["x"])
hotspot_northing = float(profile_da["y"].isel(y=hotspot_row))
hotspot_easting = float(profile_da["x"].isel(x=hotspot_col))

# Extract the transect along the hotspot row
transect = profile_da.isel(y=hotspot_row)

# Convert the hotspot's own coordinates to latitude for the title
_, hotspot_lat = transform(
    profile_da.rio.crs, "EPSG:4326", [hotspot_easting], [hotspot_northing]
)
hotspot_lat = hotspot_lat[0]

# Convert the transect's easting coordinates to longitude for the x-axis
transect_lons, _ = transform(
    profile_da.rio.crs,
    "EPSG:4326",
    transect["x"].values,
    [hotspot_northing] * transect.sizes["x"],
)

# Compute overall statistics for the profile scene
overall_min = float(profile_da.min())
overall_max = float(profile_da.max())
overall_mean = float(profile_da.mean())
overall_std = float(profile_da.std())

# Plot the profile scene and the transect through the hotspot
fig, (ax_img, ax_profile) = plt.subplots(1, 2, figsize=(16, 6))

profile_da.plot.imshow(ax=ax_img, cmap="magma", cbar_kwargs={"label": "°F"})
ax_img.axhline(hotspot_northing, color="cyan", linestyle="--", linewidth=1.5)
ax_img.set_title(f"{profile_scene_id} - Latitude: {hotspot_lat:.4f}°N")
ax_img.axis("off")

ax_profile.plot(transect_lons, transect.values, color="crimson", label="Transect")
ax_profile.axhline(overall_mean, color="black", linewidth=1, label="Mean")
ax_profile.axhspan(
    overall_mean - overall_std,
    overall_mean + overall_std,
    color="gray",
    alpha=0.3,
    label="Mean ± 1 std",
)
ax_profile.axhline(overall_min, color="steelblue", linestyle=":", label="Min")
ax_profile.axhline(overall_max, color="firebrick", linestyle=":", label="Max")
ax_profile.set_xlabel("Longitude")
ax_profile.set_ylabel("°F")
ax_profile.set_title("Thermal Profile")
ax_profile.legend()

plt.tight_layout()
plt.show()
```

![](https://3702723385-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FBgOfOR7RYuK5Fm69JAWF%2Fuploads%2Fgit-blob-da8e4d68724e73f197c9b76da5a03c08ea2fe081%2Foutput_24_0.png?alt=media)

This profile is (taken through the Warmband) shows a clear signal of activity. By extracting similar profiles through our images, or by taking pixel statistics for an AOI polygon, we can further learn about the thermal characteristics of locations and derive insights from thermal imagery.