Data Handling

1. Writing Rasters¶

In the previous chapter, we used map algebra to calculate a binary snow mask from satellite imagery. However, that result currently exists only as a temporary NumPy array in your computer’s memory. To use this snow mask in geographic information systems (GIS) software or share it with colleagues, we must write it back to the hard drive as a spatial GeoTIFF.

Writing a standard image is simple, but a GeoTIFF requires strict spatial metadata. If we just save the NumPy array, it loses all geographic context and the output file will not know where it belongs on Earth. The golden rule of raster export is to always copy the .profile (the complete metadata dictionary) from your original input raster. This guarantees that your new output raster will perfectly align in space with the original dataset.

When you perform map algebra, you often change the fundamental structure of the data. For example, our input was a 7 band image of 16 bit integers, but our newly calculated snow mask is a single band of binary values (True or False). Exporting a simple binary mask as a 32 bit float is bad practice and wastes massive amounts of disk space. Instead, we should convert the boolean mask to an 8 bit unsigned integer (uint8) so it only stores zeros and ones.

Concept Check: The Bloated Export¶

Scenario: You have calculated a binary flood mask (values 0 and 1) from a 32-bit continuous elevation raster. You copy the original profile dictionary and instantly write the mask to disk using dst.write(mask, 1). Why might your output file be unnecessarily massive in size?

A) Binary masks inherently take up more disk space than continuous data because they cannot be compressed.

B) You did not convert the mask to a GeoDataFrame before exporting.

C) You forgot to update the profile’s dtype to uint8, meaning the file is allocating 32 bits of memory for every single 0 and 1.

Before writing the file, we must update the copied profile dictionary to reflect these new data types. Additionally, to keep our project workspace organized, it is best practice to save all generated files into a dedicated folder rather than cluttering the main directory. We can use Python’s built-in os module to automatically create an outputs directory if it does not already exist.

Once the profile is updated and our output directory is ready, we open a new file connection using 'w' (write) mode, pass in the profile arguments, and use .write() to save the array to the hard drive.

2. Handling Large Datasets¶

Satellite imagery can be massive. A single full scene might be several gigabytes in size. Loading an entire file into memory with a simple .read() can easily overload your RAM and crash your computer. When working with large files, you must actively manage memory by processing the data efficiently.

• Read one band: If you only need the Near Infrared band for your analysis, avoid loading the entire multispectral stack. Explicitly use src.read(4) instead of src.read().

• Use lower precision: Arrays consume memory based on their data type. A float64 array takes up twice as much RAM as a float32 array. Use .astype("float32") to downgrade precision if extreme decimal accuracy is not required.

• Read one window: You can load a small, specific bounding box of pixels into memory without touching the rest of the file using Rasterio’s Window class.

• Process chunk by chunk: The most robust way to handle massive files is to loop through the raster in blocks. You open the file, read a small chunk into memory, process it, write that chunk directly to the new output file, and then move on to the next block.

Once this subset is securely in memory, it behaves exactly like any other NumPy array. Before applying complex algorithms or map algebra to large datasets, it is highly recommended to visualize the extracted chunk to verify you grabbed the correct geographic area.

While manually defining a single window is perfect for targeted extraction, applying map algebra to an entire multi-gigabyte scene requires automation. You can bridge this gap and process a massive file sequentially using Rasterio’s built-in block_windows() iterator.

3. Clipping¶

Often, your study area is much smaller than the satellite scene you downloaded. Clipping restricts your analysis to a specific region, vastly reducing processing time and memory usage.

Slice the Array¶

Because rasters are essentially NumPy arrays once read into memory, the fastest way to clip data is to simply slice the rows and columns using index numbers.

Use Geographic Bounds¶

Slicing by pixel index is fast, but it is not geographically accurate. If you know the exact real-world coordinate bounding box of your study area, you can use Rasterio to safely convert those coordinates into a pixel window.

Clip with a Shapefile¶

The most common and precise clipping method uses vector boundaries (polygons) to extract exact environmental features. In this example, we will clip our raster to the exact outline of the Haupapa/Tasman glacier in New Zealand (a glacier similar in size to the Aletsch glacier in Switzerland).

Rasterio provides a dedicated mask module to handle this complex geometry intersection.

4. Reprojection¶

Different regions and datasets use different Coordinate Reference Systems (CRS) to minimize geographic distortion. For example, spatial data in Switzerland is commonly projected in local coordinates (LV95, EPSG:2056) or global zones (UTM32N, EPSG:32632).

