如何用Python（非arcpy）将大栅格按矢量裁剪为6个NumPy数组？

Got it, let's work through this solution together—you can absolutely do this with open-source Python libraries, no ArcPy required. We'll use rasterio for raster handling, geopandas for polygon processing, and numpy for the final array outputs. Here's a step-by-step breakdown:

Step 1: Install Required Libraries #

First, make sure you have the necessary packages installed. Run this in your terminal:

pip install rasterio geopandas numpy

Step 2: Load Your Data #

We'll load both the large raster and your polygon data (either a Shapefile or an existing GeoPandas DataFrame):

import rasterio
from rasterio.mask import mask
import geopandas as gpd
import numpy as np

# Load the raster file
raster_path = "path/to/your/large_raster.tif"
with rasterio.open(raster_path) as src:
    raster_meta = src.meta  # Store raster metadata (CRS, dimensions, etc.)
    raster_crs = src.crs    # Get raster coordinate system

# Load polygons - skip this block if you already have a GeoDataFrame
polygon_path = "path/to/your/polygons.shp"
polygons = gpd.read_file(polygon_path)

Step 3: Align Coordinate Systems #

This is critical—if your raster and polygons don't share the same CRS, the crop won't work correctly. Let's fix that:

if polygons.crs != raster_crs:
    # Convert polygons to match the raster's CRS
    polygons = polygons.to_crs(raster_crs)

Step 4: Crop Raster for Each Polygon #

Now we'll loop through each polygon, crop the raster to its extent, and store the result as a NumPy array:

cropped_arrays = []

# Iterate over each polygon in the GeoDataFrame
for idx, polygon_row in polygons.iterrows():
    # Convert polygon geometry to a GeoJSON-like format (required by rasterio)
    polygon_geom = [polygon_row.geometry.__geo_interface__]
    
    # Crop the raster to the polygon
    # `crop=True` ensures we only get the area inside the polygon
    # `nodata=np.nan` replaces pixels outside the polygon with NaN (adjust as needed)
    cropped_raster, _ = mask(src=src, shapes=polygon_geom, crop=True, nodata=np.nan)
    
    # Remove the extra band dimension (since we're working with a single band raster)
    # If you have multi-band data, skip this step or adjust accordingly
    cropped_array = np.squeeze(cropped_raster)
    
    # Add the cropped array to our list
    cropped_arrays.append(cropped_array)

Step 5: Verify the Result #

After running the code, cropped_arrays will contain 6 NumPy arrays (one for each polygon). You can check their shapes to confirm:

for i, arr in enumerate(cropped_arrays):
    print(f"Polygon {i+1} array shape: {arr.shape}")

Key Notes to Adjust for Your Use Case #

• Multi-band rasters: If your raster has multiple bands, src.read() will return a (num_bands, height, width) array. The mask function will preserve this structure—just skip the np.squeeze() step if you want to keep all bands.
• Pixel inclusion: Use the all_touched parameter in mask() to control whether pixels partially touching the polygon are included. Set all_touched=True to include them, or False to only include pixels fully inside the polygon.
• Nodata values: Replace np.nan with your raster's actual nodata value (e.g., -9999) if needed.
• Large rasters: For extremely large rasters that don't fit in memory, consider using rasterio's windowed reading to process chunks at a time.

内容的提问来源于stack exchange，提问作者Lisa Mathew