Our current New Zealand dataset is in UTM Zone 59S (EPSG:32759). However, local environmental agencies primarily use the New Zealand Transverse Mercator (NZTM, EPSG:2193). To combine our snow mask with local government vector data, we must reproject our raster.

Reprojecting a raster is computationally intensive because the computer cannot simply stretch the image. It must recalculate the size, shape, and value of every individual pixel in a new mathematical grid—a process known as resampling. Rasterio handles this via its warp module, which acts as a Python wrapper around the powerful C++ library GDAL (Geospatial Data Abstraction Library). GDAL is the industry standard engine that powers almost all open-source geographic software (including QGIS).

Because we are warping the grid, we must first calculate the new dimensions and transformation matrix for the destination CRS before we can actually move the pixel values.

5. Exercise: Extract and Export¶

Instead of reprojecting massive raster files, a standard spatial optimization is to reproject your lightweight vector boundaries to match the raster, clip the data, and then export the much smaller result.

Imagine a colleague needs a standalone GeoTIFF of the Near Infrared data (Band 4), restricted exactly to the boundaries of Lake Mapourika.

Tasks:

1. Open the Landsat subset (LC09_L1TP_075090_20230224_20230308_02_T1_subset.tif) and load the Lake Mapourika vector (JRC_lake_mapourika.gpkg).

2. Reproject the vector dataframe to match the raster’s CRS.

3. Use rasterio.mask.mask to clip the raster using the lake’s geometry.

4. Extract just the NIR band from the clipped multiband output.

5. Update the original spatial profile with the new height, width, and transform returned by the mask function. Remember to set count=1 since we are only saving the NIR band.

6. Write the result to a new file named outputs/lake_mapourika_nir.tif.

7. Plot the nir_clipped array to verify your extraction. Hint: The mask function sets pixels outside the polygon to 0. Use np.where() to convert these zeros to NaN (NaN, Not a Number) so the background renders cleanly.

Sample Solution

Notice how much faster and cleaner this script is compared to reprojecting the entire raster grid. We let GeoPandas handle the reprojection in a single line, use Rasterio to mask the data, and immediately write the small subset to disk. We also use the optimized lake_gdf.geometry.tolist() method to extract our shapes.

import os
import rasterio
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
from rasterio.mask import mask

input_raster = "data/LC09_L1TP_075090_20230224_20230308_02_T1_subset.tif"
vector_path = "data/JRC_lake_mapourika.gpkg"
output_raster = "outputs/lake_mapourika_nir.tif"

os.makedirs("outputs", exist_ok=True)

# 1. Load vector
lake_gdf = gpd.read_file(vector_path)

with rasterio.open(input_raster) as src:
    # 2. Reproject the vector to match the raster
    lake_gdf = lake_gdf.to_crs(src.crs)
    shapes = lake_gdf.geometry.tolist()
    
    # 3. Clip the raster
    out_image, out_transform = mask(src, shapes, crop=True)
    
    # 4. Extract Band 4 (Index 3 in the out_image array)
    nir_clipped = out_image[3]
    
    # 5. Copy and update the profile
    profile = src.profile.copy()
    profile.update({
        "height": nir_clipped.shape[0],
        "width": nir_clipped.shape[1],
        "transform": out_transform,
        "count": 1
    })

# 6. Write to disk
with rasterio.open(output_raster, "w", **profile) as dst:
    dst.write(nir_clipped, 1)

print(f"Successfully extracted and saved: {output_raster}")

# 7. Visual Confirmation
plt.figure(figsize=(6, 6))

# Replace 0 (the masked background) with NaN for a clean white background
nir_clean = np.where(nir_clipped == 0, np.nan, nir_clipped)

plt.imshow(nir_clean, cmap="viridis")
plt.title("Extracted Lake Mapourika (NIR)", fontsize=15)
plt.axis("off")
plt.show()

6. Summary¶

Effectively managing how data moves in and out of your computer’s memory is crucial for functional spatial programming.

• Writing requires metadata: Arrays lack spatial context. Always copy the .profile from the source image, update the data types if your math altered them, and pass those parameters into the write function.

• Respect memory limits: Satellite files are massive. Protect your system by reading only the specific bands you need, reducing data precision, or utilizing windowed reading.

• Clip to save time: Subsetting arrays using array indices, geographic coordinates, or exact vector polygons (masks) restricts processing strictly to your area of interest.

• Reproject to align: Use rasterio.warp to resample and transform pixel grids to match local Coordinate Reference Systems (like moving from UTM to local standard projections).